New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domain.F90 in NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DOM – NEMO

source: NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DOM/domain.F90 @ 14789

Last change on this file since 14789 was 14789, checked in by mcastril, 3 years ago

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