source: NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DOM/domain.F90 @ 12326

Last change on this file since 12326 was 12326, checked in by deazer, 17 months ago

Adjusted to account for flux form advection
should now work with both types of advection and enveloping bathymetry
moved one off caluculations of depth scales to domain.F90

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