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/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM – NEMO

source: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domain.F90 @ 15006

Last change on this file since 15006 was 15006, checked in by gsamson, 3 years ago

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