source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90 @ 10009

Last change on this file since 10009 was 10009, checked in by gm, 3 years ago

#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh

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