source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/penzAM98/MY_SRC/domain.F90 @ 13562

Last change on this file since 13562 was 13562, checked in by gm, 5 months ago

zgr_pse created only for NS coast

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