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 branches/UKMO/dev_r6393_CO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_r6393_CO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 7023

Last change on this file since 7023 was 7023, checked in by deazer, 7 years ago

Include update for restart date labelling for rose cycling.

File size: 25.2 KB
Line 
1MODULE domain
2   !!==============================================================================
3   !!                       ***  MODULE domain   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6   !! History :  OPA  !  1990-10  (C. Levy - G. Madec)  Original code
7   !!                 !  1992-01  (M. Imbard) insert time step initialization
8   !!                 !  1996-06  (G. Madec) generalized vertical coordinate
9   !!                 !  1997-02  (G. Madec) creation of domwri.F
10   !!                 !  2001-05  (E.Durand - G. Madec) insert closed sea
11   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module
12   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
13   !!            3.3  !  2010-11  (G. Madec)  initialisation in C1D configuration
14   !!            3.6  !  2013     ( J. Simeon, C. Calone, G. Madec, C. Ethe ) Online coarsening of outputs
15   !!            3.7  !  2015-11  (G. Madec, A. Coward)  time varying zgr by default
16   !!----------------------------------------------------------------------
17   
18   !!----------------------------------------------------------------------
19   !!   dom_init       : initialize the space and time domain
20   !!   dom_nam        : read and contral domain namelists
21   !!   dom_ctl        : control print for the ocean domain
22   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean variables
25   USE dom_oce         ! domain: ocean
26   USE sbc_oce         ! surface boundary condition: ocean
27   USE phycst          ! physical constants
28   USE closea          ! closed seas
29   USE domhgr          ! domain: set the horizontal mesh
30   USE domzgr          ! domain: set the vertical mesh
31   USE domstp          ! domain: set the time-step
32   USE dommsk          ! domain: set the mask system
33   USE domwri          ! domain: write the meshmask file
34   USE domvvl          ! variable volume
35   USE c1d             ! 1D vertical configuration
36   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
37   !
38   USE in_out_manager  ! I/O manager
39   USE wrk_nemo        ! Memory Allocation
40   USE lib_mpp         ! distributed memory computing library
41   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
42   USE timing          ! Timing
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   dom_init   ! called by opa.F90
48
49   !!-------------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
53   !!-------------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE dom_init
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE dom_init  ***
59      !!                   
60      !! ** Purpose :   Domain initialization. Call the routines that are
61      !!              required to create the arrays which define the space
62      !!              and time domain of the ocean model.
63      !!
64      !! ** Method  : - dom_msk: compute the masks from the bathymetry file
65      !!              - dom_hgr: compute or read the horizontal grid-point position
66      !!                         and scale factors, and the coriolis factor
67      !!              - dom_zgr: define the vertical coordinate and the bathymetry
68      !!              - dom_stp: defined the model time step
69      !!              - dom_wri: create the meshmask file if nmsh=1
70      !!              - 1D configuration, move Coriolis, u and v at T-point
71      !!----------------------------------------------------------------------
72      INTEGER ::   jk          ! dummy loop indices
73      INTEGER ::   iconf = 0   ! local integers
74      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_hu_0, z1_hv_0
75      !!----------------------------------------------------------------------
76      !
77      IF( nn_timing == 1 )   CALL timing_start('dom_init')
78      !
79      IF(lwp) THEN
80         WRITE(numout,*)
81         WRITE(numout,*) 'dom_init : domain initialization'
82         WRITE(numout,*) '~~~~~~~~'
83      ENDIF
84      !
85      !                       !==  Reference coordinate system  ==!
86      !
87                     CALL dom_nam               ! read namelist ( namrun, namdom )
88                     CALL dom_clo               ! Closed seas and lake
89                     CALL dom_hgr               ! Horizontal mesh
90                     CALL dom_zgr               ! Vertical mesh and bathymetry
91                     CALL dom_msk               ! Masks
92      IF( ln_sco )   CALL dom_stiff             ! Maximum stiffness ratio/hydrostatic consistency
93      !
94      ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1)   ! Reference ocean thickness
95      hu_0(:,:) = e3u_0(:,:,1) * umask(:,:,1)
96      hv_0(:,:) = e3v_0(:,:,1) * vmask(:,:,1)
97      DO jk = 2, jpk
98         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
99         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
100         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
101      END DO
102      !
103      !              !==  time varying part of coordinate system  ==!
104      !
105      IF( ln_linssh ) THEN          ! Fix in time : set to the reference one for all
106         !       before        !          now          !       after         !
107         ;  gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
108         ;  gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
109         ;                     ;   gde3w_n = gde3w_0   !        ---          !
110         !                                                                 
111         ;    e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
112         ;    e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
113         ;    e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
114         ;                     ;     e3f_n =   e3f_0   !        ---          !
115         ;    e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
116         ;   e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
117         ;   e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
118         !
119         CALL wrk_alloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
120         !
121         z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
122         z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
123         !
124         !        before       !          now          !       after         !
125         ;                     ;      ht_n =    ht_0   !                     ! water column thickness
126         ;     hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   !
127         ;     hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   !
128         ;  r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness
129         ;  r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   !
130         !
131         CALL wrk_dealloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
132         !
133      ELSE                         ! time varying : initialize before/now/after variables
134         !
135         CALL dom_vvl_init 
136         !
137      ENDIF
138      !
139      IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
140      !
141                             CALL dom_stp       ! time step
142      IF( nmsh /= 0 .AND. .NOT. ln_iscpl )                         CALL dom_wri      ! Create a domain file
143      IF( nmsh /= 0 .AND.       ln_iscpl .AND. .NOT. ln_rstart )   CALL dom_wri      ! Create a domain file
144      IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control
145      !
146      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
147      !
148   END SUBROUTINE dom_init
149
150
151   SUBROUTINE dom_nam
152      !!----------------------------------------------------------------------
153      !!                     ***  ROUTINE dom_nam  ***
154      !!                   
155      !! ** Purpose :   read domaine namelists and print the variables.
156      !!
157      !! ** input   : - namrun namelist
158      !!              - namdom namelist
159      !!              - namnc4 namelist   ! "key_netcdf4" only
160      !!----------------------------------------------------------------------
161      USE ioipsl
162      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 &
163         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , ln_rstdate, nn_rstctl,   &
164         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
165         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     &
166         &             ln_cfmeta, ln_iscpl
167      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, &
168         &             rn_atfp , rn_rdt   , nn_closea   , ln_crs      , jphgr_msh ,                  &
169         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         &
170         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  &
171         &             ppa2, ppkth2, ppacr2
172#if defined key_netcdf4
173      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
174#endif
175      INTEGER  ::   ios                 ! Local integer output status for namelist read
176      !!----------------------------------------------------------------------
177
178      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
179      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
180901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
181
182      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
183      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
184902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
185      IF(lwm) WRITE ( numond, namrun )
186      !
187      IF(lwp) THEN                  ! control print
188         WRITE(numout,*)
189         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
190         WRITE(numout,*) '~~~~~~~ '
191         WRITE(numout,*) '   Namelist namrun'
192         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
193         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
194         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
195         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
196         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
197         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
198         WRITE(numout,*) '      restart logical                 ln_rstart  = ' , ln_rstart
199         WRITE(numout,*) '      datestamping of restarts        ln_rstdate  = ', ln_rstdate
200         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
201         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
202         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
203         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
204         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
205         WRITE(numout,*) '      initial time of day in hhmm     nn_time0   = ', nn_time0
206         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
207         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
208         IF( ln_rst_list ) THEN
209            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
210         ELSE
211            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
212         ENDIF
213         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
214         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
215         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
216         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
217         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
218         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl   = ', ln_iscpl
219      ENDIF
220
221      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
222      cexper = cn_exp
223      nrstdt = nn_rstctl
224      nit000 = nn_it000
225      nitend = nn_itend
226      ndate0 = nn_date0
227      nleapy = nn_leapy
228      ninist = nn_istate
229      nstock = nn_stock
230      nstocklist = nn_stocklist
231      nwrite = nn_write
232      neuler = nn_euler
233      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
234         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
235         CALL ctl_warn( ctmp1 )
236         neuler = 0
237      ENDIF
238
239      !                             ! control of output frequency
240      IF ( nstock == 0 .OR. nstock > nitend ) THEN
241         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
242         CALL ctl_warn( ctmp1 )
243         nstock = nitend
244      ENDIF
245      IF ( nwrite == 0 ) THEN
246         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
247         CALL ctl_warn( ctmp1 )
248         nwrite = nitend
249      ENDIF
250
251#if defined key_agrif
252      IF( Agrif_Root() ) THEN
253#endif
254      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
255      CASE (  1 ) 
256         CALL ioconf_calendar('gregorian')
257         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
258      CASE (  0 )
259         CALL ioconf_calendar('noleap')
260         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
261      CASE ( 30 )
262         CALL ioconf_calendar('360d')
263         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
264      END SELECT
265#if defined key_agrif
266      ENDIF
267#endif
268
269      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
270      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
271903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
272 
273      !
274      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
275      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
276904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
277      IF(lwm) WRITE ( numond, namdom )
278      !
279      IF(lwp) THEN
280         WRITE(numout,*)
281         WRITE(numout,*) '   Namelist namdom : space & time domain'
282         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
283         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
284         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
285         WRITE(numout,*) '      min number of ocean level (<0)       '
286         WRITE(numout,*) '      treshold to open the isf cavity   rn_isfhmin   = ', rn_isfhmin, ' (m)'
287         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
288         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
289         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
290         WRITE(numout,*) '           = 0   no file created           '
291         WRITE(numout,*) '           = 1   mesh_mask                 '
292         WRITE(numout,*) '           = 2   mesh and mask             '
293         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
294         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
295         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
296         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
297         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
298         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
299         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
300         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
301         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
302         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
303         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
304         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
305         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
306         WRITE(numout,*) '                                        ppa0            = ', ppa0
307         WRITE(numout,*) '                                        ppa1            = ', ppa1
308         WRITE(numout,*) '                                        ppkth           = ', ppkth
309         WRITE(numout,*) '                                        ppacr           = ', ppacr
310         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
311         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
312         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
313         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
314         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
315         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
316      ENDIF
317      !
318      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
319      e3zps_min = rn_e3zps_min
320      e3zps_rat = rn_e3zps_rat
321      nmsh      = nn_msh
322      atfp      = rn_atfp
323      rdt       = rn_rdt
324
325#if defined key_netcdf4
326      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
327      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
328      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
329907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
330
331      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
332      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
333908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
334      IF(lwm) WRITE( numond, namnc4 )
335
336      IF(lwp) THEN                        ! control print
337         WRITE(numout,*)
338         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
339         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
340         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
341         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
342         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
343      ENDIF
344
345      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
346      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
347      snc4set%ni   = nn_nchunks_i
348      snc4set%nj   = nn_nchunks_j
349      snc4set%nk   = nn_nchunks_k
350      snc4set%luse = ln_nc4zip
351#else
352      snc4set%luse = .FALSE.        ! No NetCDF 4 case
353#endif
354      !
355   END SUBROUTINE dom_nam
356
357
358   SUBROUTINE dom_ctl
359      !!----------------------------------------------------------------------
360      !!                     ***  ROUTINE dom_ctl  ***
361      !!
362      !! ** Purpose :   Domain control.
363      !!
364      !! ** Method  :   compute and print extrema of masked scale factors
365      !!----------------------------------------------------------------------
366      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
367      INTEGER, DIMENSION(2) ::   iloc   !
368      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
369      !!----------------------------------------------------------------------
370      !
371      IF(lk_mpp) THEN
372         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
373         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
374         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
375         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
376      ELSE
377         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
378         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
379         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
380         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
381
382         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
383         iimi1 = iloc(1) + nimpp - 1
384         ijmi1 = iloc(2) + njmpp - 1
385         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
386         iimi2 = iloc(1) + nimpp - 1
387         ijmi2 = iloc(2) + njmpp - 1
388         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
389         iima1 = iloc(1) + nimpp - 1
390         ijma1 = iloc(2) + njmpp - 1
391         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
392         iima2 = iloc(1) + nimpp - 1
393         ijma2 = iloc(2) + njmpp - 1
394      ENDIF
395      IF(lwp) THEN
396         WRITE(numout,*)
397         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
398         WRITE(numout,*) '~~~~~~~'
399         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
400         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
401         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
402         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
403      ENDIF
404      !
405   END SUBROUTINE dom_ctl
406
407
408   SUBROUTINE dom_stiff
409      !!----------------------------------------------------------------------
410      !!                  ***  ROUTINE dom_stiff  ***
411      !!                     
412      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
413      !!
414      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
415      !!                Save the maximum in the vertical direction
416      !!                (this number is only relevant in s-coordinates)
417      !!
418      !!                Haney, R. L., 1991: On the pressure gradient force
419      !!                over steep topography in sigma coordinate ocean models.
420      !!                J. Phys. Oceanogr., 21, 610???619.
421      !!----------------------------------------------------------------------
422      INTEGER  ::   ji, jj, jk 
423      REAL(wp) ::   zrxmax
424      REAL(wp), DIMENSION(4) ::   zr1
425      !!----------------------------------------------------------------------
426      rx1(:,:) = 0._wp
427      zrxmax   = 0._wp
428      zr1(:)   = 0._wp
429      !
430      DO ji = 2, jpim1
431         DO jj = 2, jpjm1
432            DO jk = 1, jpkm1
433               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
434                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
435                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
436                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
437               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
438                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
439                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
440                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
441               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
442                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
443                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
444                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
445               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
446                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
447                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
448                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
449               zrxmax = MAXVAL( zr1(1:4) )
450               rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax )
451            END DO
452         END DO
453      END DO
454      CALL lbc_lnk( rx1, 'T', 1. )
455      !
456      zrxmax = MAXVAL( rx1 )
457      !
458      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
459      !
460      IF(lwp) THEN
461         WRITE(numout,*)
462         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
463         WRITE(numout,*) '~~~~~~~~~'
464      ENDIF
465      !
466   END SUBROUTINE dom_stiff
467
468   !!======================================================================
469END MODULE domain
Note: See TracBrowser for help on using the repository browser.