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

source: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 5770

Last change on this file since 5770 was 5770, checked in by gm, 9 years ago

#1593: LDF-ADV, step II.2: phasing the improvements/simplifications of advective tracer trend (see wiki page)
  • Property svn:keywords set to Id
File size: 23.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   !!----------------------------------------------------------------------
16   
17   !!----------------------------------------------------------------------
18   !!   dom_init       : initialize the space and time domain
19   !!   dom_nam        : read and contral domain namelists
20   !!   dom_ctl        : control print for the ocean domain
21   !!----------------------------------------------------------------------
22   USE oce             ! ocean variables
23   USE dom_oce         ! domain: ocean
24   USE sbc_oce         ! surface boundary condition: ocean
25   USE phycst          ! physical constants
26   USE closea          ! closed seas
27   USE in_out_manager  ! I/O manager
28   USE lib_mpp         ! distributed memory computing library
29
30   USE domhgr          ! domain: set the horizontal mesh
31   USE domzgr          ! domain: set the vertical mesh
32   USE domstp          ! domain: set the time-step
33   USE dommsk          ! domain: set the mask system
34   USE domwri          ! domain: write the meshmask file
35   USE domvvl          ! variable volume
36   USE c1d             ! 1D vertical configuration
37   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
38   USE timing          ! Timing
39   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   dom_init   ! called by opa.F90
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48   !!-------------------------------------------------------------------------
49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
50   !! $Id$
51   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
52   !!-------------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE dom_init
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE dom_init  ***
58      !!                   
59      !! ** Purpose :   Domain initialization. Call the routines that are
60      !!              required to create the arrays which define the space
61      !!              and time domain of the ocean model.
62      !!
63      !! ** Method  : - dom_msk: compute the masks from the bathymetry file
64      !!              - dom_hgr: compute or read the horizontal grid-point position
65      !!                         and scale factors, and the coriolis factor
66      !!              - dom_zgr: define the vertical coordinate and the bathymetry
67      !!              - dom_stp: defined the model time step
68      !!              - dom_wri: create the meshmask file if nmsh=1
69      !!              - 1D configuration, move Coriolis, u and v at T-point
70      !!----------------------------------------------------------------------
71      INTEGER ::   jk          ! dummy loop argument
72      INTEGER ::   iconf = 0   ! local integers
73      !!----------------------------------------------------------------------
74      !
75      IF( nn_timing == 1 )   CALL timing_start('dom_init')
76      !
77      IF(lwp) THEN
78         WRITE(numout,*)
79         WRITE(numout,*) 'dom_init : domain initialization'
80         WRITE(numout,*) '~~~~~~~~'
81      ENDIF
82      !
83                             CALL dom_nam      ! read namelist ( namrun, namdom )
84                             CALL dom_clo      ! Closed seas and lake
85                             CALL dom_hgr      ! Horizontal mesh
86                             CALL dom_zgr      ! Vertical mesh and bathymetry
87                             CALL dom_msk      ! Masks
88      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
89      !
90      ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at T-points
91      hu_0(:,:) = 0.0_wp                       ! Reference ocean depth at U-points
92      hv_0(:,:) = 0.0_wp                       ! Reference ocean depth at V-points
93      DO jk = 1, jpk
94         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
95         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
96         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
97      END DO
98      !
99      IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh
100      !
101      IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
102      !
103      !
104      hu(:,:) = 0._wp                          ! Ocean depth at U-points
105      hv(:,:) = 0._wp                          ! Ocean depth at V-points
106      ht(:,:) = 0._wp                          ! Ocean depth at T-points
107      DO jk = 1, jpkm1
108         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
109         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
110         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
111      END DO
112      !                                        ! Inverse of the local depth
113      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
114      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
115
116                             CALL dom_stp      ! time step
117      IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file
118      IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control
119      !
120      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
121      !
122   END SUBROUTINE dom_init
123
124
125   SUBROUTINE dom_nam
126      !!----------------------------------------------------------------------
127      !!                     ***  ROUTINE dom_nam  ***
128      !!                   
129      !! ** Purpose :   read domaine namelists and print the variables.
130      !!
131      !! ** input   : - namrun namelist
132      !!              - namdom namelist
133      !!              - namnc4 namelist   ! "key_netcdf4" only
134      !!----------------------------------------------------------------------
135      USE ioipsl
136      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
137         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
138         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
139         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler
140      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   &
141         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  &
142         &             rn_rdtmax, rn_rdth     , nn_closea , ln_crs,    &
143         &             jphgr_msh, &
144         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
145         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
146         &             ppa2, ppkth2, ppacr2
147#if defined key_netcdf4
148      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
149#endif
150      INTEGER  ::   ios                 ! Local integer output status for namelist read
151      !!----------------------------------------------------------------------
152
153      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
154      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
155901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
156
157      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
158      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
159902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
160      IF(lwm) WRITE ( numond, namrun )
161      !
162      IF(lwp) THEN                  ! control print
163         WRITE(numout,*)
164         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
165         WRITE(numout,*) '~~~~~~~ '
166         WRITE(numout,*) '   Namelist namrun'
167         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
168         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
169         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
170         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
171         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
172         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
173         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
174         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
175         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
176         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
177         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
178         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
179         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
180         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
181         IF( ln_rst_list ) THEN
182            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
183         ELSE
184            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
185         ENDIF
186         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
187         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
188         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
189         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
190         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
191         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
192      ENDIF
193
194      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
195      cexper = cn_exp
196      nrstdt = nn_rstctl
197      nit000 = nn_it000
198      nitend = nn_itend
199      ndate0 = nn_date0
200      nleapy = nn_leapy
201      ninist = nn_istate
202      nstock = nn_stock
203      nstocklist = nn_stocklist
204      nwrite = nn_write
205      neuler = nn_euler
206      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
207         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
208         CALL ctl_warn( ctmp1 )
209         neuler = 0
210      ENDIF
211
212      !                             ! control of output frequency
213      IF ( nstock == 0 .OR. nstock > nitend ) THEN
214         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
215         CALL ctl_warn( ctmp1 )
216         nstock = nitend
217      ENDIF
218      IF ( nwrite == 0 ) THEN
219         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
220         CALL ctl_warn( ctmp1 )
221         nwrite = nitend
222      ENDIF
223
224#if defined key_agrif
225      IF( Agrif_Root() ) THEN
226#endif
227      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
228      CASE (  1 ) 
229         CALL ioconf_calendar('gregorian')
230         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
231      CASE (  0 )
232         CALL ioconf_calendar('noleap')
233         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
234      CASE ( 30 )
235         CALL ioconf_calendar('360d')
236         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
237      END SELECT
238#if defined key_agrif
239      ENDIF
240#endif
241
242      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
243      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
244903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
245 
246      !
247      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
248      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
249904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
250      IF(lwm) WRITE ( numond, namdom )
251
252      IF(lwp) THEN
253         WRITE(numout,*)
254         WRITE(numout,*) '   Namelist namdom : space & time domain'
255         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
256         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
257         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
258         WRITE(numout,*) '      min number of ocean level (<0)       '
259         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
260         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
261         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
262         WRITE(numout,*) '           = 0   no file created           '
263         WRITE(numout,*) '           = 1   mesh_mask                 '
264         WRITE(numout,*) '           = 2   mesh and mask             '
265         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
266         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
267         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
268         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
269         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
270         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
271         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
272         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
273         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
274         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
275         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
276         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
277         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
278         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
279         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
280         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
281         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
282         WRITE(numout,*) '                                        ppa0            = ', ppa0
283         WRITE(numout,*) '                                        ppa1            = ', ppa1
284         WRITE(numout,*) '                                        ppkth           = ', ppkth
285         WRITE(numout,*) '                                        ppacr           = ', ppacr
286         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
287         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
288         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
289         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
290         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
291         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
292      ENDIF
293
294      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
295      e3zps_min = rn_e3zps_min
296      e3zps_rat = rn_e3zps_rat
297      nmsh      = nn_msh
298      nacc      = nn_acc
299      atfp      = rn_atfp
300      rdt       = rn_rdt
301      rdtmin    = rn_rdtmin
302      rdtmax    = rn_rdtmin
303      rdth      = rn_rdth
304
305#if defined key_netcdf4
306      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
307      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
308      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
309907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
310
311      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
312      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
313908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
314      IF(lwm) WRITE( numond, namnc4 )
315
316      IF(lwp) THEN                        ! control print
317         WRITE(numout,*)
318         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
319         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
320         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
321         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
322         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
323      ENDIF
324
325      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
326      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
327      snc4set%ni   = nn_nchunks_i
328      snc4set%nj   = nn_nchunks_j
329      snc4set%nk   = nn_nchunks_k
330      snc4set%luse = ln_nc4zip
331#else
332      snc4set%luse = .FALSE.        ! No NetCDF 4 case
333#endif
334      !
335   END SUBROUTINE dom_nam
336
337
338   SUBROUTINE dom_ctl
339      !!----------------------------------------------------------------------
340      !!                     ***  ROUTINE dom_ctl  ***
341      !!
342      !! ** Purpose :   Domain control.
343      !!
344      !! ** Method  :   compute and print extrema of masked scale factors
345      !!----------------------------------------------------------------------
346      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
347      INTEGER, DIMENSION(2) ::   iloc   !
348      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
349      !!----------------------------------------------------------------------
350      !
351      IF(lk_mpp) THEN
352         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
353         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
354         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
355         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
356      ELSE
357         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
358         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
359         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
360         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
361
362         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
363         iimi1 = iloc(1) + nimpp - 1
364         ijmi1 = iloc(2) + njmpp - 1
365         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
366         iimi2 = iloc(1) + nimpp - 1
367         ijmi2 = iloc(2) + njmpp - 1
368         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
369         iima1 = iloc(1) + nimpp - 1
370         ijma1 = iloc(2) + njmpp - 1
371         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
372         iima2 = iloc(1) + nimpp - 1
373         ijma2 = iloc(2) + njmpp - 1
374      ENDIF
375      IF(lwp) THEN
376         WRITE(numout,*)
377         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
378         WRITE(numout,*) '~~~~~~~'
379         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
380         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
381         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
382         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
383      ENDIF
384      !
385   END SUBROUTINE dom_ctl
386
387   SUBROUTINE dom_stiff
388      !!----------------------------------------------------------------------
389      !!                  ***  ROUTINE dom_stiff  ***
390      !!                     
391      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
392      !!
393      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
394      !!                Save the maximum in the vertical direction
395      !!                (this number is only relevant in s-coordinates)
396      !!
397      !!                Haney, R. L., 1991: On the pressure gradient force
398      !!                over steep topography in sigma coordinate ocean models.
399      !!                J. Phys. Oceanogr., 21, 610???619.
400      !!----------------------------------------------------------------------
401      INTEGER  ::   ji, jj, jk 
402      REAL(wp) ::   zrxmax
403      REAL(wp), DIMENSION(4) :: zr1
404      !!----------------------------------------------------------------------
405      rx1(:,:) = 0.e0
406      zrxmax   = 0.e0
407      zr1(:)   = 0.e0
408     
409      DO ji = 2, jpim1
410         DO jj = 2, jpjm1
411            DO jk = 1, jpkm1
412               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji-1,jj  ,jk  )  & 
413                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1)) &
414                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji-1,jj  ,jk  )  &
415                    &                         -gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1) + rsmall) )
416               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw_0(ji+1,jj  ,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
417                    &                         +gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
418                    &                        /(gdepw_0(ji+1,jj  ,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
419                    &                         -gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
420               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw_0(ji  ,jj+1,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
421                    &                         +gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
422                    &                        /(gdepw_0(ji  ,jj+1,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
423                    &                         -gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
424               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji  ,jj-1,jk  )  &
425                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1)) &
426                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji  ,jj-1,jk  )  &
427                    &                         -gdepw_0(ji,  jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1) + rsmall) )
428               zrxmax = MAXVAL(zr1(1:4))
429               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
430            END DO
431         END DO
432      END DO
433
434      CALL lbc_lnk( rx1, 'T', 1. )
435
436      zrxmax = MAXVAL(rx1)
437
438      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
439
440      IF(lwp) THEN
441         WRITE(numout,*)
442         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
443         WRITE(numout,*) '~~~~~~~~~'
444      ENDIF
445
446   END SUBROUTINE dom_stiff
447
448
449
450   !!======================================================================
451END MODULE domain
Note: See TracBrowser for help on using the repository browser.