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/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domain.F90 @ 14018

Last change on this file since 14018 was 14018, checked in by techene, 3 years ago

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