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.
domain.F90 in NEMO/trunk/src/OCE/DOM – NEMO

source: NEMO/trunk/src/OCE/DOM/domain.F90 @ 15267

Last change on this file since 15267 was 15267, checked in by smasson, 3 years ago

trunk: new nogather nolding, #2724
  • Property svn:keywords set to Id
File size: 40.2 KB
Line 
1MODULE domain
2   !!==============================================================================
3   !!                       ***  MODULE domain   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6   !! History :  OPA  !  1990-10  (C. Levy - G. Madec)  Original code
7   !!                 !  1992-01  (M. Imbard) insert time step initialization
8   !!                 !  1996-06  (G. Madec) generalized vertical coordinate
9   !!                 !  1997-02  (G. Madec) creation of domwri.F
10   !!                 !  2001-05  (E.Durand - G. Madec) insert closed sea
11   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module
12   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
13   !!            3.3  !  2010-11  (G. Madec)  initialisation in C1D configuration
14   !!            3.6  !  2013     ( J. Simeon, C. Calone, G. Madec, C. Ethe ) Online coarsening of outputs
15   !!            3.7  !  2015-11  (G. Madec, A. Coward)  time varying zgr by default
16   !!            4.0  !  2016-10  (G. Madec, S. Flavoni)  domain configuration / user defined interface
17   !!            4.1  !  2020-02  (G. Madec, S. Techene)  introduce ssh to h0 ratio
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   dom_init      : initialize the space and time domain
22   !!   dom_glo       : initialize global domain <--> local domain indices
23   !!   dom_nam       : read and contral domain namelists
24   !!   dom_ctl       : control print for the ocean domain
25   !!   domain_cfg    : read the global domain size in domain configuration file
26   !!   cfg_write     : create the domain configuration file
27   !!----------------------------------------------------------------------
28   USE oce            ! ocean variables
29   USE dom_oce        ! domain: ocean
30   USE domtile        ! tiling utilities
31#if defined key_qco
32   USE domqco         ! quasi-eulerian coord.
33#elif defined key_linssh
34   !                  ! fix in time coord.
35#else
36   USE domvvl         ! variable volume coord.
37#endif
38#if defined key_agrif
39   USE agrif_oce_interp, ONLY : Agrif_istate_ssh ! ssh interpolated from parent
40#endif
41   USE sbc_oce        ! surface boundary condition: ocean
42   USE trc_oce        ! shared ocean & passive tracers variab
43   USE phycst         ! physical constants
44   USE domhgr         ! domain: set the horizontal mesh
45   USE domzgr         ! domain: set the vertical mesh
46   USE dommsk         ! domain: set the mask system
47   USE domwri         ! domain: write the meshmask file
48   USE wet_dry , ONLY : ll_wd     ! wet & drying flag
49   USE closea  , ONLY : dom_clo   ! closed seas routine
50   USE c1d
51   !
52   USE in_out_manager ! I/O manager
53   USE iom            ! I/O library
54   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
55   USE lib_mpp        ! distributed memory computing library
56   USE restart        ! only for lrst_oce and rst_read_ssh
57
58   IMPLICIT NONE
59   PRIVATE
60
61   PUBLIC   dom_init     ! called by nemogcm.F90
62   PUBLIC   domain_cfg   ! called by nemogcm.F90
63
64   !! * Substitutions
65#  include "do_loop_substitute.h90"
66   !!-------------------------------------------------------------------------
67   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
68   !! $Id$
69   !! Software governed by the CeCILL license (see ./LICENSE)
70   !!-------------------------------------------------------------------------
71CONTAINS
72
73   SUBROUTINE dom_init( Kbb, Kmm, Kaa )
74      !!----------------------------------------------------------------------
75      !!                  ***  ROUTINE dom_init  ***
76      !!
77      !! ** Purpose :   Domain initialization. Call the routines that are
78      !!              required to create the arrays which define the space
79      !!              and time domain of the ocean model.
80      !!
81      !! ** Method  : - dom_msk: compute the masks from the bathymetry file
82      !!              - dom_hgr: compute or read the horizontal grid-point position
83      !!                         and scale factors, and the coriolis factor
84      !!              - dom_zgr: define the vertical coordinate and the bathymetry
85      !!              - dom_wri: create the meshmask file (ln_meshmask=T)
86      !!              - 1D configuration, move Coriolis, u and v at T-point
87      !!----------------------------------------------------------------------
88      INTEGER          , INTENT(in) :: Kbb, Kmm, Kaa          ! ocean time level indices
89      !
90      INTEGER ::   ji, jj, jk, jt   ! dummy loop indices
91      INTEGER ::   iconf = 0    ! local integers
92      REAL(wp)::   zrdt
93      CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))"
94      INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level
95      REAL(wp), DIMENSION(jpi,jpj) ::   z1_hu_0, z1_hv_0
96      !!----------------------------------------------------------------------
97      !
98      IF(lwp) THEN         ! Ocean domain Parameters (control print)
99         WRITE(numout,*)
100         WRITE(numout,*) 'dom_init : domain initialization'
101         WRITE(numout,*) '~~~~~~~~'
102         !
103         WRITE(numout,*)     '   Domain info'
104         WRITE(numout,*)     '      dimension of model:'
105         WRITE(numout,*)     '             Local domain      Global domain       Data domain '
106         WRITE(numout,cform) '        ','   jpi     : ', jpi, '   jpiglo  : ', jpiglo
107         WRITE(numout,cform) '        ','   jpj     : ', jpj, '   jpjglo  : ', jpjglo
108         WRITE(numout,cform) '        ','   jpk     : ', jpk, '   jpkglo  : ', jpkglo
109         WRITE(numout,cform) '       ' ,'   jpij    : ', jpij
110         WRITE(numout,*)     '      mpp local domain info (mpp):'
111         WRITE(numout,*)     '              jpni    : ', jpni, '   nn_hls  : ', nn_hls
112         WRITE(numout,*)     '              jpnj    : ', jpnj, '   nn_hls  : ', nn_hls
113         WRITE(numout,*)     '              jpnij   : ', jpnij
114         WRITE(numout,*)     '      lateral boundary of the Global domain:'
115         WRITE(numout,*)     '              cyclic east-west             :', l_Iperio
116         WRITE(numout,*)     '              cyclic north-south           :', l_Jperio
117         WRITE(numout,*)     '              North Pole folding           :', l_NFold
118         WRITE(numout,*)     '                 type of North pole Folding:', c_NFtype
119         WRITE(numout,*)     '      Ocean model configuration used:'
120         WRITE(numout,*)     '              cn_cfg = ', TRIM( cn_cfg ), '   nn_cfg = ', nn_cfg
121      ENDIF
122
123      !
124      !           !==  Reference coordinate system  ==!
125      !
126!      CALL dom_glo                      ! global domain versus local domain
127      CALL dom_nam                      ! read namelist ( namrun, namdom )
128      CALL dom_tile_init                ! Tile domain
129
130      IF( ln_c1d )   CALL     c1d_init                    ! 1D column configuration
131      !
132      CALL dom_hgr                      ! Horizontal mesh
133
134      IF( ln_closea ) CALL dom_clo      ! Read in masks to define closed seas and lakes
135
136      CALL dom_zgr( ik_top, ik_bot )    ! Vertical mesh and bathymetry (return top and bottom ocean t-level indices)
137
138      CALL dom_msk( ik_top, ik_bot )    ! Masks
139      !
140      ht_0(:,:) = 0._wp  ! Reference ocean thickness
141      hu_0(:,:) = 0._wp
142      hv_0(:,:) = 0._wp
143      hf_0(:,:) = 0._wp
144      DO jk = 1, jpkm1
145         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
146         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
147         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
148      END DO
149      !
150      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
151         hf_0(ji,jj) = hf_0(ji,jj) + e3f_0(ji,jj,jk)*vmask(ji,jj,jk)*vmask(ji+1,jj,jk)
152      END_3D
153      CALL lbc_lnk('domain', hf_0, 'F', 1._wp)
154      !
155      IF( lk_SWE ) THEN      ! SWE case redefine hf_0
156         hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,1) * ssfmask(:,:)
157      ENDIF
158      !
159      r1_ht_0(:,:) = ssmask (:,:) / ( ht_0(:,:) + 1._wp -  ssmask (:,:) )
160      r1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp -  ssumask(:,:) )
161      r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp -  ssvmask(:,:) )
162      r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp -  ssfmask(:,:) )
163      !
164      IF( ll_wd ) THEN       ! wet and drying (check ht_0 >= 0)
165         DO_2D( 1, 1, 1, 1 )
166            IF( ht_0(ji,jj) < 0._wp .AND. ssmask(ji,jj) == 1._wp ) THEN
167               CALL ctl_stop( 'dom_init : ht_0 must be positive at potentially wet points' )
168            ENDIF
169         END_2D
170      ENDIF
171      !
172      !           !==  initialisation of time varying coordinate  ==!
173      !
174      !                                 != ssh initialization
175      !
176      IF( l_SAS ) THEN        !* No ocean dynamics calculation : set to 0
177         ssh(:,:,:) = 0._wp
178#if defined key_agrif
179      ELSEIF( .NOT.Agrif_root() .AND.    &
180         &     ln_init_chfrpar ) THEN        !* Interpolate initial ssh from parent
181         CALL Agrif_istate_ssh( Kbb, Kmm, Kaa )
182#endif
183      ELSE                                   !* Read in restart file or set by user
184         CALL rst_read_ssh( Kbb, Kmm, Kaa )
185      ENDIF
186      !     
187#if defined key_qco
188      !                                 != Quasi-Euerian coordinate case
189      !
190      IF( .NOT.l_offline )   CALL dom_qco_init( Kbb, Kmm, Kaa )
191#elif defined key_linssh
192      !                                 != Fix in time : key_linssh case, set through domzgr_substitute.h90
193#else
194      !
195      IF( ln_linssh ) THEN              != Fix in time : set to the reference one for all
196         !
197         DO jt = 1, jpt                         ! depth of t- and w-grid-points
198            gdept(:,:,:,jt) = gdept_0(:,:,:)
199            gdepw(:,:,:,jt) = gdepw_0(:,:,:)
200         END DO
201            gde3w(:,:,:)    = gde3w_0(:,:,:)    ! = gdept as the sum of e3t
202         !
203         DO jt = 1, jpt                         ! vertical scale factors
204            e3t (:,:,:,jt) =  e3t_0(:,:,:)
205            e3u (:,:,:,jt) =  e3u_0(:,:,:)
206            e3v (:,:,:,jt) =  e3v_0(:,:,:)
207            e3w (:,:,:,jt) =  e3w_0(:,:,:)
208            e3uw(:,:,:,jt) = e3uw_0(:,:,:)
209            e3vw(:,:,:,jt) = e3vw_0(:,:,:)
210         END DO
211            e3f (:,:,:)    =  e3f_0(:,:,:)
212         !
213         DO jt = 1, jpt                         ! water column thickness and its inverse
214               hu(:,:,jt) =    hu_0(:,:)
215               hv(:,:,jt) =    hv_0(:,:)
216            r1_hu(:,:,jt) = r1_hu_0(:,:)
217            r1_hv(:,:,jt) = r1_hv_0(:,:)
218         END DO
219               ht   (:,:) =    ht_0(:,:)
220         !
221      ELSE                              != Time varying : initialize before/now/after variables
222         !
223         IF( .NOT.l_offline )   CALL dom_vvl_init( Kbb, Kmm, Kaa )
224         !
225      ENDIF
226#endif
227
228      !
229
230#if defined key_agrif
231      IF( .NOT. Agrif_Root() ) CALL Agrif_Init_Domain( Kbb, Kmm, Kaa )
232#endif
233      IF( ln_meshmask    )   CALL dom_wri       ! Create a domain file
234      IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control
235      !
236      IF( ln_write_cfg   )   CALL cfg_write     ! create the configuration file
237      !
238      IF(lwp) THEN
239         WRITE(numout,*)
240         WRITE(numout,*) 'dom_init :   ==>>>   END of domain initialization'
241         WRITE(numout,*) '~~~~~~~~'
242         WRITE(numout,*)
243      ENDIF
244      !
245   END SUBROUTINE dom_init
246
247
248   SUBROUTINE dom_glo
249      !!----------------------------------------------------------------------
250      !!                     ***  ROUTINE dom_glo  ***
251      !!
252      !! ** Purpose :   initialization of global domain <--> local domain indices
253      !!
254      !! ** Method  :
255      !!
256      !! ** Action  : - mig , mjg : local  domain indices ==> global domain, including halos, indices
257      !!              - mig0, mjg0: local  domain indices ==> global domain, excluding halos, indices
258      !!              - mi0 , mi1 : global domain indices ==> local  domain indices
259      !!              - mj0 , mj1   (if global point not in the local domain ==> mi0>mi1 and/or mj0>mj1)
260      !!----------------------------------------------------------------------
261      INTEGER ::   ji, jj   ! dummy loop argument
262      !!----------------------------------------------------------------------
263      !
264      DO ji = 1, jpi                 ! local domain indices ==> global domain indices, including halos
265        mig(ji) = ji + nimpp - 1
266      END DO
267      DO jj = 1, jpj
268        mjg(jj) = jj + njmpp - 1
269      END DO
270      !                              ! local domain indices ==> global domain indices, excluding halos
271      !
272      mig0(:) = mig(:) - nn_hls
273      mjg0(:) = mjg(:) - nn_hls
274      !                              ! global domain, including halos, indices ==> local domain indices
275      !                                   ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the
276      !                                   ! local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain.
277      DO ji = 1, jpiglo
278        mi0(ji) = MAX( 1 , MIN( ji - nimpp + 1, jpi+1 ) )
279        mi1(ji) = MAX( 0 , MIN( ji - nimpp + 1, jpi   ) )
280      END DO
281      DO jj = 1, jpjglo
282        mj0(jj) = MAX( 1 , MIN( jj - njmpp + 1, jpj+1 ) )
283        mj1(jj) = MAX( 0 , MIN( jj - njmpp + 1, jpj   ) )
284      END DO
285      IF(lwp) THEN                   ! control print
286         WRITE(numout,*)
287         WRITE(numout,*) 'dom_glo : domain: global <<==>> local '
288         WRITE(numout,*) '~~~~~~~ '
289         WRITE(numout,*) '   global domain:   jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpkglo = ', jpkglo
290         WRITE(numout,*) '   local  domain:   jpi    = ', jpi   , ' jpj    = ', jpj   , ' jpk    = ', jpk
291         WRITE(numout,*)
292      ENDIF
293      !
294   END SUBROUTINE dom_glo
295
296
297   SUBROUTINE dom_nam
298      !!----------------------------------------------------------------------
299      !!                     ***  ROUTINE dom_nam  ***
300      !!
301      !! ** Purpose :   read domaine namelists and print the variables.
302      !!
303      !! ** input   : - namrun namelist
304      !!              - namdom namelist
305      !!              - namtile namelist
306      !!              - namnc4 namelist   ! "key_netcdf4" only
307      !!----------------------------------------------------------------------
308      USE ioipsl
309      !!
310      INTEGER ::   ios   ! Local integer
311      REAL(wp)::   zrdt
312      !!----------------------------------------------------------------------
313      !
314      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 &
315         &             nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     &
316         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
317         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, ln_1st_euler  , &
318         &             ln_cfmeta, ln_xios_read, nn_wxios
319      NAMELIST/namdom/ ln_linssh, rn_Dt, rn_atfp, ln_crs, ln_c1d, ln_meshmask
320      NAMELIST/namtile/ ln_tile, nn_ltile_i, nn_ltile_j
321#if defined key_netcdf4
322      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
323#endif
324      !!----------------------------------------------------------------------
325      !
326      IF(lwp) THEN
327         WRITE(numout,*)
328         WRITE(numout,*) 'dom_nam : domain initialization through namelist read'
329         WRITE(numout,*) '~~~~~~~ '
330      ENDIF
331      !
332      !                       !=======================!
333      !                       !==  namelist namdom  ==!
334      !                       !=======================!
335      !
336      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
337903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist' )
338      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
339904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist' )
340      IF(lwm) WRITE( numond, namdom )
341      !
342#if defined key_linssh
343      ln_linssh = lk_linssh      ! overwrite ln_linssh with the logical associated with key_linssh
344#endif
345      !
346#if defined key_agrif
347      IF( .NOT. Agrif_Root() ) THEN    ! AGRIF child, subdivide the Parent timestep
348         rn_Dt = Agrif_Parent (rn_Dt ) / Agrif_Rhot()
349      ENDIF
350#endif
351      !
352      IF(lwp) THEN
353         WRITE(numout,*)
354         WRITE(numout,*) '   Namelist : namdom   ---   space & time domain'
355         WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh
356         WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask
357         WRITE(numout,*) '      ocean time step                         rn_Dt       = ', rn_Dt
358         WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp
359         WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs
360         WRITE(numout,*) '      single column domain (1x1pt)            ln_c1d      = ', ln_c1d
361      ENDIF
362      !
363      ! set current model timestep rDt = 2*rn_Dt if MLF or rDt = rn_Dt if RK3
364      rDt   = 2._wp * rn_Dt
365      r1_Dt = 1._wp / rDt
366      !
367      IF( l_SAS .AND. .NOT.ln_linssh ) THEN
368         CALL ctl_warn( 'SAS requires linear ssh : force ln_linssh = T' )
369         ln_linssh = .TRUE.
370      ENDIF
371      !
372#if defined key_qco
373      IF( ln_linssh )   CALL ctl_stop( 'STOP','domain: key_qco and ln_linssh=T or key_linssh are incompatible' )
374#endif
375      !
376      !                       !=======================!
377      !                       !==  namelist namrun  ==!
378      !                       !=======================!
379      !
380      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
381901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist' )
382      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
383902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist' )
384      IF(lwm) WRITE ( numond, namrun )
385
386#if defined key_agrif
387      IF( .NOT. Agrif_Root() ) THEN
388            nn_it000 = (Agrif_Parent(nn_it000)-1)*Agrif_IRhot() + 1
389            nn_itend =  Agrif_Parent(nn_itend)   *Agrif_IRhot()
390      ENDIF
391#endif
392      !
393      IF(lwp) THEN                  ! control print
394         WRITE(numout,*) '   Namelist : namrun   ---   run parameters'
395         WRITE(numout,*) '      Assimilation cycle              nn_no           = ', nn_no
396         WRITE(numout,*) '      experiment name for output      cn_exp          = ', TRIM( cn_exp           )
397         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in    = ', TRIM( cn_ocerst_in     )
398         WRITE(numout,*) '      restart input directory         cn_ocerst_indir = ', TRIM( cn_ocerst_indir  )
399         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out   = ', TRIM( cn_ocerst_out    )
400         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', TRIM( cn_ocerst_outdir )
401         WRITE(numout,*) '      restart logical                 ln_rstart       = ', ln_rstart
402         WRITE(numout,*) '      start with forward time step    ln_1st_euler    = ', ln_1st_euler
403         WRITE(numout,*) '      control of time step            nn_rstctl       = ', nn_rstctl
404         WRITE(numout,*) '      number of the first time step   nn_it000        = ', nn_it000
405         WRITE(numout,*) '      number of the last time step    nn_itend        = ', nn_itend
406         WRITE(numout,*) '      initial calendar date aammjj    nn_date0        = ', nn_date0
407         WRITE(numout,*) '      initial time of day in hhmm     nn_time0        = ', nn_time0
408         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy        = ', nn_leapy
409         WRITE(numout,*) '      initial state output            nn_istate       = ', nn_istate
410         IF( ln_rst_list ) THEN
411            WRITE(numout,*) '      list of restart dump times      nn_stocklist    =', nn_stocklist
412         ELSE
413            WRITE(numout,*) '      frequency of restart file       nn_stock        = ', nn_stock
414         ENDIF
415#if ! defined key_xios
416         WRITE(numout,*) '      frequency of output file        nn_write        = ', nn_write
417#endif
418         WRITE(numout,*) '      mask land points                ln_mskland      = ', ln_mskland
419         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta       = ', ln_cfmeta
420         WRITE(numout,*) '      overwrite an existing file      ln_clobber      = ', ln_clobber
421         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz      = ', nn_chunksz
422         IF( TRIM(Agrif_CFixed()) == '0' ) THEN
423            WRITE(numout,*) '      READ restart for a single file using XIOS ln_xios_read =', ln_xios_read
424            WRITE(numout,*) '      Write restart using XIOS        nn_wxios   = ', nn_wxios
425         ELSE
426            WRITE(numout,*) "      AGRIF: nn_wxios will be ingored. See setting for parent"
427            WRITE(numout,*) "      AGRIF: ln_xios_read will be ingored. See setting for parent"
428         ENDIF
429      ENDIF
430
431      cexper = cn_exp         ! conversion DOCTOR names into model names (this should disappear soon)
432      nrstdt = nn_rstctl
433      nit000 = nn_it000
434      nitend = nn_itend
435      ndate0 = nn_date0
436      nleapy = nn_leapy
437      ninist = nn_istate
438      !
439      !                                        !==  Set parameters for restart reading using xIOS  ==!
440      !
441      IF( TRIM(Agrif_CFixed()) == '0' ) THEN
442         lrxios = ln_xios_read .AND. ln_rstart
443         IF( nn_wxios > 0 )   lwxios = .TRUE.           !* set output file type for XIOS based on NEMO namelist
444         nxioso = nn_wxios
445      ENDIF
446      !                                        !==  Check consistency between ln_rstart and ln_1st_euler  ==!   (i.e. set l_1st_euler)
447      l_1st_euler = ln_1st_euler
448      !
449      IF( ln_rstart ) THEN                              !*  Restart case
450         !
451         IF(lwp) WRITE(numout,*)
452         IF(lwp) WRITE(numout,*) '   open the restart file'
453         CALL rst_read_open                                              !- Open the restart file
454         !
455         IF( iom_varid( numror, 'rdt', ldstop = .FALSE. ) > 0 ) THEN     !- Check time-step consistency and force Euler restart if changed
456            CALL iom_get( numror, 'rdt', zrdt )
457            IF( zrdt /= rn_Dt ) THEN
458               IF(lwp) WRITE( numout,*)
459               IF(lwp) WRITE( numout,*) '   rn_Dt = ', rn_Dt,' not equal to the READ one rdt = ', zrdt
460               IF(lwp) WRITE( numout,*)
461               IF(lwp) WRITE( numout,*) '      ==>>>   forced euler first time-step'
462               l_1st_euler =  .TRUE.
463            ENDIF
464         ENDIF
465         !
466         IF( .NOT.l_SAS .AND. iom_varid( numror, 'sshb', ldstop = .FALSE. ) <= 0 ) THEN   !- Check absence of one of the Kbb field (here sshb)
467            !                                                                             !  (any Kbb field is missing ==> all Kbb fields are missing)
468            IF( .NOT.l_1st_euler ) THEN
469               CALL ctl_warn('dom_nam : ssh at Kbb not found in restart files ',   &
470                  &                        'l_1st_euler forced to .true. and ' ,   &
471                  &                        'ssh(Kbb) = ssh(Kmm) '                  )
472               l_1st_euler = .TRUE.
473            ENDIF
474         ENDIF
475      ELSEIF( .NOT.l_1st_euler ) THEN                   !*  Initialization case
476         IF(lwp) WRITE(numout,*)
477         IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)'
478         IF(lwp) WRITE(numout,*)'           an Euler initial time step is used : l_1st_euler is forced to .true. '
479         l_1st_euler = .TRUE.
480      ENDIF
481      !
482      !                                        !==  control of output frequency  ==!
483      !
484      IF( .NOT. ln_rst_list ) THEN   ! we use nn_stock
485         IF( nn_stock == -1 )   CALL ctl_warn( 'nn_stock = -1 --> no restart will be done' )
486         IF( nn_stock == 0 .OR. nn_stock > nitend ) THEN
487            WRITE(ctmp1,*) 'nn_stock = ', nn_stock, ' it is forced to ', nitend
488            CALL ctl_warn( ctmp1 )
489            nn_stock = nitend
490         ENDIF
491      ENDIF
492#if ! defined key_xios
493      IF( nn_write == -1 )   CALL ctl_warn( 'nn_write = -1 --> no output files will be done' )
494      IF ( nn_write == 0 ) THEN
495         WRITE(ctmp1,*) 'nn_write = ', nn_write, ' it is forced to ', nitend
496         CALL ctl_warn( ctmp1 )
497         nn_write = nitend
498      ENDIF
499#endif
500
501      IF( Agrif_Root() ) THEN
502         IF(lwp) WRITE(numout,*)
503         SELECT CASE ( nleapy )                !==  Choose calendar for IOIPSL  ==!
504         CASE (  1 )
505            CALL ioconf_calendar('gregorian')
506            IF(lwp) WRITE(numout,*) '   ==>>>   The IOIPSL calendar is "gregorian", i.e. leap year'
507         CASE (  0 )
508            CALL ioconf_calendar('noleap')
509            IF(lwp) WRITE(numout,*) '   ==>>>   The IOIPSL calendar is "noleap", i.e. no leap year'
510         CASE ( 30 )
511            CALL ioconf_calendar('360d')
512            IF(lwp) WRITE(numout,*) '   ==>>>   The IOIPSL calendar is "360d", i.e. 360 days in a year'
513         END SELECT
514      ENDIF
515      !
516      !                       !========================!
517      !                       !==  namelist namtile  ==!
518      !                       !========================!
519      !
520      READ  ( numnam_ref, namtile, IOSTAT = ios, ERR = 905 )
521905   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtile in reference namelist' )
522      READ  ( numnam_cfg, namtile, IOSTAT = ios, ERR = 906 )
523906   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtile in configuration namelist' )
524      IF(lwm) WRITE( numond, namtile )
525
526      IF(lwp) THEN
527         WRITE(numout,*)
528         WRITE(numout,*)    '   Namelist : namtile   ---   Domain tiling decomposition'
529         WRITE(numout,*)    '      Tiling (T) or not (F)                ln_tile    = ', ln_tile
530         WRITE(numout,*)    '      Length of tile in i                  nn_ltile_i = ', nn_ltile_i
531         WRITE(numout,*)    '      Length of tile in j                  nn_ltile_j = ', nn_ltile_j
532         WRITE(numout,*)
533         IF( ln_tile ) THEN
534            WRITE(numout,*) '      The domain will be decomposed into tiles of size', nn_ltile_i, 'x', nn_ltile_j
535         ELSE
536            WRITE(numout,*) '      Domain tiling will NOT be used'
537         ENDIF
538      ENDIF
539      !
540#if defined key_netcdf4
541      !                       !=======================!
542      !                       !==  namelist namnc4  ==!   NetCDF 4 case   ("key_netcdf4" defined)
543      !                       !=======================!
544      !
545      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
546907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namnc4 in reference namelist' )
547      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
548908   IF( ios >  0 )   CALL ctl_nam ( ios , 'namnc4 in configuration namelist' )
549      IF(lwm) WRITE( numond, namnc4 )
550
551      IF(lwp) THEN                        ! control print
552         WRITE(numout,*)
553         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters ("key_netcdf4" defined)'
554         WRITE(numout,*) '      number of chunks in i-dimension             nn_nchunks_i = ', nn_nchunks_i
555         WRITE(numout,*) '      number of chunks in j-dimension             nn_nchunks_j = ', nn_nchunks_j
556         WRITE(numout,*) '      number of chunks in k-dimension             nn_nchunks_k = ', nn_nchunks_k
557         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression   ln_nc4zip    = ', ln_nc4zip
558      ENDIF
559
560      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
561      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
562      snc4set%ni   = nn_nchunks_i
563      snc4set%nj   = nn_nchunks_j
564      snc4set%nk   = nn_nchunks_k
565      snc4set%luse = ln_nc4zip
566#else
567      snc4set%luse = .FALSE.        ! No NetCDF 4 case
568#endif
569      !
570   END SUBROUTINE dom_nam
571
572
573   SUBROUTINE dom_ctl
574      !!----------------------------------------------------------------------
575      !!                     ***  ROUTINE dom_ctl  ***
576      !!
577      !! ** Purpose :   Domain control.
578      !!
579      !! ** Method  :   compute and print extrema of masked scale factors
580      !!----------------------------------------------------------------------
581      LOGICAL, DIMENSION(jpi,jpj) ::   llmsk
582      INTEGER, DIMENSION(2)       ::   imil, imip, imi1, imi2, imal, imap, ima1, ima2
583      REAL(wp)                    ::   zglmin, zglmax, zgpmin, zgpmax, ze1min, ze1max, ze2min, ze2max
584      !!----------------------------------------------------------------------
585      !
586      llmsk = tmask_h(:,:) == 1._wp
587      !
588      CALL mpp_minloc( 'domain', glamt(:,:), llmsk, zglmin, imil )
589      CALL mpp_minloc( 'domain', gphit(:,:), llmsk, zgpmin, imip )
590      CALL mpp_minloc( 'domain',   e1t(:,:), llmsk, ze1min, imi1 )
591      CALL mpp_minloc( 'domain',   e2t(:,:), llmsk, ze2min, imi2 )
592      CALL mpp_maxloc( 'domain', glamt(:,:), llmsk, zglmax, imal )
593      CALL mpp_maxloc( 'domain', gphit(:,:), llmsk, zgpmax, imap )
594      CALL mpp_maxloc( 'domain',   e1t(:,:), llmsk, ze1max, ima1 )
595      CALL mpp_maxloc( 'domain',   e2t(:,:), llmsk, ze2max, ima2 )
596      !
597      IF(lwp) THEN
598         WRITE(numout,*)
599         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
600         WRITE(numout,*) '~~~~~~~'
601         WRITE(numout,"(14x,'glamt mini: ',1f10.2,' at i = ',i5,' j= ',i5)") zglmin, imil(1), imil(2)
602         WRITE(numout,"(14x,'glamt maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") zglmax, imal(1), imal(2)
603         WRITE(numout,"(14x,'gphit mini: ',1f10.2,' at i = ',i5,' j= ',i5)") zgpmin, imip(1), imip(2)
604         WRITE(numout,"(14x,'gphit maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") zgpmax, imap(1), imap(2)
605         WRITE(numout,"(14x,'  e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, imi1(1), imi1(2)
606         WRITE(numout,"(14x,'  e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, ima1(1), ima1(2)
607         WRITE(numout,"(14x,'  e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, imi2(1), imi2(2)
608         WRITE(numout,"(14x,'  e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, ima2(1), ima2(2)
609      ENDIF
610      !
611   END SUBROUTINE dom_ctl
612
613
614   SUBROUTINE domain_cfg( cd_cfg, kk_cfg, kpi, kpj, kpk, ldIperio, ldJperio, ldNFold, cdNFtype )
615      !!----------------------------------------------------------------------
616      !!                     ***  ROUTINE domain_cfg  ***
617      !!
618      !! ** Purpose :   read the domain size in domain configuration file
619      !!
620      !! ** Method  :   read the cn_domcfg NetCDF file
621      !!----------------------------------------------------------------------
622      CHARACTER(len=*), INTENT(out) ::   cd_cfg               ! configuration name
623      INTEGER         , INTENT(out) ::   kk_cfg               ! configuration resolution
624      INTEGER         , INTENT(out) ::   kpi, kpj, kpk        ! global domain sizes
625      LOGICAL         , INTENT(out) ::   ldIperio, ldJperio   ! i- and j- periodicity
626      LOGICAL         , INTENT(out) ::   ldNFold              ! North pole folding
627      CHARACTER(len=1), INTENT(out) ::   cdNFtype             ! Folding type: T or F
628      !
629      CHARACTER(len=7) ::   catt                  ! 'T', 'F', '-' or 'UNKNOWN'
630      INTEGER ::   inum, iperio, iatt             ! local integer
631      REAL(wp) ::   zorca_res                     ! local scalars
632      REAL(wp) ::   zperio                        !   -      -
633      INTEGER, DIMENSION(4) ::   idvar, idimsz    ! size   of dimensions
634      !!----------------------------------------------------------------------
635      !
636      IF(lwp) THEN
637         WRITE(numout,*) '           '
638         WRITE(numout,*) 'domain_cfg : domain size read in ', TRIM( cn_domcfg ), ' file'
639         WRITE(numout,*) '~~~~~~~~~~ '
640      ENDIF
641      !
642      CALL iom_open( cn_domcfg, inum )
643      !
644      CALL iom_getatt( inum,  'CfgName', cd_cfg )   ! returns 'UNKNOWN' if not found
645      CALL iom_getatt( inum, 'CfgIndex', kk_cfg )   ! returns      -999 if not found
646      !
647      ! ------- keep compatibility with OLD VERSION... start -------
648      IF( cd_cfg == 'UNKNOWN' .AND. kk_cfg == -999 ) THEN
649         IF(  iom_varid( inum, 'ORCA'       , ldstop = .FALSE. ) > 0  .AND.  &
650            & iom_varid( inum, 'ORCA_index' , ldstop = .FALSE. ) > 0    ) THEN
651            !
652            cd_cfg = 'ORCA'
653            CALL iom_get( inum, 'ORCA_index', zorca_res )   ;   kk_cfg = NINT( zorca_res )
654            !
655         ELSE
656            CALL iom_getatt( inum, 'cn_cfg', cd_cfg )  ! returns 'UNKNOWN' if not found
657            CALL iom_getatt( inum, 'nn_cfg', kk_cfg )  ! returns      -999 if not found
658         ENDIF
659      ENDIF
660      ! ------- keep compatibility with OLD VERSION... end -------
661      !
662      idvar = iom_varid( inum, 'e3t_0', kdimsz = idimsz )   ! use e3t_0, that must exist, to get jp(ijk)glo
663      kpi = idimsz(1)
664      kpj = idimsz(2)
665      kpk = idimsz(3)
666      !
667      CALL iom_getatt( inum, 'Iperio', iatt )   ;   ldIperio = iatt == 1   ! returns      -999 if not found -> default = .false.
668      CALL iom_getatt( inum, 'Jperio', iatt )   ;   ldJperio = iatt == 1   ! returns      -999 if not found -> default = .false.
669      CALL iom_getatt( inum,  'NFold', iatt )   ;   ldNFold  = iatt == 1   ! returns      -999 if not found -> default = .false.
670      CALL iom_getatt( inum, 'NFtype', catt )                              ! returns 'UNKNOWN' if not found
671      IF( LEN_TRIM(catt) == 1 ) THEN   ;   cdNFtype = TRIM(catt)
672      ELSE                             ;   cdNFtype = '-'
673      ENDIF
674      !
675      ! ------- keep compatibility with OLD VERSION... start -------
676      IF( iatt == -999 .AND. catt == 'UNKNOWN' .AND. iom_varid( inum, 'jperio', ldstop = .FALSE. ) > 0 ) THEN                   
677         CALL iom_get( inum, 'jperio', zperio )   ;   iperio = NINT( zperio )
678         ldIperio = iperio == 1  .OR. iperio == 4 .OR. iperio == 6 .OR. iperio == 7   ! i-periodicity
679         ldJperio = iperio == 2  .OR. iperio == 7                                     ! j-periodicity
680         ldNFold  = iperio >= 3 .AND. iperio <= 6                                     ! North pole folding
681         IF(     iperio == 3 .OR. iperio == 4 ) THEN   ;   cdNFtype = 'T'             !    folding at T point
682         ELSEIF( iperio == 5 .OR. iperio == 6 ) THEN   ;   cdNFtype = 'F'             !    folding at F point
683         ELSE                                          ;   cdNFtype = '-'             !    default value
684         ENDIF
685      ENDIF
686      ! ------- keep compatibility with OLD VERSION... end -------
687      !
688      CALL iom_close( inum )
689      !
690      IF(lwp) THEN
691         WRITE(numout,*) '   .'
692         WRITE(numout,*) '   ==>>>   ', TRIM(cn_cfg), ' configuration '
693         WRITE(numout,*) '   .'
694         WRITE(numout,*) '      nn_cfg = ', kk_cfg
695         WRITE(numout,*) '      Ni0glo = ', kpi
696         WRITE(numout,*) '      Nj0glo = ', kpj
697         WRITE(numout,*) '      jpkglo = ', kpk
698      ENDIF
699      !
700   END SUBROUTINE domain_cfg
701
702
703   SUBROUTINE cfg_write
704      !!----------------------------------------------------------------------
705      !!                  ***  ROUTINE cfg_write  ***
706      !!
707      !! ** Purpose :   Create the "cn_domcfg_out" file, a NetCDF file which
708      !!              contains all the ocean domain informations required to
709      !!              define an ocean configuration.
710      !!
711      !! ** Method  :   Write in a file all the arrays required to set up an
712      !!              ocean configuration.
713      !!
714      !! ** output file :   domcfg_out.nc : domain size, characteristics, horizontal
715      !!                       mesh, Coriolis parameter, and vertical scale factors
716      !!                    NB: also contain ORCA family information
717      !!----------------------------------------------------------------------
718      INTEGER           ::   ji, jj, jk   ! dummy loop indices
719      INTEGER           ::   inum     ! local units
720      CHARACTER(len=21) ::   clnam    ! filename (mesh and mask informations)
721      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! workspace
722      !!----------------------------------------------------------------------
723      !
724      IF(lwp) WRITE(numout,*)
725      IF(lwp) WRITE(numout,*) 'cfg_write : create the domain configuration file (', TRIM(cn_domcfg_out),'.nc)'
726      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
727      !
728      !                       ! ============================= !
729      !                       !  create 'domcfg_out.nc' file  !
730      !                       ! ============================= !
731      !
732      clnam = cn_domcfg_out  ! filename (configuration information)
733      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE. )
734      !
735      !                             !==  Configuration specificities  ==!
736      !
737      CALL iom_putatt( inum,  'CfgName', TRIM(cn_cfg) )
738      CALL iom_putatt( inum, 'CfgIndex',      nn_cfg  )
739      !
740      !                             !==  domain characteristics  ==!
741      !
742      !                                   ! lateral boundary of the global domain
743      CALL iom_putatt( inum, 'Iperio', COUNT( (/l_Iperio/) ) )
744      CALL iom_putatt( inum, 'Jperio', COUNT( (/l_Jperio/) ) )
745      CALL iom_putatt( inum,  'NFold', COUNT( (/l_NFold /) ) )
746      CALL iom_putatt( inum, 'NFtype',          c_NFtype     )
747
748      !                                   ! type of vertical coordinate
749      IF(ln_zco)   CALL iom_putatt( inum, 'VertCoord', 'zco' )
750      IF(ln_zps)   CALL iom_putatt( inum, 'VertCoord', 'zps' )
751      IF(ln_sco)   CALL iom_putatt( inum, 'VertCoord', 'sco' )
752     
753      !                                   ! ocean cavities under iceshelves
754      CALL iom_putatt( inum, 'IsfCav', COUNT( (/ln_isfcav/) ) )
755      !
756      !                             !==  horizontal mesh  !
757      !
758      CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 )   ! latitude
759      CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
760      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
761      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
762      !
763      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )   ! longitude
764      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
765      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
766      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
767      !
768      CALL iom_rstput( 0, 0, inum, 'e1t'  , e1t  , ktype = jp_r8 )   ! i-scale factors (e1.)
769      CALL iom_rstput( 0, 0, inum, 'e1u'  , e1u  , ktype = jp_r8 )
770      CALL iom_rstput( 0, 0, inum, 'e1v'  , e1v  , ktype = jp_r8 )
771      CALL iom_rstput( 0, 0, inum, 'e1f'  , e1f  , ktype = jp_r8 )
772      !
773      CALL iom_rstput( 0, 0, inum, 'e2t'  , e2t  , ktype = jp_r8 )   ! j-scale factors (e2.)
774      CALL iom_rstput( 0, 0, inum, 'e2u'  , e2u  , ktype = jp_r8 )
775      CALL iom_rstput( 0, 0, inum, 'e2v'  , e2v  , ktype = jp_r8 )
776      CALL iom_rstput( 0, 0, inum, 'e2f'  , e2f  , ktype = jp_r8 )
777      !
778      CALL iom_rstput( 0, 0, inum, 'ff_f' , ff_f , ktype = jp_r8 )   ! coriolis factor
779      CALL iom_rstput( 0, 0, inum, 'ff_t' , ff_t , ktype = jp_r8 )
780      !
781      !                             !==  vertical mesh  ==!
782      !
783      CALL iom_rstput( 0, 0, inum, 'e3t_1d'  , e3t_1d , ktype = jp_r8 )   ! reference 1D-coordinate
784      CALL iom_rstput( 0, 0, inum, 'e3w_1d'  , e3w_1d , ktype = jp_r8 )
785      !
786      CALL iom_rstput( 0, 0, inum, 'e3t_0'   , e3t_0  , ktype = jp_r8 )   ! vertical scale factors
787      CALL iom_rstput( 0, 0, inum, 'e3u_0'   , e3u_0  , ktype = jp_r8 )
788      CALL iom_rstput( 0, 0, inum, 'e3v_0'   , e3v_0  , ktype = jp_r8 )
789      CALL iom_rstput( 0, 0, inum, 'e3f_0'   , e3f_0  , ktype = jp_r8 )
790      CALL iom_rstput( 0, 0, inum, 'e3w_0'   , e3w_0  , ktype = jp_r8 )
791      CALL iom_rstput( 0, 0, inum, 'e3uw_0'  , e3uw_0 , ktype = jp_r8 )
792      CALL iom_rstput( 0, 0, inum, 'e3vw_0'  , e3vw_0 , ktype = jp_r8 )
793      !
794      !                             !==  wet top and bottom level  ==!   (caution: multiplied by ssmask)
795      !
796      CALL iom_rstput( 0, 0, inum, 'top_level'    , REAL( mikt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points (ISF)
797      CALL iom_rstput( 0, 0, inum, 'bottom_level' , REAL( mbkt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points
798      !
799      IF( ln_sco ) THEN             ! s-coordinate: store grid stiffness ratio  (Not required anyway)
800         CALL dom_stiff( z2d )
801         CALL iom_rstput( 0, 0, inum, 'stiffness', z2d )        !    ! Max. grid stiffness ratio
802      ENDIF
803      !
804      IF( ll_wd ) THEN              ! wetting and drying domain
805         CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 )
806      ENDIF
807      !                       ! ============================ !
808      !                       !        close the files
809      !                       ! ============================ !
810      CALL iom_close( inum )
811      !
812   END SUBROUTINE cfg_write
813
814   !!======================================================================
815END MODULE domain
Note: See TracBrowser for help on using the repository browser.