source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90 @ 11277

Last change on this file since 11277 was 11277, checked in by kingr, 17 months ago

Merged Juan's changes for running AMM15 woth wave coupling.
Corrected minor logic error to allow AMM7-uncoupled to reproduce earlier results.
Few line spacing changes to allow merging with OBS br and trunk rvn 5518.

File size: 41.4 KB
Line 
1MODULE nemogcm
2   !!======================================================================
3   !!                       ***  MODULE nemogcm   ***
4   !! Ocean system   : NEMO GCM (ocean dynamics, on-line tracers, biochemistry and sea-ice)
5   !!======================================================================
6   !! History :  OPA  ! 1990-10  (C. Levy, G. Madec)  Original code
7   !!            7.0  ! 1991-11  (M. Imbard, C. Levy, G. Madec)
8   !!            7.1  ! 1993-03  (M. Imbard, C. Levy, G. Madec, O. Marti, M. Guyon, A. Lazar,
9   !!                             P. Delecluse, C. Perigaud, G. Caniaux, B. Colot, C. Maes) release 7.1
10   !!             -   ! 1992-06  (L.Terray)  coupling implementation
11   !!             -   ! 1993-11  (M.A. Filiberti) IGLOO sea-ice
12   !!            8.0  ! 1996-03  (M. Imbard, C. Levy, G. Madec, O. Marti, M. Guyon, A. Lazar,
13   !!                             P. Delecluse, L.Terray, M.A. Filiberti, J. Vialar, A.M. Treguier, M. Levy) release 8.0
14   !!            8.1  ! 1997-06  (M. Imbard, G. Madec)
15   !!            8.2  ! 1999-11  (M. Imbard, H. Goosse)  LIM sea-ice model
16   !!                 ! 1999-12  (V. Thierry, A-M. Treguier, M. Imbard, M-A. Foujols)  OPEN-MP
17   !!                 ! 2000-07  (J-M Molines, M. Imbard)  Open Boundary Conditions  (CLIPPER)
18   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and modules
19   !!             -   ! 2004-06  (R. Redler, NEC CCRLE, Germany) add OASIS[3/4] coupled interfaces
20   !!             -   ! 2004-08  (C. Talandier) New trends organization
21   !!             -   ! 2005-06  (C. Ethe) Add the 1D configuration possibility
22   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
23   !!             -   ! 2006-03  (L. Debreu, C. Mazauric)  Agrif implementation
24   !!             -   ! 2006-04  (G. Madec, R. Benshila)  Step reorganization
25   !!             -   ! 2007-07  (J. Chanut, A. Sellar) Unstructured open boundaries (BDY)
26   !!            3.2  ! 2009-08  (S. Masson)  open/write in the listing file in mpp
27   !!            3.3  ! 2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
28   !!             -   ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
29   !!            3.3.1! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation
30   !!            3.4  ! 2011-11  (C. Harris) decomposition changes for running with CICE
31   !!                 ! 2012-05  (C. Calone, J. Simeon, G. Madec, C. Ethe) Add grid coarsening
32   !!----------------------------------------------------------------------
33
34   !!----------------------------------------------------------------------
35   !!   nemo_gcm       : solve ocean dynamics, tracer, biogeochemistry and/or sea-ice
36   !!   nemo_init      : initialization of the NEMO system
37   !!   nemo_ctl       : initialisation of the contol print
38   !!   nemo_closefile : close remaining open files
39   !!   nemo_alloc     : dynamical allocation
40   !!   nemo_partition : calculate MPP domain decomposition
41   !!   factorise      : calculate the factors of the no. of MPI processes
42   !!----------------------------------------------------------------------
43   USE step_oce        ! module used in the ocean time stepping module
44   USE cla             ! cross land advection               (tra_cla routine)
45   USE domcfg          ! domain configuration               (dom_cfg routine)
46   USE mppini          ! shared/distributed memory setting (mpp_init routine)
47   USE domain          ! domain initialization             (dom_init routine)
48#if defined key_nemocice_decomp
49   USE ice_domain_size, only: nx_global, ny_global
50#endif
51   USE tideini         ! tidal components initialization   (tide_ini routine)
52   USE bdyini          ! open boundary cond. setting       (bdy_init routine)
53   USE bdydta          ! open boundary cond. setting   (bdy_dta_init routine)
54   USE bdytides        ! open boundary cond. setting   (bdytide_init routine)
55   USE istate          ! initial state setting          (istate_init routine)
56   USE ldfdyn          ! lateral viscosity setting      (ldfdyn_init routine)
57   USE ldftra          ! lateral diffusivity setting    (ldftra_init routine)
58   USE zdfini          ! vertical physics setting          (zdf_init routine)
59   USE phycst          ! physical constant                  (par_cst routine)
60   USE trdini          ! dyn/tra trends initialization     (trd_init routine)
61   USE asminc          ! assimilation increments     
62   USE asmbkg          ! writing out state trajectory
63   USE asmbgc          ! biogeochemical assimilation increments
64   USE diaptr          ! poleward transports           (dia_ptr_init routine)
65   USE diadct          ! sections transports           (dia_dct_init routine)
66   USE diaobs          ! Observation diagnostics       (dia_obs_init routine)
67   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
68   USE step            ! NEMO time-stepping                 (stp     routine)
69   USE icbini          ! handle bergs, initialisation
70   USE icbstp          ! handle bergs, calving, themodynamics and transport
71   USE sbccpl
72   USE cpl_oasis3      ! OASIS3 coupling
73   USE c1d             ! 1D configuration
74   USE step_c1d        ! Time stepping loop for the 1D configuration
75   USE dyndmp          ! Momentum damping
76#if defined key_top
77   USE trcini          ! passive tracer initialisation
78#endif
79   USE lib_mpp         ! distributed memory computing
80#if defined key_iomput
81   USE xios
82#endif
83   USE sbctide, ONLY: lk_tide
84   USE crsini          ! initialise grid coarsening utility
85   USE lbcnfd, ONLY: isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges
86   USE sbc_oce, ONLY: lk_oasis
87   USE stopar
88   USE stopts
89   USE diatmb          ! Top,middle,bottom output
90   USE dia25h          ! 25h mean output
91   USE diaopfoam       ! FOAM operational output
92   USE diurnal_bulk    ! diurnal bulk SST
93
94   IMPLICIT NONE
95   PRIVATE
96
97   PUBLIC   nemo_gcm    ! called by model.F90
98   PUBLIC   nemo_init   ! needed by AGRIF
99   PUBLIC   nemo_alloc  ! needed by TAM
100
101   CHARACTER(lc) ::   cform_aaa="( /, 'AAAAAAAA', / ) "     ! flag for output listing
102
103   !!----------------------------------------------------------------------
104   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
105   !! $Id$
106   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
107   !!----------------------------------------------------------------------
108CONTAINS
109
110   SUBROUTINE nemo_gcm
111      !!----------------------------------------------------------------------
112      !!                     ***  ROUTINE nemo_gcm  ***
113      !!
114      !! ** Purpose :   NEMO solves the primitive equations on an orthogonal
115      !!              curvilinear mesh on the sphere.
116      !!
117      !! ** Method  : - model general initialization
118      !!              - launch the time-stepping (stp routine)
119      !!              - finalize the run by closing files and communications
120      !!
121      !! References : Madec, Delecluse, Imbard, and Levy, 1997:  internal report, IPSL.
122      !!              Madec, 2008, internal report, IPSL.
123      !!----------------------------------------------------------------------
124      INTEGER ::   istp       ! time step index
125      !!----------------------------------------------------------------------
126      !
127#if defined key_agrif
128      CALL Agrif_Init_Grids()      ! AGRIF: set the meshes
129#endif
130
131      !                            !-----------------------!
132      CALL nemo_init               !==  Initialisations  ==!
133      !                            !-----------------------!
134#if defined key_agrif
135      CALL Agrif_Declare_Var_dom   ! AGRIF: set the meshes for DOM
136      CALL Agrif_Declare_Var       !  "      "   "   "      "  DYN/TRA
137# if defined key_top
138      CALL Agrif_Declare_Var_top   !  "      "   "   "      "  TOP
139# endif
140# if defined key_lim2
141      CALL Agrif_Declare_Var_lim2  !  "      "   "   "      "  LIM
142# endif
143#endif
144      ! check that all process are still there... If some process have an error,
145      ! they will never enter in step and other processes will wait until the end of the cpu time!
146      IF( lk_mpp )   CALL mpp_max( nstop )
147
148      IF(lwp) WRITE(numout,cform_aaa)   ! Flag AAAAAAA
149
150      !                            !-----------------------!
151      !                            !==   time stepping   ==!
152      !                            !-----------------------!
153      istp = nit000
154#if defined key_c1d
155         DO WHILE ( istp <= nitend .AND. nstop == 0 )
156            CALL stp_c1d( istp )
157            istp = istp + 1
158         END DO
159#else
160          IF( lk_asminc ) THEN
161             IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 )    ! Output background fields
162             IF( ln_asmdin ) THEN                        ! Direct initialization
163                IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 )    ! Tracers
164                IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics
165                IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH
166                IF( lk_bgcinc ) CALL bgc_asm_inc( nit000 - 1 )    ! BGC
167             ENDIF
168          ENDIF
169
170#if defined key_agrif
171          CALL Agrif_Regrid()
172#endif
173
174         DO WHILE ( istp <= nitend .AND. nstop == 0 )
175#if defined key_agrif
176            CALL stp                         ! AGRIF: time stepping
177#else
178            IF (lk_oasis) CALL sbc_cpl_snd( istp )  ! Coupling to atmos   
179            CALL stp( istp )   
180            ! We don't couple on the final timestep because   
181            ! our restart file has already been written   
182            ! and contains all the necessary data for a   
183            ! restart. sbc_cpl_snd could be called here   
184            ! but it would require   
185            ! a) A test to ensure it was not performed   
186            !    on the very last time-step   
187            ! b) the presence of another call to   
188            !    sbc_cpl_snd call prior to the main DO loop   
189            ! This solution produces identical results   
190            ! with fewer lines of code. 
191#endif
192            istp = istp + 1
193            IF( lk_mpp )   CALL mpp_max( nstop )
194         END DO
195#endif
196
197      IF( lk_diaobs   )   CALL dia_obs_wri
198      !
199      IF( ( lk_asminc ).AND.( ln_balwri ) ) CALL asm_bgc_bal_wri( nitend )  ! Output balancing increments
200      !
201      IF( ln_icebergs )   CALL icb_end( nitend )
202
203      !                            !------------------------!
204      !                            !==  finalize the run  ==!
205      !                            !------------------------!
206      IF(lwp) WRITE(numout,cform_aaa)   ! Flag AAAAAAA
207      !
208      IF( nstop /= 0 .AND. lwp ) THEN   ! error print
209         WRITE(numout,cform_err)
210         WRITE(numout,*) nstop, ' error have been found'
211      ENDIF
212      !
213#if defined key_agrif
214      IF( .NOT. Agrif_Root() ) THEN
215         CALL Agrif_ParentGrid_To_ChildGrid()
216         IF( lk_diaobs ) CALL dia_obs_wri
217         IF( nn_timing == 1 )   CALL timing_finalize
218         CALL Agrif_ChildGrid_To_ParentGrid()
219      ENDIF
220#endif
221      IF( nn_timing == 1 )   CALL timing_finalize
222      !
223      CALL nemo_closefile
224      !
225#if defined key_iomput
226      CALL xios_finalize                ! end mpp communications with xios
227      IF( lk_oasis ) CALL cpl_finalize    ! end coupling and mpp communications with OASIS
228#else
229      IF( lk_oasis ) THEN
230         CALL cpl_finalize              ! end coupling and mpp communications with OASIS
231      ELSE
232         IF( lk_mpp )   CALL mppstop    ! end mpp communications
233      ENDIF
234#endif
235      !
236   END SUBROUTINE nemo_gcm
237
238
239   SUBROUTINE nemo_init
240      !!----------------------------------------------------------------------
241      !!                     ***  ROUTINE nemo_init  ***
242      !!
243      !! ** Purpose :   initialization of the NEMO GCM
244      !!----------------------------------------------------------------------
245      INTEGER ::   ji            ! dummy loop indices
246      INTEGER ::   ilocal_comm   ! local integer
247      INTEGER ::   ios
248      CHARACTER(len=80), DIMENSION(16) ::   cltxt
249      !
250      NAMELIST/namctl/ ln_ctl  , nn_print, nn_ictls, nn_ictle,   &
251         &             nn_isplt, nn_jsplt, nn_jctls, nn_jctle,   &
252         &             nn_bench, nn_timing
253      NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
254         &             jpizoom, jpjzoom, jperio, ln_use_jattr
255      !!----------------------------------------------------------------------
256      !
257      cltxt = ''
258      cxios_context = 'nemo'
259      !
260      !                             ! Open reference namelist and configuration namelist files
261      CALL ctl_opn( numnam_ref, 'namelist_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. )
262      CALL ctl_opn( numnam_cfg, 'namelist_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. )
263      !
264      REWIND( numnam_ref )              ! Namelist namctl in reference namelist : Control prints & Benchmark
265      READ  ( numnam_ref, namctl, IOSTAT = ios, ERR = 901 )
266901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namctl in reference namelist', .TRUE. )
267
268      REWIND( numnam_cfg )              ! Namelist namctl in confguration namelist : Control prints & Benchmark
269      READ  ( numnam_cfg, namctl, IOSTAT = ios, ERR = 902 )
270902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namctl in configuration namelist', .TRUE. )
271
272      !
273      REWIND( numnam_ref )              ! Namelist namcfg in reference namelist : Control prints & Benchmark
274      READ  ( numnam_ref, namcfg, IOSTAT = ios, ERR = 903 )
275903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in reference namelist', .TRUE. )
276
277      REWIND( numnam_cfg )              ! Namelist namcfg in confguration namelist : Control prints & Benchmark
278      READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 )
279904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. )   
280
281! Force values for AGRIF zoom (cf. agrif_user.F90)
282#if defined key_agrif
283   IF( .NOT. Agrif_Root() ) THEN
284      jpiglo  = nbcellsx + 2 + 2*nbghostcells
285      jpjglo  = nbcellsy + 2 + 2*nbghostcells
286      jpi     = ( jpiglo-2*jpreci + (jpni-1+0) ) / jpni + 2*jpreci
287      jpj     = ( jpjglo-2*jprecj + (jpnj-1+0) ) / jpnj + 2*jprecj
288      jpidta  = jpiglo
289      jpjdta  = jpjglo
290      jpizoom = 1
291      jpjzoom = 1
292      nperio  = 0
293      jperio  = 0
294      ln_use_jattr = .false.
295   ENDIF
296#endif
297      !
298      !                             !--------------------------------------------!
299      !                             !  set communicator & select the local node  !
300      !                             !  NB: mynode also opens output.namelist.dyn !
301      !                             !      on unit number numond on first proc   !
302      !                             !--------------------------------------------!
303#if defined key_iomput
304      IF( Agrif_Root() ) THEN
305         IF( lk_oasis ) THEN
306            CALL cpl_init( "oceanx", ilocal_comm )                     ! nemo local communicator given by oasis
307            CALL xios_initialize( "not used",local_comm=ilocal_comm )    ! send nemo communicator to xios
308         ELSE
309            CALL  xios_initialize( "for_xios_mpi_id",return_comm=ilocal_comm )    ! nemo local communicator given by xios
310         ENDIF
311      ENDIF
312      ! Nodes selection (control print return in cltxt)
313      narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )
314#else
315      IF( lk_oasis ) THEN
316         IF( Agrif_Root() ) THEN
317            CALL cpl_init( "oceanx", ilocal_comm )                      ! nemo local communicator given by oasis
318         ENDIF
319         ! Nodes selection (control print return in cltxt)
320         narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )
321      ELSE
322         ilocal_comm = 0
323         ! Nodes selection (control print return in cltxt)
324         narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop )
325      ENDIF
326#endif
327      narea = narea + 1                                     ! mynode return the rank of proc (0 --> jpnij -1 )
328
329      lwm = (narea == 1)                                    ! control of output namelists
330      lwp = (narea == 1) .OR. ln_ctl                        ! control of all listing output print
331
332      IF(lwm) THEN
333         ! write merged namelists from earlier to output namelist now that the
334         ! file has been opened in call to mynode. nammpp has already been
335         ! written in mynode (if lk_mpp_mpi)
336         WRITE( numond, namctl )
337         WRITE( numond, namcfg )
338      ENDIF
339
340      ! If dimensions of processor grid weren't specified in the namelist file
341      ! then we calculate them here now that we have our communicator size
342      IF( (jpni < 1) .OR. (jpnj < 1) )THEN
343#if   defined key_mpp_mpi
344         IF( Agrif_Root() ) CALL nemo_partition(mppsize)
345#else
346         jpni  = 1
347         jpnj  = 1
348         jpnij = jpni*jpnj
349#endif
350      END IF
351
352      ! Calculate domain dimensions given calculated jpni and jpnj
353      ! This used to be done in par_oce.F90 when they were parameters rather
354      ! than variables
355      IF( Agrif_Root() ) THEN
356#if defined key_nemocice_decomp
357         jpi = ( nx_global+2-2*jpreci + (jpni-1) ) / jpni + 2*jpreci ! first  dim.
358         jpj = ( ny_global+2-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj ! second dim.
359#else
360         jpi = ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci   ! first  dim.
361         jpj = ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj   ! second dim.
362#endif
363      ENDIF         
364         jpk = jpkdta                                             ! third dim
365#if defined key_agrif
366         ! simple trick to use same vertical grid as parent
367         ! but different number of levels:
368         ! Save maximum number of levels in jpkdta, then define all vertical grids
369         ! with this number.
370         ! Suppress once vertical online interpolation is ok
371         IF(.NOT.Agrif_Root()) jpkdta = Agrif_Parent(jpkdta)
372#endif
373         jpim1 = jpi-1                                            ! inner domain indices
374         jpjm1 = jpj-1                                            !   "           "
375         jpkm1 = jpk-1                                            !   "           "
376         jpij  = jpi*jpj                                          !  jpi x j
377
378      IF(lwp) THEN                            ! open listing units
379         !
380         CALL ctl_opn( numout, 'ocean.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
381         !
382         WRITE(numout,*)
383         WRITE(numout,*) '   CNRS - NERC - Met OFFICE - MERCATOR-ocean - INGV - CMCC'
384         WRITE(numout,*) '                       NEMO team'
385         WRITE(numout,*) '            Ocean General Circulation Model'
386         WRITE(numout,*) '                  version 3.6  (2015) '
387         WRITE(numout,*)
388         WRITE(numout,*)
389         DO ji = 1, SIZE(cltxt)
390            IF( TRIM(cltxt(ji)) /= '' )   WRITE(numout,*) cltxt(ji)      ! control print of mynode
391         END DO
392         WRITE(numout,cform_aaa)                                         ! Flag AAAAAAA
393         !
394      ENDIF
395
396      ! Now we know the dimensions of the grid and numout has been set we can
397      ! allocate arrays
398      CALL nemo_alloc()
399
400      !                             !-------------------------------!
401      !                             !  NEMO general initialization  !
402      !                             !-------------------------------!
403
404      CALL nemo_ctl                          ! Control prints & Benchmark
405
406      !                                      ! Domain decomposition
407      IF( jpni*jpnj == jpnij ) THEN   ;   CALL mpp_init      ! standard cutting out
408      ELSE                            ;   CALL mpp_init2     ! eliminate land processors
409      ENDIF
410      !
411      IF( nn_timing == 1 )  CALL timing_init
412      !
413      !                                      ! General initialization
414                            CALL     phy_cst    ! Physical constants
415                            CALL     eos_init   ! Equation of state
416      IF( lk_c1d        )   CALL     c1d_init   ! 1D column configuration
417                            CALL     dom_cfg    ! Domain configuration
418                            CALL     dom_init   ! Domain
419
420      IF( ln_nnogather )    CALL nemo_northcomms   ! Initialise the northfold neighbour lists (must be done after the masks are defined)
421
422      IF( ln_ctl        )   CALL prt_ctl_init   ! Print control
423
424                            CALL  istate_init   ! ocean initial state (Dynamics and tracers)
425
426      CALL diurnal_sst_bulk_init                ! diurnal sst
427      IF ( ln_diurnal ) CALL diurnal_sst_coolskin_init   ! cool skin 
428
429      IF( lk_tide       )   CALL    tide_init( nit000 )    ! Initialisation of the tidal harmonics
430
431                            CALL     sbc_init   ! Forcings : surface module (clem: moved here for bdy purpose)
432
433      IF( lk_bdy        )   CALL     bdy_init   ! Open boundaries initialisation
434      IF( lk_bdy        )   CALL bdy_dta_init   ! Open boundaries initialisation of external data arrays
435      IF( lk_bdy .AND. lk_tide )   &
436         &                  CALL bdytide_init   ! Open boundaries initialisation of tidal harmonic forcing
437
438                            CALL dyn_nept_init  ! simplified form of Neptune effect
439      !     
440      IF( ln_crs        )   CALL     crs_init   ! Domain initialization of coarsened grid
441      !
442                                ! Ocean physics
443      !                                         ! Vertical physics
444                            CALL     zdf_init      ! namelist read
445                            CALL zdf_bfr_init      ! bottom friction
446      IF( lk_zdfric     )   CALL zdf_ric_init      ! Richardson number dependent Kz
447      IF( lk_zdftke     )   CALL zdf_tke_init      ! TKE closure scheme
448      IF( lk_zdfgls     )   CALL zdf_gls_init      ! GLS closure scheme
449      IF( lk_zdfkpp     )   CALL zdf_kpp_init      ! KPP closure scheme
450      IF( lk_zdftmx     )   CALL zdf_tmx_init      ! tidal vertical mixing
451      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
452         &                  CALL zdf_ddm_init      ! double diffusive mixing
453      !                                         ! Lateral physics
454                            CALL ldf_tra_init      ! Lateral ocean tracer physics
455                            CALL ldf_dyn_init      ! Lateral ocean momentum physics
456      IF( lk_ldfslp     )   CALL ldf_slp_init      ! slope of lateral mixing
457
458      !                                     ! Active tracers
459                            CALL tra_qsr_init   ! penetrative solar radiation qsr
460                            CALL tra_bbc_init   ! bottom heat flux
461      IF( lk_trabbl     )   CALL tra_bbl_init   ! advective (and/or diffusive) bottom boundary layer scheme
462                            CALL tra_dmp_init   ! internal damping trends- tracers
463                            CALL tra_adv_init   ! horizontal & vertical advection
464                            CALL tra_ldf_init   ! lateral mixing
465                            CALL tra_zdf_init   ! vertical mixing and after tracer fields
466
467      !                                     ! Dynamics
468      IF( lk_c1d        )   CALL dyn_dmp_init   ! internal damping trends- momentum
469                            CALL dyn_adv_init   ! advection (vector or flux form)
470                            CALL dyn_vor_init   ! vorticity term including Coriolis
471                            CALL dyn_ldf_init   ! lateral mixing
472                            CALL dyn_hpg_init   ! horizontal gradient of Hydrostatic pressure
473                            CALL dyn_zdf_init   ! vertical diffusion
474                            CALL dyn_spg_init   ! surface pressure gradient
475
476      !                                     ! Misc. options
477      IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_init       ! Cross Land Advection
478                            CALL icb_init( rdt, nit000)   ! initialise icebergs instance
479                            CALL sto_par_init   ! Stochastic parametrization
480      IF( ln_sto_eos     )  CALL sto_pts_init   ! RRandom T/S fluctuations
481     
482#if defined key_top
483      !                                     ! Passive tracers
484                            CALL     trc_init
485#endif
486      !                                     ! Diagnostics
487      IF( lk_floats     )   CALL     flo_init   ! drifting Floats
488      IF( lk_diaar5     )   CALL dia_ar5_init   ! ar5 diag
489                            CALL dia_ptr_init   ! Poleward TRansports initialization
490      IF( lk_diadct     )   CALL dia_dct_init   ! Sections tranports
491                            CALL dia_hsb_init   ! heat content, salt content and volume budgets
492                            CALL     trd_init   ! Mixed-layer/Vorticity/Integral constraints trends
493      IF( lk_diaobs     ) THEN                  ! Observation & model comparison
494                            CALL dia_obs_init            ! Initialize observational data
495                            CALL dia_obs( nit000 - 1 )   ! Observation operator for restart
496      ENDIF
497
498      !                                     ! Assimilation increments
499      IF( lk_asminc ) THEN 
500#if defined key_shelf 
501         CALL  zdf_mxl_tref()     ! Initialization of hmld_tref
502#endif
503         CALL asm_inc_init     ! Initialize assimilation increments
504      ENDIF
505
506      IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler
507                            CALL dia_tmb_init  ! TMB outputs
508                            CALL dia_25h_init  ! 25h mean  outputs
509                            CALL dia_diaopfoam_init  ! FOAM operational output
510      !
511   END SUBROUTINE nemo_init
512
513
514   SUBROUTINE nemo_ctl
515      !!----------------------------------------------------------------------
516      !!                     ***  ROUTINE nemo_ctl  ***
517      !!
518      !! ** Purpose :   control print setting
519      !!
520      !! ** Method  : - print namctl information and check some consistencies
521      !!----------------------------------------------------------------------
522      !
523      IF(lwp) THEN                  ! control print
524         WRITE(numout,*)
525         WRITE(numout,*) 'nemo_ctl: Control prints & Benchmark'
526         WRITE(numout,*) '~~~~~~~ '
527         WRITE(numout,*) '   Namelist namctl'
528         WRITE(numout,*) '      run control (for debugging)     ln_ctl     = ', ln_ctl
529         WRITE(numout,*) '      level of print                  nn_print   = ', nn_print
530         WRITE(numout,*) '      Start i indice for SUM control  nn_ictls   = ', nn_ictls
531         WRITE(numout,*) '      End i indice for SUM control    nn_ictle   = ', nn_ictle
532         WRITE(numout,*) '      Start j indice for SUM control  nn_jctls   = ', nn_jctls
533         WRITE(numout,*) '      End j indice for SUM control    nn_jctle   = ', nn_jctle
534         WRITE(numout,*) '      number of proc. following i     nn_isplt   = ', nn_isplt
535         WRITE(numout,*) '      number of proc. following j     nn_jsplt   = ', nn_jsplt
536         WRITE(numout,*) '      benchmark parameter (0/1)       nn_bench   = ', nn_bench
537         WRITE(numout,*) '      timing activated    (0/1)       nn_timing  = ', nn_timing
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
549      IF(lwp) THEN                  ! control print
550         WRITE(numout,*)
551         WRITE(numout,*) 'namcfg  : configuration initialization through namelist read'
552         WRITE(numout,*) '~~~~~~~ '
553         WRITE(numout,*) '   Namelist namcfg'
554         WRITE(numout,*) '      configuration name              cp_cfg      = ', TRIM(cp_cfg)
555         WRITE(numout,*) '      configuration zoom name         cp_cfz      = ', TRIM(cp_cfz)
556         WRITE(numout,*) '      configuration resolution        jp_cfg      = ', jp_cfg
557         WRITE(numout,*) '      1st lateral dimension ( >= jpi ) jpidta     = ', jpidta
558         WRITE(numout,*) '      2nd    "         "    ( >= jpj ) jpjdta     = ', jpjdta
559         WRITE(numout,*) '      3nd    "         "               jpkdta     = ', jpkdta
560         WRITE(numout,*) '      1st dimension of global domain in i jpiglo  = ', jpiglo
561         WRITE(numout,*) '      2nd    -                  -    in j jpjglo  = ', jpjglo
562         WRITE(numout,*) '      left bottom i index of the zoom (in data domain) jpizoom = ', jpizoom
563         WRITE(numout,*) '      left bottom j index of the zoom (in data domain) jpizoom = ', jpjzoom
564         WRITE(numout,*) '      lateral cond. type (between 0 and 6) jperio = ', jperio   
565         WRITE(numout,*) '      use file attribute if exists as i/p j-start ln_use_jattr = ', ln_use_jattr
566      ENDIF
567      !                             ! Parameter control
568      !
569      IF( ln_ctl ) THEN                 ! sub-domain area indices for the control prints
570         IF( lk_mpp .AND. jpnij > 1 ) THEN
571            isplt = jpni   ;   jsplt = jpnj   ;   ijsplt = jpni*jpnj   ! the domain is forced to the real split domain
572         ELSE
573            IF( isplt == 1 .AND. jsplt == 1  ) THEN
574               CALL ctl_warn( ' - isplt & jsplt are equal to 1',   &
575                  &           ' - the print control will be done over the whole domain' )
576            ENDIF
577            ijsplt = isplt * jsplt            ! total number of processors ijsplt
578         ENDIF
579         IF(lwp) WRITE(numout,*)'          - The total number of processors over which the'
580         IF(lwp) WRITE(numout,*)'            print control will be done is ijsplt : ', ijsplt
581         !
582         !                              ! indices used for the SUM control
583         IF( nictls+nictle+njctls+njctle == 0 )   THEN    ! print control done over the default area
584            lsp_area = .FALSE.
585         ELSE                                             ! print control done over a specific  area
586            lsp_area = .TRUE.
587            IF( nictls < 1 .OR. nictls > jpiglo )   THEN
588               CALL ctl_warn( '          - nictls must be 1<=nictls>=jpiglo, it is forced to 1' )
589               nictls = 1
590            ENDIF
591            IF( nictle < 1 .OR. nictle > jpiglo )   THEN
592               CALL ctl_warn( '          - nictle must be 1<=nictle>=jpiglo, it is forced to jpiglo' )
593               nictle = jpiglo
594            ENDIF
595            IF( njctls < 1 .OR. njctls > jpjglo )   THEN
596               CALL ctl_warn( '          - njctls must be 1<=njctls>=jpjglo, it is forced to 1' )
597               njctls = 1
598            ENDIF
599            IF( njctle < 1 .OR. njctle > jpjglo )   THEN
600               CALL ctl_warn( '          - njctle must be 1<=njctle>=jpjglo, it is forced to jpjglo' )
601               njctle = jpjglo
602            ENDIF
603         ENDIF
604      ENDIF
605      !
606      IF( nbench == 1 ) THEN              ! Benchmark
607         SELECT CASE ( cp_cfg )
608         CASE ( 'gyre' )   ;   CALL ctl_warn( ' The Benchmark is activated ' )
609         CASE DEFAULT      ;   CALL ctl_stop( ' The Benchmark is based on the GYRE configuration:',   &
610            &                                 ' cp_cfg = "gyre" in namelist &namcfg or set nbench = 0' )
611         END SELECT
612      ENDIF
613      !
614      IF( 1_wp /= SIGN(1._wp,-0._wp)  )   CALL ctl_stop( 'nemo_ctl: The intrinsec SIGN function follows ',  &
615         &                                               'f2003 standard. '                              ,  &
616         &                                               'Compile with key_nosignedzero enabled' )
617      !
618   END SUBROUTINE nemo_ctl
619
620
621   SUBROUTINE nemo_closefile
622      !!----------------------------------------------------------------------
623      !!                     ***  ROUTINE nemo_closefile  ***
624      !!
625      !! ** Purpose :   Close the files
626      !!----------------------------------------------------------------------
627      !
628      IF( lk_mpp )   CALL mppsync
629      !
630      CALL iom_close                                 ! close all input/output files managed by iom_*
631      !
632      IF( numstp          /= -1 )   CLOSE( numstp          )   ! time-step file
633      IF( numsol          /= -1 )   CLOSE( numsol          )   ! solver file
634      IF( numnam_ref      /= -1 )   CLOSE( numnam_ref      )   ! oce reference namelist
635      IF( numnam_cfg      /= -1 )   CLOSE( numnam_cfg      )   ! oce configuration namelist
636      IF( lwm.AND.numond  /= -1 )   CLOSE( numond          )   ! oce output namelist
637      IF( numnam_ice_ref  /= -1 )   CLOSE( numnam_ice_ref  )   ! ice reference namelist
638      IF( numnam_ice_cfg  /= -1 )   CLOSE( numnam_ice_cfg  )   ! ice configuration namelist
639      IF( lwm.AND.numoni  /= -1 )   CLOSE( numoni          )   ! ice output namelist
640      IF( numevo_ice      /= -1 )   CLOSE( numevo_ice      )   ! ice variables (temp. evolution)
641      IF( numout          /=  6 )   CLOSE( numout          )   ! standard model output file
642      IF( numdct_vol      /= -1 )   CLOSE( numdct_vol      )   ! volume transports
643      IF( numdct_heat     /= -1 )   CLOSE( numdct_heat     )   ! heat transports
644      IF( numdct_salt     /= -1 )   CLOSE( numdct_salt     )   ! salt transports
645
646      !
647      numout = 6                                     ! redefine numout in case it is used after this point...
648      !
649   END SUBROUTINE nemo_closefile
650
651
652   SUBROUTINE nemo_alloc
653      !!----------------------------------------------------------------------
654      !!                     ***  ROUTINE nemo_alloc  ***
655      !!
656      !! ** Purpose :   Allocate all the dynamic arrays of the OPA modules
657      !!
658      !! ** Method  :
659      !!----------------------------------------------------------------------
660      USE diawri    , ONLY: dia_wri_alloc
661      USE dom_oce   , ONLY: dom_oce_alloc
662      USE ldfdyn_oce, ONLY: ldfdyn_oce_alloc
663      USE ldftra_oce, ONLY: ldftra_oce_alloc
664      USE trc_oce   , ONLY: trc_oce_alloc
665      USE diainsitutem, ONLY: insitu_tem_alloc
666#if defined key_diadct 
667      USE diadct    , ONLY: diadct_alloc 
668#endif 
669#if defined key_bdy
670      USE bdy_oce   , ONLY: bdy_oce_alloc
671#endif
672      !
673      INTEGER :: ierr
674      !!----------------------------------------------------------------------
675      !
676      ierr =        oce_alloc       ()          ! ocean
677      ierr = ierr + dia_wri_alloc   ()
678      ierr = ierr + dom_oce_alloc   ()          ! ocean domain
679      ierr = ierr + ldfdyn_oce_alloc()          ! ocean lateral  physics : dynamics
680      ierr = ierr + ldftra_oce_alloc()          ! ocean lateral  physics : tracers
681      ierr = ierr + zdf_oce_alloc   ()          ! ocean vertical physics
682      ierr = ierr + insitu_tem_alloc()
683      !
684      ierr = ierr + trc_oce_alloc   ()          ! shared TRC / TRA arrays
685      !
686#if defined key_diadct 
687      ierr = ierr + diadct_alloc    ()          !
688#endif 
689#if defined key_bdy
690      ierr = ierr + bdy_oce_alloc   ()          ! bdy masks (incl. initialization)
691#endif
692      !
693      IF( lk_mpp    )   CALL mpp_sum( ierr )
694      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'nemo_alloc : unable to allocate standard ocean arrays' )
695      !
696   END SUBROUTINE nemo_alloc
697
698
699   SUBROUTINE nemo_partition( num_pes )
700      !!----------------------------------------------------------------------
701      !!                 ***  ROUTINE nemo_partition  ***
702      !!
703      !! ** Purpose :
704      !!
705      !! ** Method  :
706      !!----------------------------------------------------------------------
707      INTEGER, INTENT(in) ::   num_pes   ! The number of MPI processes we have
708      !
709      INTEGER, PARAMETER :: nfactmax = 20
710      INTEGER :: nfact ! The no. of factors returned
711      INTEGER :: ierr  ! Error flag
712      INTEGER :: ji
713      INTEGER :: idiff, mindiff, imin ! For choosing pair of factors that are closest in value
714      INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors
715      !!----------------------------------------------------------------------
716      !
717      ierr = 0
718      !
719      CALL factorise( ifact, nfactmax, nfact, num_pes, ierr )
720      !
721      IF( nfact <= 1 ) THEN
722         WRITE (numout, *) 'WARNING: factorisation of number of PEs failed'
723         WRITE (numout, *) '       : using grid of ',num_pes,' x 1'
724         jpnj = 1
725         jpni = num_pes
726      ELSE
727         ! Search through factors for the pair that are closest in value
728         mindiff = 1000000
729         imin    = 1
730         DO ji = 1, nfact-1, 2
731            idiff = ABS( ifact(ji) - ifact(ji+1) )
732            IF( idiff < mindiff ) THEN
733               mindiff = idiff
734               imin = ji
735            ENDIF
736         END DO
737         jpnj = ifact(imin)
738         jpni = ifact(imin + 1)
739      ENDIF
740      !
741      jpnij = jpni*jpnj
742      !
743   END SUBROUTINE nemo_partition
744
745
746   SUBROUTINE factorise( kfax, kmaxfax, knfax, kn, kerr )
747      !!----------------------------------------------------------------------
748      !!                     ***  ROUTINE factorise  ***
749      !!
750      !! ** Purpose :   return the prime factors of n.
751      !!                knfax factors are returned in array kfax which is of
752      !!                maximum dimension kmaxfax.
753      !! ** Method  :
754      !!----------------------------------------------------------------------
755      INTEGER                    , INTENT(in   ) ::   kn, kmaxfax
756      INTEGER                    , INTENT(  out) ::   kerr, knfax
757      INTEGER, DIMENSION(kmaxfax), INTENT(  out) ::   kfax
758      !
759      INTEGER :: ifac, jl, inu
760      INTEGER, PARAMETER :: ntest = 14
761      INTEGER, DIMENSION(ntest) :: ilfax
762      !
763      ! ilfax contains the set of allowed factors.
764      ilfax(:) = (/(2**jl,jl=ntest,1,-1)/)
765      !!----------------------------------------------------------------------
766      ! ilfax contains the set of allowed factors.
767      ilfax(:) = (/(2**jl,jl=ntest,1,-1)/)
768
769      ! Clear the error flag and initialise output vars
770      kerr = 0
771      kfax = 1
772      knfax = 0
773
774      ! Find the factors of n.
775      IF( kn == 1 )   GOTO 20
776
777      ! nu holds the unfactorised part of the number.
778      ! knfax holds the number of factors found.
779      ! l points to the allowed factor list.
780      ! ifac holds the current factor.
781
782      inu   = kn
783      knfax = 0
784
785      DO jl = ntest, 1, -1
786         !
787         ifac = ilfax(jl)
788         IF( ifac > inu )   CYCLE
789
790         ! Test whether the factor will divide.
791
792         IF( MOD(inu,ifac) == 0 ) THEN
793            !
794            knfax = knfax + 1            ! Add the factor to the list
795            IF( knfax > kmaxfax ) THEN
796               kerr = 6
797               write (*,*) 'FACTOR: insufficient space in factor array ', knfax
798               return
799            ENDIF
800            kfax(knfax) = ifac
801            ! Store the other factor that goes with this one
802            knfax = knfax + 1
803            kfax(knfax) = inu / ifac
804            !WRITE (*,*) 'ARPDBG, factors ',knfax-1,' & ',knfax,' are ', kfax(knfax-1),' and ',kfax(knfax)
805         ENDIF
806         !
807      END DO
808
809   20 CONTINUE      ! Label 20 is the exit point from the factor search loop.
810      !
811   END SUBROUTINE factorise
812
813#if defined key_mpp_mpi
814
815   SUBROUTINE nemo_northcomms
816      !!======================================================================
817      !!                     ***  ROUTINE  nemo_northcomms  ***
818      !! nemo_northcomms    :  Setup for north fold exchanges with explicit
819      !!                       point-to-point messaging
820      !!=====================================================================
821      !!----------------------------------------------------------------------
822      !!
823      !! ** Purpose :   Initialization of the northern neighbours lists.
824      !!----------------------------------------------------------------------
825      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)
826      !!    2.0  ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
827      !!----------------------------------------------------------------------
828
829      INTEGER  ::   sxM, dxM, sxT, dxT, jn
830      INTEGER  ::   njmppmax
831
832      njmppmax = MAXVAL( njmppt )
833   
834      !initializes the north-fold communication variables
835      isendto(:) = 0
836      nsndto = 0
837
838      !if I am a process in the north
839      IF ( njmpp == njmppmax ) THEN
840          !sxM is the first point (in the global domain) needed to compute the
841          !north-fold for the current process
842          sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
843          !dxM is the last point (in the global domain) needed to compute the
844          !north-fold for the current process
845          dxM = jpiglo - nimppt(narea) + 2
846
847          !loop over the other north-fold processes to find the processes
848          !managing the points belonging to the sxT-dxT range
849 
850          DO jn = 1, jpni
851                !sxT is the first point (in the global domain) of the jn
852                !process
853                sxT = nfiimpp(jn, jpnj)
854                !dxT is the last point (in the global domain) of the jn
855                !process
856                dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1
857                IF ((sxM .gt. sxT) .AND. (sxM .lt. dxT)) THEN
858                   nsndto = nsndto + 1
859                     isendto(nsndto) = jn
860                ELSEIF ((sxM .le. sxT) .AND. (dxM .ge. dxT)) THEN
861                   nsndto = nsndto + 1
862                     isendto(nsndto) = jn
863                ELSEIF ((dxM .lt. dxT) .AND. (sxT .lt. dxM)) THEN
864                   nsndto = nsndto + 1
865                     isendto(nsndto) = jn
866                END IF
867          END DO
868          nfsloop = 1
869          nfeloop = nlci
870          DO jn = 2,jpni-1
871           IF(nfipproc(jn,jpnj) .eq. (narea - 1)) THEN
872              IF (nfipproc(jn - 1 ,jpnj) .eq. -1) THEN
873                 nfsloop = nldi
874              ENDIF
875              IF (nfipproc(jn + 1,jpnj) .eq. -1) THEN
876                 nfeloop = nlei
877              ENDIF
878           ENDIF
879        END DO
880
881      ENDIF
882      l_north_nogather = .TRUE.
883   END SUBROUTINE nemo_northcomms
884#else
885   SUBROUTINE nemo_northcomms      ! Dummy routine
886      WRITE(*,*) 'nemo_northcomms: You should not have seen this print! error?'
887   END SUBROUTINE nemo_northcomms
888#endif
889
890   !!======================================================================
891END MODULE nemogcm
892
893
Note: See TracBrowser for help on using the repository browser.