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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 4366

Last change on this file since 4366 was 4366, checked in by jchanut, 10 years ago

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