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/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domain.F90 @ 13895

Last change on this file since 13895 was 13895, checked in by techene, 4 years ago

#2574 => led to restart reorganization and some cleaning

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