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 @ 5883

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

#1613: vvl by default: TRA/TRC remove optimization associated with linear free surface

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