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/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 5862

Last change on this file since 5862 was 5862, checked in by gm, 8 years ago

#1613: vvl by default: non-vvl: initialize _b,n,a scale factors with _0 arrays

  • Property svn:keywords set to Id
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   !!             -   ! 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 argument
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) * tmask(:,:,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      IF( lk_vvl ) THEN       !==  time varying coordinate system  ==!
104         !
105         CALL dom_vvl_init                      ! set before/now/after variables
106         !
107      ELSE                    !==  Fix in time  ==!   set everything to the reference one for all
108         !
109         !       before        !          now          !       after         !
110         ;  gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
111         ;  gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
112         ;                     ;   gde3w_n = gde3w_0   !        ---          !
113         !                                                                 
114         ;    e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
115         ;    e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
116         ;    e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
117         ;                     ;     e3f_n =   e3f_0   !        ---          !
118         ;    e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
119         ;   e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
120         ;   e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
121         !
122         CALL wrk_alloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
123         !
124         z1_hu_0(:,:) = umask_i(:,:) / ( hu_0(:,:) + 1._wp - umask_i(:,:) )     ! _i mask due to ISF
125         z1_hv_0(:,:) = vmask_i(:,:) / ( hv_0(:,:) + 1._wp - vmask_i(:,:) )
126         !
127         !        before       !          now          !       after         !
128         ;                     ;      ht_n =    ht_0   !                     ! water column thickness
129         ;     hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   !
130         ;     hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   !
131         ;  r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness
132         ;  r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   !
133         !
134         CALL wrk_dealloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
135      ENDIF
136      !
137      IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
138      !
139                             CALL dom_stp       ! time step
140      IF( nmsh /= 0      )   CALL dom_wri       ! Create a domain file
141      IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control
142      !
143      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
144      !
145   END SUBROUTINE dom_init
146
147
148   SUBROUTINE dom_nam
149      !!----------------------------------------------------------------------
150      !!                     ***  ROUTINE dom_nam  ***
151      !!                   
152      !! ** Purpose :   read domaine namelists and print the variables.
153      !!
154      !! ** input   : - namrun namelist
155      !!              - namdom namelist
156      !!              - namnc4 namelist   ! "key_netcdf4" only
157      !!----------------------------------------------------------------------
158      USE ioipsl
159      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
160         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
161         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
162         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler
163      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   &
164         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  &
165         &             rn_rdtmax, rn_rdth     , nn_closea , ln_crs,    &
166         &             jphgr_msh, &
167         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
168         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
169         &             ppa2, ppkth2, ppacr2
170#if defined key_netcdf4
171      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
172#endif
173      INTEGER  ::   ios                 ! Local integer output status for namelist read
174      !!----------------------------------------------------------------------
175
176      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
177      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
178901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
179
180      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
181      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
182902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
183      IF(lwm) WRITE ( numond, namrun )
184      !
185      IF(lwp) THEN                  ! control print
186         WRITE(numout,*)
187         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
188         WRITE(numout,*) '~~~~~~~ '
189         WRITE(numout,*) '   Namelist namrun'
190         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
191         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
192         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
193         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
194         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
195         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
196         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
197         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
198         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
199         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
200         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
201         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
202         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
203         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
204         IF( ln_rst_list ) THEN
205            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
206         ELSE
207            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
208         ENDIF
209         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
210         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
211         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
212         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
213         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
214         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
215      ENDIF
216
217      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
218      cexper = cn_exp
219      nrstdt = nn_rstctl
220      nit000 = nn_it000
221      nitend = nn_itend
222      ndate0 = nn_date0
223      nleapy = nn_leapy
224      ninist = nn_istate
225      nstock = nn_stock
226      nstocklist = nn_stocklist
227      nwrite = nn_write
228      neuler = nn_euler
229      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
230         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
231         CALL ctl_warn( ctmp1 )
232         neuler = 0
233      ENDIF
234
235      !                             ! control of output frequency
236      IF ( nstock == 0 .OR. nstock > nitend ) THEN
237         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
238         CALL ctl_warn( ctmp1 )
239         nstock = nitend
240      ENDIF
241      IF ( nwrite == 0 ) THEN
242         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
243         CALL ctl_warn( ctmp1 )
244         nwrite = nitend
245      ENDIF
246
247#if defined key_agrif
248      IF( Agrif_Root() ) THEN
249#endif
250      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
251      CASE (  1 ) 
252         CALL ioconf_calendar('gregorian')
253         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
254      CASE (  0 )
255         CALL ioconf_calendar('noleap')
256         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
257      CASE ( 30 )
258         CALL ioconf_calendar('360d')
259         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
260      END SELECT
261#if defined key_agrif
262      ENDIF
263#endif
264
265      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
266      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
267903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
268 
269      !
270      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
271      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
272904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
273      IF(lwm) WRITE ( numond, namdom )
274      !
275      IF(lwp) THEN
276         WRITE(numout,*)
277         WRITE(numout,*) '   Namelist namdom : space & time domain'
278         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
279         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
280         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
281         WRITE(numout,*) '      min number of ocean level (<0)       '
282         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
283         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
284         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
285         WRITE(numout,*) '           = 0   no file created           '
286         WRITE(numout,*) '           = 1   mesh_mask                 '
287         WRITE(numout,*) '           = 2   mesh and mask             '
288         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
289         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
290         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
291         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
292         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
293         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
294         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
295         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
296         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
297         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
298         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
299         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
300         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
301         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
302         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
303         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
304         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
305         WRITE(numout,*) '                                        ppa0            = ', ppa0
306         WRITE(numout,*) '                                        ppa1            = ', ppa1
307         WRITE(numout,*) '                                        ppkth           = ', ppkth
308         WRITE(numout,*) '                                        ppacr           = ', ppacr
309         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
310         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
311         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
312         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
313         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
314         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
315      ENDIF
316      !
317      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
318      e3zps_min = rn_e3zps_min
319      e3zps_rat = rn_e3zps_rat
320      nmsh      = nn_msh
321      nacc      = nn_acc
322      atfp      = rn_atfp
323      rdt       = rn_rdt
324      rdtmin    = rn_rdtmin
325      rdtmax    = rn_rdtmin
326      rdth      = rn_rdth
327
328#if defined key_netcdf4
329      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
330      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
331      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
332907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
333
334      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
335      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
336908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
337      IF(lwm) WRITE( numond, namnc4 )
338
339      IF(lwp) THEN                        ! control print
340         WRITE(numout,*)
341         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
342         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
343         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
344         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
345         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
346      ENDIF
347
348      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
349      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
350      snc4set%ni   = nn_nchunks_i
351      snc4set%nj   = nn_nchunks_j
352      snc4set%nk   = nn_nchunks_k
353      snc4set%luse = ln_nc4zip
354#else
355      snc4set%luse = .FALSE.        ! No NetCDF 4 case
356#endif
357      !
358   END SUBROUTINE dom_nam
359
360
361   SUBROUTINE dom_ctl
362      !!----------------------------------------------------------------------
363      !!                     ***  ROUTINE dom_ctl  ***
364      !!
365      !! ** Purpose :   Domain control.
366      !!
367      !! ** Method  :   compute and print extrema of masked scale factors
368      !!----------------------------------------------------------------------
369      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
370      INTEGER, DIMENSION(2) ::   iloc   !
371      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
372      !!----------------------------------------------------------------------
373      !
374      IF(lk_mpp) THEN
375         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
376         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
377         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
378         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
379      ELSE
380         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
381         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
382         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
383         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
384
385         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
386         iimi1 = iloc(1) + nimpp - 1
387         ijmi1 = iloc(2) + njmpp - 1
388         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
389         iimi2 = iloc(1) + nimpp - 1
390         ijmi2 = iloc(2) + njmpp - 1
391         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
392         iima1 = iloc(1) + nimpp - 1
393         ijma1 = iloc(2) + njmpp - 1
394         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
395         iima2 = iloc(1) + nimpp - 1
396         ijma2 = iloc(2) + njmpp - 1
397      ENDIF
398      IF(lwp) THEN
399         WRITE(numout,*)
400         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
401         WRITE(numout,*) '~~~~~~~'
402         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
403         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
404         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
405         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
406      ENDIF
407      !
408   END SUBROUTINE dom_ctl
409
410
411   SUBROUTINE dom_stiff
412      !!----------------------------------------------------------------------
413      !!                  ***  ROUTINE dom_stiff  ***
414      !!                     
415      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
416      !!
417      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
418      !!                Save the maximum in the vertical direction
419      !!                (this number is only relevant in s-coordinates)
420      !!
421      !!                Haney, R. L., 1991: On the pressure gradient force
422      !!                over steep topography in sigma coordinate ocean models.
423      !!                J. Phys. Oceanogr., 21, 610???619.
424      !!----------------------------------------------------------------------
425      INTEGER  ::   ji, jj, jk 
426      REAL(wp) ::   zrxmax
427      REAL(wp), DIMENSION(4) :: zr1
428      !!----------------------------------------------------------------------
429      rx1(:,:) = 0._wp
430      zrxmax   = 0._wp
431      zr1(:)   = 0._wp
432      !
433      DO ji = 2, jpim1
434         DO jj = 2, jpjm1
435            DO jk = 1, jpkm1
436               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
437                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
438                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
439                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
440               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
441                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
442                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
443                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
444               zr1(3) =ABS(  (  gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
445                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
446                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
447                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
448               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
449                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
450                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
451                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
452               zrxmax = MAXVAL( zr1(1:4) )
453               rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax )
454            END DO
455         END DO
456      END DO
457      CALL lbc_lnk( rx1, 'T', 1. )
458      !
459      zrxmax = MAXVAL( rx1 )
460      !
461      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
462      !
463      IF(lwp) THEN
464         WRITE(numout,*)
465         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
466         WRITE(numout,*) '~~~~~~~~~'
467      ENDIF
468      !
469   END SUBROUTINE dom_stiff
470
471   !!======================================================================
472END MODULE domain
Note: See TracBrowser for help on using the repository browser.