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/trunk/src/OCE/DOM – NEMO

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

Last change on this file since 14171 was 14171, checked in by jchanut, 3 years ago

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