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

source: branches/UKMO/r5936_hadgem3_mct/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 7127

Last change on this file since 7127 was 7127, checked in by jcastill, 7 years ago

Remove svn keywords

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