New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
nemogcm_tam.F90 in NEMO/branches/NERC/dev_release-3.4_NEMOTAM_consolidated/NEMOGCM/NEMO/OPATAM_SRC – NEMO

source: NEMO/branches/NERC/dev_release-3.4_NEMOTAM_consolidated/NEMOGCM/NEMO/OPATAM_SRC/nemogcm_tam.F90 @ 15454

Last change on this file since 15454 was 15454, checked in by smueller, 8 months ago

Addition of tangent-linear and adjoint model time-stepping loops and of a supporting structure for the optional inclusion of external NEMOTAM-related source code

File size: 35.8 KB
Line 
1MODULE nemogcm_tam
2#if defined key_tam
3   !!======================================================================
4   !!                       ***  MODULE nemogcm   ***
5   !! Ocean system   : NEMO GCM (ocean dynamics, on-line tracers, biochemistry and sea-ice)
6   !!======================================================================
7   !! History :  OPA  ! 1990-10  (C. Levy, G. Madec)  Original code
8   !!            7.0  ! 1991-11  (M. Imbard, C. Levy, G. Madec)
9   !!            7.1  ! 1993-03  (M. Imbard, C. Levy, G. Madec, O. Marti, M. Guyon, A. Lazar,
10   !!                             P. Delecluse, C. Perigaud, G. Caniaux, B. Colot, C. Maes) release 7.1
11   !!             -   ! 1992-06  (L.Terray)  coupling implementation
12   !!             -   ! 1993-11  (M.A. Filiberti) IGLOO sea-ice
13   !!            8.0  ! 1996-03  (M. Imbard, C. Levy, G. Madec, O. Marti, M. Guyon, A. Lazar,
14   !!                             P. Delecluse, L.Terray, M.A. Filiberti, J. Vialar, A.M. Treguier, M. Levy) release 8.0
15   !!            8.1  ! 1997-06  (M. Imbard, G. Madec)
16   !!            8.2  ! 1999-11  (M. Imbard, H. Goosse)  LIM sea-ice model
17   !!                 ! 1999-12  (V. Thierry, A-M. Treguier, M. Imbard, M-A. Foujols)  OPEN-MP
18   !!                 ! 2000-07  (J-M Molines, M. Imbard)  Open Boundary Conditions  (CLIPPER)
19   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and modules
20   !!             -   ! 2004-06  (R. Redler, NEC CCRLE, Germany) add OASIS[3/4] coupled interfaces
21   !!             -   ! 2004-08  (C. Talandier) New trends organization
22   !!             -   ! 2005-06  (C. Ethe) Add the 1D configuration possibility
23   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
24   !!             -   ! 2006-03  (L. Debreu, C. Mazauric)  Agrif implementation
25   !!             -   ! 2006-04  (G. Madec, R. Benshila)  Step reorganization
26   !!             -   ! 2007-07  (J. Chanut, A. Sellar) Unstructured open boundaries (BDY)
27   !!            3.2  ! 2009-08  (S. Masson)  open/write in the listing file in mpp
28   !!            3.3  ! 2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
29   !!             -   ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
30   !!            3.3.1! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation
31   !!            3.4  ! 2011-11  (C. Harris) decomposition changes for running with CICE
32   !! History of TAM:
33   !!            3.4  ! 2012-07  (P.-A. Bouttier) Phasing with 3.4 version
34   !!----------------------------------------------------------------------
35
36   !!----------------------------------------------------------------------
37   !!   nemo_gcm       : solve ocean dynamics, tracer, biogeochemistry and/or sea-ice
38   !!   nemo_init      : initialization of the NEMO system
39   !!   nemo_ctl       : initialisation of the contol print
40   !!   nemo_alloc     : dynamical allocation
41   !!   nemo_partition : calculate MPP domain decomposition
42   !!   factorise      : calculate the factors of the no. of MPI processes
43   !!----------------------------------------------------------------------
44   USE step_oce        ! module used in the ocean time stepping module
45   USE sbc_oce         ! surface boundary condition: ocean
46   USE cla             ! cross land advection               (tra_cla routine)
47   USE domcfg          ! domain configuration               (dom_cfg routine)
48   USE mppini          ! shared/distributed memory setting (mpp_init routine)
49   USE domain          ! domain initialization             (dom_init routine)
50   USE obcini          ! open boundary cond. initialization (obc_ini routine)
51   USE bdyini          ! open boundary cond. initialization (bdy_init routine)
52   USE bdydta          ! open boundary cond. initialization (bdy_dta_init routine)
53   USE bdytides        ! open boundary cond. initialization (tide_init routine)
54   USE istate          ! initial state setting          (istate_init routine)
55   USE ldfdyn          ! lateral viscosity setting      (ldfdyn_init routine)
56   USE ldftra          ! lateral diffusivity setting    (ldftra_init routine)
57   USE zdfini          ! vertical physics setting          (zdf_init routine)
58   USE phycst          ! physical constant                  (par_cst routine)
59   USE trdmod          ! momentum/tracers trends       (trd_mod_init routine)
60   USE diaptr          ! poleward transports           (dia_ptr_init routine)
61   USE diadct          ! sections transports           (dia_dct_init routine)
62   USE diaobs          ! Observation diagnostics       (dia_obs_init routine)
63   USE step            ! NEMO time-stepping                 (stp     routine)
64   USE tradmp
65   USE trabbl
66#if defined key_oasis3
67   USE cpl_oasis3      ! OASIS3 coupling
68#elif defined key_oasis4
69   USE cpl_oasis4      ! OASIS4 coupling (not working)
70#endif
71   USE c1d             ! 1D configuration
72   USE step_c1d        ! Time stepping loop for the 1D configuration
73#if defined key_top
74   USE trcini          ! passive tracer initialisation
75#endif
76   USE lib_mpp         ! distributed memory computing
77#if defined key_iomput
78   USE mod_ioclient
79#endif
80   USE nemogcm
81   USE step_tam
82   USE sbcssr_tam
83   USE step_oce_tam
84   USE zdf_oce_tam
85   USE trabbl_tam
86   USE tamtst
87   USE tamctl
88   USE lib_mpp_tam
89   USE paresp
90   !USE tamtrj
91   USE trj_tam
92   USE iom
93   USE wrk_nemo
94   IMPLICIT NONE
95   PRIVATE
96
97   PUBLIC   nemo_gcm_tam    ! called by model.F90
98   PUBLIC   nemo_init_tam   ! needed by AGRIF
99
100   CHARACTER(lc) ::   cform_aaa="( /, 'AAAAAAAA', / ) "     ! flag for output listing
101
102   !!----------------------------------------------------------------------
103   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
104   !! $Id$
105   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
106   !!----------------------------------------------------------------------
107CONTAINS
108
109   SUBROUTINE nemo_gcm_tam
110      !!----------------------------------------------------------------------
111      !!                     ***  ROUTINE nemo_gcm_tam  ***
112      !!
113      !! ** Purpose :   NEMO solves the primitive equations on an orthogonal
114      !!              curvilinear mesh on the sphere.
115      !!
116      !! ** Method  : - model general initialization
117      !!              - launch the time-stepping (stp routine)
118      !!              - finalize the run by closing files and communications
119      !!
120      !! References : Madec, Delecluse, Imbard, and Levy, 1997:  internal report, IPSL.
121      !!              Madec, 2008, internal report, IPSL.
122      !!----------------------------------------------------------------------
123      !                            !-----------------------!
124      !                            !==  Initialisations  ==!
125      CALL nemo_init               !-----------------------!
126      CALL nemo_init_tam         
127      !
128      ! check that all process are still there... If some process have an error,
129      ! they will never enter in step and other processes will wait until the end of the cpu time!
130      IF( lk_mpp )   CALL mpp_max( nstop )
131
132      IF(lwp) WRITE(numout,cform_aaa)   ! Flag AAAAAAA
133
134      !                            !-----------------------!
135      !                            !==   time stepping   ==!
136      !                            !-----------------------!
137      !
138      ! NEMOTAM options
139      ! ---------------
140      !
141      ! nn_tam = 1: Tangent-linear model time stepping
142      ! nn_tam = 2: Adjoint model time stepping
143      ! nn_tam = 3: Call external code (see below)
144      ! nn_tam = 0: Call tam_tst
145      !
146      ! External NEMOTAM code
147      ! ---------------------
148      !
149      ! When key_tam_ext is defined, two include files are required (e.g., in
150      ! the MY_SRC directory): tam_ext_subr.h90 and tam_ext_var.h90. All of the
151      ! subroutines listed below have to be provided in include file
152      ! tam_ext_subr.h90
153      !
154      !  +---------------------+--------------------+----------+------------------------------------+
155      !  | Subroutine          | called if nn_tam = | argument | position of call                   |
156      !  +---------------------+--------------------+----------+------------------------------------+
157      !  | tam_ext_tan_init    |                  1 |          | precedes istate_init_tan call      |
158      !  | tam_ext_tan_prestp  |                  1 |     istp | precedes stp_tan         call      |
159      !  | tam_ext_tan_poststp |                  1 |     istp | follows  stp_tan         call      |
160      !  | tam_ext_tan_final   |                  1 |          | follows  final time step           |
161      !  | tam_ext_adj_init    |                  2 |          | precedes first time step           |
162      !  | tam_ext_adj_prestp  |                  2 |     istp | precedes stp_adj call              |
163      !  | tam_ext_adj_poststp |                  2 |     istp | follows stp_adj call               |
164      !  | tam_ext_adj_final   |                  2 |          | follows istate_init_adj call       |
165      !  | tam_ext_main        |                  3 |          | n/a                                |
166      !  +---------------------+--------------------+----------+------------------------------------+
167      !
168      ! and tam_ext_vars.h90 can be empty or used to declare addtional
169      ! TAM-specific variables (declared in module tamctl).
170      !
171      IF ( nn_tam == 1 ) THEN      ! TL mode
172
173         ! Initialize energy weights
174         CALL par_esp
175
176         ! Time-step loop
177         CALL tam_tsloop_tan
178
179         IF (lwp) THEN
180            WRITE(numout,*)
181            WRITE(numout,*) ' TAM: TL finished'
182            WRITE(numout,*) ' ----------------'
183            WRITE(numout,*)
184         ENDIF
185         CALL flush(numout)
186
187      ELSEIF ( nn_tam == 2 ) THEN  ! ADJ mode
188
189         ! Initialize energy weights
190         CALL par_esp
191
192         ! Time-step loop
193         CALL tam_tsloop_adj
194
195         IF (lwp) THEN
196            WRITE(numout,*)
197            WRITE(numout,*) ' TAM: ADJ finished'
198            WRITE(numout,*) ' -----------------'
199            WRITE(numout,*)
200         ENDIF
201         CALL flush(numout)
202
203      ELSEIF ( nn_tam == 3 ) THEN
204
205         ! Initialize energy weights
206         CALL par_esp
207
208         CALL tam_ext_main
209
210      ELSE                         ! Test mode
211
212         CALL tam_tst
213
214         IF (lwp) THEN
215            WRITE(numout,*)
216            WRITE(numout,*) ' tamtst: Finished testing operators'
217            WRITE(numout,*) ' ------'
218            WRITE(numout,*)
219         ENDIF
220
221      ENDIF
222      !                            !------------------------!
223      !                            !==  finalize the run  ==!
224      !                            !------------------------!
225      IF(lwp) WRITE(numout,cform_aaa)   ! Flag AAAAAAA
226      !
227      IF( nstop /= 0 .AND. lwp ) THEN   ! error print
228         WRITE(numout,cform_err)
229         WRITE(numout,*) nstop, ' error have been found'
230      ENDIF
231      !
232      IF( nn_timing == 1 )   CALL timing_finalize
233      !!
234      CALL nemo_closefile
235      IF( lk_mpp )   CALL mppstop       ! end mpp communications
236      !
237   END SUBROUTINE nemo_gcm_tam
238
239   SUBROUTINE tam_tsloop_tan(ld_input, ld_output)
240      !!----------------------------------------------------------------------
241      !!                     *** ROUTINE tam_tsloop_tan ***
242      !!
243      !! ** Purpose : Default NEMOTAM tangent-linear time-stepping loop
244      !!
245      !!----------------------------------------------------------------------
246
247      LOGICAL, INTENT(in), OPTIONAL            ::   ld_input, ld_output   ! optional arguments (default .TRUE.)
248      LOGICAL                                  ::   ll_input, ll_output   !    - ld_input/ll_input:   input  of initial state
249                                                                          !    - ld_output/ll_output: output of final state
250      INTEGER                                  ::   istp, inum            ! timestep index, file unit
251      REAL(KIND=wp), POINTER, DIMENSION(:,:)   ::   z2d
252      REAL(KIND=wp), POINTER, DIMENSION(:,:,:) ::   z3d
253
254      ! Optional arguments
255      ll_input  = .TRUE. ; IF ( PRESENT( ld_input  ) ) ll_input   = ld_input
256      ll_output = .TRUE. ; IF ( PRESENT( ld_output ) ) ll_output  = ld_output
257
258      ! Update trajectory
259      istp = nit000 - 1
260      CALL trj_rea( istp, 1, .true. )
261
262      ! Call various initialisation subroutines
263      CALL     oce_tam_init( 1 )
264      CALL sbc_oce_tam_init( 1 )
265      CALL trc_oce_tam_init( 1 )   
266      CALL sol_oce_tam_init( 1 )
267
268      ! Initialise tangent-linear state variables
269      un_tl  (:,:,:)        = 0.0_wp
270      vn_tl  (:,:,:)        = 0.0_wp
271      tsn_tl (:,:,:,jp_tem) = 0.0_wp
272      tsn_tl (:,:,:,jp_sal) = 0.0_wp
273      sshn_tl(  :,:)        = 0.0_wp
274
275      ! Load initial tangent-linear state
276      IF ( ll_input ) THEN
277         CALL wrk_alloc( jpi, jpj,      z2d )
278         CALL wrk_alloc( jpi, jpj, jpk, z3d )
279         CALL iom_open( 'tan_input.nc', inum, kiolib = jprstlib )
280         CALL iom_get( inum, jpdom_autoglo, 'tn_tl',   z3d )
281         tsn_tl(:,:,:,jp_tem) = z3d
282         CALL iom_get( inum, jpdom_autoglo, 'sn_tl',   z3d )
283         tsn_tl(:,:,:,jp_sal) = z3d
284         CALL iom_get( inum, jpdom_autoglo, 'un_tl',   z3d )
285         un_tl(:,:,:)         = z3d
286         CALL iom_get( inum, jpdom_autoglo, 'vn_tl',   z3d ) 
287         vn_tl(:,:,:)         = z3d
288         CALL iom_get( inum, jpdom_autoglo, 'sshn_tl', z2d )
289         sshn_tl(:,:)         = z2d
290         CALL iom_close( inum )
291         CALL wrk_dealloc( jpi, jpj,      z2d )
292         CALL wrk_dealloc( jpi, jpj, jpk, z3d )
293      END IF
294
295      CALL tam_ext_tan_init
296
297      ! Tangent-linear version of istate
298      CALL istate_init_tan
299
300      ! Call day_tam, as call to day_init commented out in istate_init_tan
301      CALL day_tam( nit000, 0 )
302
303      ! Time step
304      DO istp = nit000, nitend, 1
305
306         CALL tam_ext_tan_prestp( istp )
307
308         CALL stp_tan( istp )
309
310         CALL tam_ext_tan_poststp( istp )
311
312      END DO
313
314      CALL tam_ext_tan_final
315
316      ! Output final tangent-linear state
317      IF ( ll_output ) THEN
318         CALL iom_open( 'tan_output.nc', inum, ldwrt = .TRUE., kiolib = jprstlib )
319         CALL iom_rstput( istp, istp, inum, 'tn_tl',   tsn_tl(:,:,:,jp_tem) )
320         CALL iom_rstput( istp, istp, inum, 'sn_tl',   tsn_tl(:,:,:,jp_sal) )
321         CALL iom_rstput( istp, istp, inum, 'un_tl',   un_tl                )
322         CALL iom_rstput( istp, istp, inum, 'vn_tl',   vn_tl                )
323         CALL iom_rstput( istp, istp, inum, 'sshn_tl', sshn_tl              )
324         CALL iom_close( inum )
325      END IF
326
327   END SUBROUTINE tam_tsloop_tan
328
329   SUBROUTINE tam_tsloop_adj(ld_input, ld_output)
330      !!----------------------------------------------------------------------
331      !!                     *** ROUTINE tam_tsloop_adj ***
332      !!
333      !! ** Purpose : Default NEMOTAM adjoint time-stepping loop
334      !!
335      !!----------------------------------------------------------------------
336
337      LOGICAL, INTENT(in), OPTIONAL            ::   ld_input, ld_output   ! optional arguments (default .TRUE.)
338      LOGICAL                                  ::   ll_input, ll_output   !    - ld_input/ll_input:   input  of initial state
339                                                                          !    - ld_output/ll_output: output of final state
340      INTEGER                                  ::   istp, inum            ! timestep index, file unit
341      REAL(KIND=wp), POINTER, DIMENSION(:,:)   :: z2d
342      REAL(KIND=wp), POINTER, DIMENSION(:,:,:) :: z3d
343
344      ! Optional arguments
345      ll_input  = .TRUE. ; IF ( PRESENT( ld_input  ) ) ll_input   = ld_input
346      ll_output = .TRUE. ; IF ( PRESENT( ld_output ) ) ll_output  = ld_output
347
348      ! Update calendar (go to final time step one at a time: calendar
349      ! is initialised at nit000, and subroutine 'day_tam' only
350      ! increments/decrements one time step each call)
351      DO istp = nit000, nitend
352         CALL day_tam( istp, 0 )
353      END DO
354 
355      ! Update trajectory
356      CALL trj_rea( istp-1, -1, .true. )
357
358      ! Call various initialisation subroutines
359      CALL oce_tam_init(2)
360      CALL sbc_oce_tam_init(2)
361      CALL trc_oce_tam_init(2)
362      CALL sol_oce_tam_init(2)
363
364      ! Initialise adjoint state variables
365      un_ad(:,:,:)    = 0.0_wp
366      vn_ad(:,:,:)    = 0.0_wp
367      sshn_ad(:,:)    = 0.0_wp
368      tsn_ad(:,:,:,:) = 0.0_wp
369      ub_ad(:,:,:)    = 0.0_wp
370      vb_ad(:,:,:)    = 0.0_wp
371      sshb_ad(:,:)    = 0.0_wp
372      tsb_ad(:,:,:,:) = 0.0_wp
373
374      ! Load initial adjoint state and only permit values in interior domain
375      IF ( ll_input ) THEN
376         CALL wrk_alloc( jpi, jpj,      z2d )
377         CALL wrk_alloc( jpi, jpj, jpk, z3d )
378         CALL iom_open( 'adj_input.nc', inum, kiolib = jprstlib )
379         CALL iom_get( inum, jpdom_autoglo, 'tn_ad',   z3d )
380         tsn_ad(:,:,:,jp_tem) = z3d * tmsk_i
381         CALL iom_get( inum, jpdom_autoglo, 'sn_ad',   z3d )
382         tsn_ad(:,:,:,jp_sal) = z3d * tmsk_i
383         CALL iom_get( inum, jpdom_autoglo, 'un_ad',   z3d )
384         un_ad(:,:,:)         = z3d * umsk_i
385         CALL iom_get( inum, jpdom_autoglo, 'vn_ad',   z3d ) 
386         vn_ad(:,:,:)         = z3d * vmsk_i
387         CALL iom_get( inum, jpdom_autoglo, 'sshn_ad', z2d )
388         sshn_ad(:,:)         = z2d * tmsk_i(:,:,1)
389         CALL iom_close( inum )
390         CALL wrk_dealloc( jpi, jpj,      z2d )
391         CALL wrk_dealloc( jpi, jpj, jpk, z3d )
392      END IF
393
394      CALL tam_ext_adj_init
395
396      ! Adjoint time steps (backwards in time)
397      DO istp = nitend, nit000, -1
398
399         CALL tam_ext_adj_prestp( istp )
400
401         ! Carry out adjoint time step
402         CALL stp_adj( istp )
403           
404         CALL tam_ext_adj_poststp( istp )
405
406      END DO
407
408      ! Adjoint version of istate
409      CALL istate_init_adj
410
411      ! Output final adjoint state
412      IF ( ll_output ) THEN
413         CALL iom_open( 'adj_output.nc', inum, ldwrt = .TRUE., kiolib = jprstlib )
414         CALL iom_rstput( istp, istp, inum, 'tn_ad',   tsn_ad(:,:,:,jp_tem) )
415         CALL iom_rstput( istp, istp, inum, 'sn_ad',   tsn_ad(:,:,:,jp_sal) )
416         CALL iom_rstput( istp, istp, inum, 'un_ad',   un_ad                )
417         CALL iom_rstput( istp, istp, inum, 'vn_ad',   vn_ad                )
418         CALL iom_rstput( istp, istp, inum, 'sshn_ad', sshn_ad              )
419         CALL iom_close( inum )
420      END IF
421
422      CALL tam_ext_adj_final
423
424   END SUBROUTINE tam_tsloop_adj
425
426#if defined key_tam_ext
427#include "tam_ext_subr.h90"
428#else
429   SUBROUTINE tam_ext_tan_init
430   END SUBROUTINE tam_ext_tan_init
431   SUBROUTINE tam_ext_tan_prestp( kstp )
432      INTEGER :: kstp
433   END SUBROUTINE tam_ext_tan_prestp
434   SUBROUTINE tam_ext_tan_poststp( kstp )
435      INTEGER :: kstp
436   END SUBROUTINE tam_ext_tan_poststp
437   SUBROUTINE tam_ext_tan_final
438   END SUBROUTINE tam_ext_tan_final
439   SUBROUTINE tam_ext_adj_init
440   END SUBROUTINE tam_ext_adj_init
441   SUBROUTINE tam_ext_adj_prestp( kstp )
442      INTEGER :: kstp
443   END SUBROUTINE tam_ext_adj_prestp
444   SUBROUTINE tam_ext_adj_poststp( kstp )
445      INTEGER :: kstp
446   END SUBROUTINE tam_ext_adj_poststp
447   SUBROUTINE tam_ext_adj_final
448   END SUBROUTINE tam_ext_adj_final
449   SUBROUTINE tam_ext_main
450   END SUBROUTINE tam_ext_main
451#endif
452
453   SUBROUTINE nemo_init_tam
454      !!----------------------------------------------------------------------
455      !!                     ***  ROUTINE nemo_init  ***
456      !!
457      !! ** Purpose :   initialization of the NEMO GCM
458      !!----------------------------------------------------------------------
459      INTEGER ::   ji            ! dummy loop indices
460      INTEGER ::   ilocal_comm   ! local integer
461      CHARACTER(len=80), DIMENSION(16) ::   cltxt
462      !!
463      NAMELIST/namctl/ ln_ctl  , nn_print, nn_ictls, nn_ictle,   &
464         &             nn_isplt, nn_jsplt, nn_jctls, nn_jctle,   &
465         &             nn_bench, nn_timing
466      NAMELIST/namtam/ nn_tam
467      !!----------------------------------------------------------------------
468      !
469      ! Read 'namtam' namelist
470      nn_tam = 1                  ! Default mode (tangent-linear)
471      REWIND( numnam )
472      READ  ( numnam, namtam )
473      IF (lwp) THEN
474         WRITE(numout,*)
475         WRITE(numout,*) 'tam : NEMOTAM mode'
476         WRITE(numout,*) '~~~'
477         WRITE(numout,*) '   Namelist namtam'
478         WRITE(numout,*) '      tam mode (0 - TEST, 1 - TL (default), 2 - ADJ, 3 - external)  nn_tam = ', nn_tam
479      END IF
480      !
481      IF( ln_rnf        )   CALL sbc_rnf_init
482      !!!!!!!!!!!!! TAM initialisation !!!!!!!!!!!!!!!!!!!!!!!!!!!
483      CALL nemo_alloc_tam
484      CALL nemo_ctl_tam                          ! Control prints & Benchmark
485
486                            CALL  istate_init_tan   ! ocean initial state (Dynamics and tracers)
487                            CALL  istate_init_adj   ! ocean initial state (Dynamics and tracers)
488      !                                     ! Ocean physics
489                            CALL     sbc_init_tam   ! Forcings : surface module
490                            CALL     sbc_ssr_ini_tam   ! Forcings : surface module
491            !                                     ! Active tracers
492                            CALL tra_qsr_init_tam   ! penetrative solar radiation qsr
493      IF( lk_trabbl     )   CALL tra_bbl_init_tam   ! advective (and/or diffusive) bottom boundary layer scheme
494      IF( ln_tradmp     )   CALL tra_dmp_init_tam   ! internal damping trends
495                            CALL tra_adv_init_tam   ! horizontal & vertical advection
496                            CALL tra_ldf_init_tam   ! lateral mixing
497                            CALL tra_zdf_init_tam   ! vertical mixing and after tracer fields
498
499      !                                     ! Dynamics
500                            CALL dyn_adv_init_tam   ! advection (vector or flux form)
501                            CALL dyn_vor_init_tam   ! vorticity term including Coriolis
502                            CALL dyn_ldf_init_tam   ! lateral mixing
503                            CALL dyn_hpg_init_tam   ! horizontal gradient of Hydrostatic pressure
504                            CALL dyn_zdf_init_tam   ! vertical diffusion
505                            CALL dyn_spg_init_tam   ! surface pressure gradient
506
507      !                                     ! Misc. options
508      IF( nn_cla == 1   )   CALL cla_init_tam       ! Cross Land Advection
509                            CALL sbc_rnf_init_tam
510
511      IF( nn_tam == 0 )     CALL tam_tst_init
512                            CALL tl_trj_ini
513   END SUBROUTINE nemo_init_tam
514
515   SUBROUTINE nemo_ctl_tam
516      !!----------------------------------------------------------------------
517      !!                     ***  ROUTINE nemo_ctl  ***
518      !!
519      !! ** Purpose :   control print setting
520      !!
521      !! ** Method  : - print namctl information and check some consistencies
522      !!----------------------------------------------------------------------
523      !
524      IF(lwp) THEN                  ! control print
525         WRITE(numout,*)
526         WRITE(numout,*) 'nemo_ctl_tam: Control prints & Benchmark'
527         WRITE(numout,*) '~~~~~~~ '
528         WRITE(numout,*) '   Namelist namctl'
529         WRITE(numout,*) '      run control (for debugging)     ln_ctl     = ', ln_ctl
530         WRITE(numout,*) '      level of print                  nn_print   = ', nn_print
531         WRITE(numout,*) '      Start i indice for SUM control  nn_ictls   = ', nn_ictls
532         WRITE(numout,*) '      End i indice for SUM control    nn_ictle   = ', nn_ictle
533         WRITE(numout,*) '      Start j indice for SUM control  nn_jctls   = ', nn_jctls
534         WRITE(numout,*) '      End j indice for SUM control    nn_jctle   = ', nn_jctle
535         WRITE(numout,*) '      number of proc. following i     nn_isplt   = ', nn_isplt
536         WRITE(numout,*) '      number of proc. following j     nn_jsplt   = ', nn_jsplt
537         WRITE(numout,*) '      benchmark parameter (0/1)       nn_bench   = ', nn_bench
538      ENDIF
539      !
540      nprint    = nn_print          ! convert DOCTOR namelist names into OLD names
541      nictls    = nn_ictls
542      nictle    = nn_ictle
543      njctls    = nn_jctls
544      njctle    = nn_jctle
545      isplt     = nn_isplt
546      jsplt     = nn_jsplt
547      nbench    = nn_bench
548      !                             ! Parameter control
549      !
550      IF( ln_ctl ) THEN                 ! sub-domain area indices for the control prints
551         IF( lk_mpp .AND. jpnij > 1 ) THEN
552            isplt = jpni   ;   jsplt = jpnj   ;   ijsplt = jpni*jpnj   ! the domain is forced to the real split domain
553         ELSE
554            IF( isplt == 1 .AND. jsplt == 1  ) THEN
555               CALL ctl_warn( ' - isplt & jsplt are equal to 1',   &
556                  &           ' - the print control will be done over the whole domain' )
557            ENDIF
558            ijsplt = isplt * jsplt            ! total number of processors ijsplt
559         ENDIF
560         IF(lwp) WRITE(numout,*)'          - The total number of processors over which the'
561         IF(lwp) WRITE(numout,*)'            print control will be done is ijsplt : ', ijsplt
562         !
563         !                              ! indices used for the SUM control
564         IF( nictls+nictle+njctls+njctle == 0 )   THEN    ! print control done over the default area
565            lsp_area = .FALSE.
566         ELSE                                             ! print control done over a specific  area
567            lsp_area = .TRUE.
568            IF( nictls < 1 .OR. nictls > jpiglo )   THEN
569               CALL ctl_warn( '          - nictls must be 1<=nictls>=jpiglo, it is forced to 1' )
570               nictls = 1
571            ENDIF
572            IF( nictle < 1 .OR. nictle > jpiglo )   THEN
573               CALL ctl_warn( '          - nictle must be 1<=nictle>=jpiglo, it is forced to jpiglo' )
574               nictle = jpiglo
575            ENDIF
576            IF( njctls < 1 .OR. njctls > jpjglo )   THEN
577               CALL ctl_warn( '          - njctls must be 1<=njctls>=jpjglo, it is forced to 1' )
578               njctls = 1
579            ENDIF
580            IF( njctle < 1 .OR. njctle > jpjglo )   THEN
581               CALL ctl_warn( '          - njctle must be 1<=njctle>=jpjglo, it is forced to jpjglo' )
582               njctle = jpjglo
583            ENDIF
584         ENDIF
585      ENDIF
586      !
587      IF( nbench == 1 ) THEN              ! Benchmark
588         SELECT CASE ( cp_cfg )
589         CASE ( 'gyre' )   ;   CALL ctl_warn( ' The Benchmark is activated ' )
590         CASE DEFAULT      ;   CALL ctl_stop( ' The Benchmark is based on the GYRE configuration:',   &
591            &                                 ' key_gyre must be used or set nbench = 0' )
592         END SELECT
593      ENDIF
594      !
595      IF( lk_c1d .AND. .NOT.lk_iomput )   CALL ctl_stop( 'nemo_ctl_tam: The 1D configuration must be used ',   &
596         &                                               'with the IOM Input/Output manager. '         ,   &
597         &                                               'Compile with key_iomput enabled' )
598      !
599   END SUBROUTINE nemo_ctl_tam
600
601   SUBROUTINE nemo_alloc_tam
602      !!----------------------------------------------------------------------
603      !!                     ***  ROUTINE nemo_alloc  ***
604      !!
605      !! ** Purpose :   Allocate all the dynamic arrays of the OPA modules
606      !!
607      !! ** Method  :
608      !!----------------------------------------------------------------------
609      !
610      INTEGER :: ierr
611      !!----------------------------------------------------------------------
612      !
613      ierr =        oce_alloc_tam       ( 0 )          ! ocean
614      ierr = ierr + zdf_oce_alloc_tam   (   )          ! ocean vertical physics
615      !
616      ierr = ierr + lib_mpp_alloc_adj   (numout)    ! mpp exchanges
617      ierr = ierr + trc_oce_alloc_tam   ( 0 )          ! shared TRC / TRA arrays
618      ierr = ierr + sbc_oce_alloc_tam   ( 0 )          ! shared TRC / TRA arrays
619      ierr = ierr + sol_oce_alloc_tam   ( 0 )          ! shared TRC / TRA arrays
620      !
621      IF( lk_mpp    )   CALL mpp_sum( ierr )
622      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'nemo_alloc_tam : unable to allocate standard ocean arrays' )
623      !
624   END SUBROUTINE nemo_alloc_tam
625
626   SUBROUTINE factorise( kfax, kmaxfax, knfax, kn, kerr )
627      !!----------------------------------------------------------------------
628      !!                     ***  ROUTINE factorise  ***
629      !!
630      !! ** Purpose :   return the prime factors of n.
631      !!                knfax factors are returned in array kfax which is of
632      !!                maximum dimension kmaxfax.
633      !! ** Method  :
634      !!----------------------------------------------------------------------
635      INTEGER                    , INTENT(in   ) ::   kn, kmaxfax
636      INTEGER                    , INTENT(  out) ::   kerr, knfax
637      INTEGER, DIMENSION(kmaxfax), INTENT(  out) ::   kfax
638      !
639      INTEGER :: ifac, jl, inu
640      INTEGER, PARAMETER :: ntest = 14
641      INTEGER :: ilfax(ntest)
642
643      ! lfax contains the set of allowed factors.
644      data (ilfax(jl),jl=1,ntest) / 16384, 8192, 4096, 2048, 1024, 512, 256,  &
645         &                            128,   64,   32,   16,    8,   4,   2  /
646      !!----------------------------------------------------------------------
647
648      ! Clear the error flag and initialise output vars
649      kerr = 0
650      kfax = 1
651      knfax = 0
652
653      ! Find the factors of n.
654      IF( kn == 1 )   GOTO 20
655
656      ! nu holds the unfactorised part of the number.
657      ! knfax holds the number of factors found.
658      ! l points to the allowed factor list.
659      ! ifac holds the current factor.
660
661      inu   = kn
662      knfax = 0
663
664      DO jl = ntest, 1, -1
665         !
666         ifac = ilfax(jl)
667         IF( ifac > inu )   CYCLE
668
669         ! Test whether the factor will divide.
670
671         IF( MOD(inu,ifac) == 0 ) THEN
672            !
673            knfax = knfax + 1            ! Add the factor to the list
674            IF( knfax > kmaxfax ) THEN
675               kerr = 6
676               write (*,*) 'FACTOR: insufficient space in factor array ', knfax
677               return
678            ENDIF
679            kfax(knfax) = ifac
680            ! Store the other factor that goes with this one
681            knfax = knfax + 1
682            kfax(knfax) = inu / ifac
683            !WRITE (*,*) 'ARPDBG, factors ',knfax-1,' & ',knfax,' are ', kfax(knfax-1),' and ',kfax(knfax)
684         ENDIF
685         !
686      END DO
687
688   20 CONTINUE      ! Label 20 is the exit point from the factor search loop.
689      !
690   END SUBROUTINE factorise
691
692#if defined key_mpp_mpi
693   SUBROUTINE nemo_northcomms
694      !!======================================================================
695      !!                     ***  ROUTINE  nemo_northcomms  ***
696      !! nemo_northcomms    :  Setup for north fold exchanges with explicit peer to peer messaging
697      !!=====================================================================
698      !!----------------------------------------------------------------------
699      !!
700      !! ** Purpose :   Initialization of the northern neighbours lists.
701      !!----------------------------------------------------------------------
702      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)
703      !!----------------------------------------------------------------------
704
705      INTEGER ::   ji, jj, jk, ij, jtyp    ! dummy loop indices
706      INTEGER ::   ijpj                    ! number of rows involved in north-fold exchange
707      INTEGER ::   northcomms_alloc        ! allocate return status
708      REAL(wp), ALLOCATABLE, DIMENSION ( :,: ) ::   znnbrs     ! workspace
709      LOGICAL,  ALLOCATABLE, DIMENSION ( : )   ::   lrankset   ! workspace
710
711      IF(lwp) WRITE(numout,*)
712      IF(lwp) WRITE(numout,*) 'nemo_northcomms : Initialization of the northern neighbours lists'
713      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
714
715      !!----------------------------------------------------------------------
716      ALLOCATE( znnbrs(jpi,jpj), stat = northcomms_alloc )
717      ALLOCATE( lrankset(jpnij), stat = northcomms_alloc )
718      IF( northcomms_alloc /= 0 ) THEN
719         WRITE(numout,cform_war)
720         WRITE(numout,*) 'northcomms_alloc : failed to allocate arrays'
721         CALL ctl_stop( 'STOP', 'nemo_northcomms : unable to allocate temporary arrays' )
722      ENDIF
723      nsndto = 0
724      isendto = -1
725      ijpj   = 4
726      !
727      ! This routine has been called because ln_nnogather has been set true ( nammpp )
728      ! However, these first few exchanges have to use the mpi_allgather method to
729      ! establish the neighbour lists to use in subsequent peer to peer exchanges.
730      ! Consequently, set l_north_nogather to be false here and set it true only after
731      ! the lists have been established.
732      !
733      l_north_nogather = .FALSE.
734      !
735      ! Exchange and store ranks on northern rows
736
737      DO jtyp = 1,4
738
739         lrankset = .FALSE.
740         znnbrs = narea
741         SELECT CASE (jtyp)
742            CASE(1)
743               CALL lbc_lnk( znnbrs, 'T', 1. )      ! Type 1: T,W-points
744            CASE(2)
745               CALL lbc_lnk( znnbrs, 'U', 1. )      ! Type 2: U-point
746            CASE(3)
747               CALL lbc_lnk( znnbrs, 'V', 1. )      ! Type 3: V-point
748            CASE(4)
749               CALL lbc_lnk( znnbrs, 'F', 1. )      ! Type 4: F-point
750         END SELECT
751
752         IF ( njmppt(narea) .EQ. MAXVAL( njmppt ) ) THEN
753            DO jj = nlcj-ijpj+1, nlcj
754               ij = jj - nlcj + ijpj
755               DO ji = 1,jpi
756                  IF ( INT(znnbrs(ji,jj)) .NE. 0 .AND. INT(znnbrs(ji,jj)) .NE. narea ) &
757               &     lrankset(INT(znnbrs(ji,jj))) = .true.
758               END DO
759            END DO
760
761            DO jj = 1,jpnij
762               IF ( lrankset(jj) ) THEN
763                  nsndto(jtyp) = nsndto(jtyp) + 1
764                  IF ( nsndto(jtyp) .GT. jpmaxngh ) THEN
765                     CALL ctl_stop( ' Too many neighbours in nemo_northcomms ', &
766                  &                 ' jpmaxngh will need to be increased ')
767                  ENDIF
768                  isendto(nsndto(jtyp),jtyp) = jj-1   ! narea converted to MPI rank
769               ENDIF
770            END DO
771         ENDIF
772
773      END DO
774
775      !
776      ! Type 5: I-point
777      !
778      ! ICE point exchanges may involve some averaging. The neighbours list is
779      ! built up using two exchanges to ensure that the whole stencil is covered.
780      ! lrankset should not be reset between these 'J' and 'K' point exchanges
781
782      jtyp = 5
783      lrankset = .FALSE.
784      znnbrs = narea
785      CALL lbc_lnk( znnbrs, 'J', 1. ) ! first ice U-V point
786
787      IF ( njmppt(narea) .EQ. MAXVAL( njmppt ) ) THEN
788         DO jj = nlcj-ijpj+1, nlcj
789            ij = jj - nlcj + ijpj
790            DO ji = 1,jpi
791               IF ( INT(znnbrs(ji,jj)) .NE. 0 .AND. INT(znnbrs(ji,jj)) .NE. narea ) &
792            &     lrankset(INT(znnbrs(ji,jj))) = .true.
793         END DO
794        END DO
795      ENDIF
796
797      znnbrs = narea
798      CALL lbc_lnk( znnbrs, 'K', 1. ) ! second ice U-V point
799
800      IF ( njmppt(narea) .EQ. MAXVAL( njmppt )) THEN
801         DO jj = nlcj-ijpj+1, nlcj
802            ij = jj - nlcj + ijpj
803            DO ji = 1,jpi
804               IF ( INT(znnbrs(ji,jj)) .NE. 0 .AND.  INT(znnbrs(ji,jj)) .NE. narea ) &
805            &       lrankset( INT(znnbrs(ji,jj))) = .true.
806            END DO
807         END DO
808
809         DO jj = 1,jpnij
810            IF ( lrankset(jj) ) THEN
811               nsndto(jtyp) = nsndto(jtyp) + 1
812               IF ( nsndto(jtyp) .GT. jpmaxngh ) THEN
813                  CALL ctl_stop( ' Too many neighbours in nemo_northcomms ', &
814               &                 ' jpmaxngh will need to be increased ')
815               ENDIF
816               isendto(nsndto(jtyp),jtyp) = jj-1   ! narea converted to MPI rank
817            ENDIF
818         END DO
819         !
820         ! For northern row areas, set l_north_nogather so that all subsequent exchanges
821         ! can use peer to peer communications at the north fold
822         !
823         l_north_nogather = .TRUE.
824         !
825      ENDIF
826      DEALLOCATE( znnbrs )
827      DEALLOCATE( lrankset )
828
829   END SUBROUTINE nemo_northcomms
830#else
831   SUBROUTINE nemo_northcomms      ! Dummy routine
832      WRITE(*,*) 'nemo_northcomms: You should not have seen this print! error?'
833   END SUBROUTINE nemo_northcomms
834#endif
835#else
836CONTAINS
837   SUBROUTINE nemo_gcm_tam
838      WRITE(*,*) 'nemo_gcm_tam: You should not have seen this print! error?'
839   END SUBROUTINE nemo_gcm_tam
840#endif
841   !!======================================================================
842END MODULE nemogcm_tam
Note: See TracBrowser for help on using the repository browser.