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

Last change on this file was 13558, checked in by smasson, 2 months ago

trunk: pass sette tests with debugging option -init=arrays,snan,huge, see #2535

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