source: NEMO/trunk/src/SWE/domain.F90 @ 13295

Last change on this file since 13295 was 13295, checked in by acc, 10 months ago

Replace do-loop macros in the trunk with alternative forms with greater flexibility for extra halo applications. This alters a lot of routines but does not change any behaviour or results. do_loop_substitute.h90 is greatly simplified by this change. SETTE results are identical to those with the previous revision

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