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

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 4152

Last change on this file since 4152 was 4152, checked in by cetlod, 10 years ago

merge in dev_LOCEAN_2013 the 2nd development branch dev_r3940_CNRS4_IOCRS, see ticket #1169
  • Property svn:keywords set to Id
File size: 23.0 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      ! 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(:,:,jk) * umask(:,:,jk)
97         hv(:,:) = hv(:,:) + fse3v(:,:,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_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,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
231         WRITE(numout,*) '      min number of ocean level (<0)       '
232         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
233         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
234         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
235         WRITE(numout,*) '           = 0   no file created           '
236         WRITE(numout,*) '           = 1   mesh_mask                 '
237         WRITE(numout,*) '           = 2   mesh and mask             '
238         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
239         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
240         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
241         WRITE(numout,*) '      time-splitting: nb of sub time-step   nn_baro   = ', nn_baro
242         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
243         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
244         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
245         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
246         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
247         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
248         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
249         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
250         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
251         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
252         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
253         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
254         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
255         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
256         WRITE(numout,*) '                                        ppa0            = ', ppa0
257         WRITE(numout,*) '                                        ppa1            = ', ppa1
258         WRITE(numout,*) '                                        ppkth           = ', ppkth
259         WRITE(numout,*) '                                        ppacr           = ', ppacr
260         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
261         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
262         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
263         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
264         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
265         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
266      ENDIF
267
268      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
269      e3zps_min = rn_e3zps_min
270      e3zps_rat = rn_e3zps_rat
271      nmsh      = nn_msh
272      nacc      = nn_acc
273      atfp      = rn_atfp
274      rdt       = rn_rdt
275      rdtmin    = rn_rdtmin
276      rdtmax    = rn_rdtmin
277      rdth      = rn_rdth
278
279      REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection
280      READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)
281905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )
282
283      REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection
284      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
285906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
286      WRITE( numond, namcla )
287
288      IF(lwp) THEN
289         WRITE(numout,*)
290         WRITE(numout,*) '   Namelist namcla'
291         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
292      ENDIF
293      IF ( nn_cla .EQ. 1 ) THEN
294         IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2
295            CONTINUE
296         ELSE
297            CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' )
298         ENDIF
299      ENDIF
300
301#if defined key_netcdf4
302      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
303      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
304      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
305907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
306
307      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
308      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
309908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
310      WRITE( numond, namnc4 )
311
312      IF(lwp) THEN                        ! control print
313         WRITE(numout,*)
314         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
315         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
316         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
317         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
318         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
319      ENDIF
320
321      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
322      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
323      snc4set%ni   = nn_nchunks_i
324      snc4set%nj   = nn_nchunks_j
325      snc4set%nk   = nn_nchunks_k
326      snc4set%luse = ln_nc4zip
327#else
328      snc4set%luse = .FALSE.        ! No NetCDF 4 case
329#endif
330      !
331   END SUBROUTINE dom_nam
332
333
334   SUBROUTINE dom_ctl
335      !!----------------------------------------------------------------------
336      !!                     ***  ROUTINE dom_ctl  ***
337      !!
338      !! ** Purpose :   Domain control.
339      !!
340      !! ** Method  :   compute and print extrema of masked scale factors
341      !!----------------------------------------------------------------------
342      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
343      INTEGER, DIMENSION(2) ::   iloc   !
344      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
345      !!----------------------------------------------------------------------
346      !
347      IF(lk_mpp) THEN
348         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
349         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
350         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
351         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
352      ELSE
353         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )   
354         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )   
355         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )   
356         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )   
357
358         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
359         iimi1 = iloc(1) + nimpp - 1
360         ijmi1 = iloc(2) + njmpp - 1
361         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
362         iimi2 = iloc(1) + nimpp - 1
363         ijmi2 = iloc(2) + njmpp - 1
364         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
365         iima1 = iloc(1) + nimpp - 1
366         ijma1 = iloc(2) + njmpp - 1
367         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
368         iima2 = iloc(1) + nimpp - 1
369         ijma2 = iloc(2) + njmpp - 1
370      ENDIF
371      IF(lwp) THEN
372         WRITE(numout,*)
373         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
374         WRITE(numout,*) '~~~~~~~'
375         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
376         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
377         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
378         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
379      ENDIF
380      !
381   END SUBROUTINE dom_ctl
382
383   SUBROUTINE dom_stiff
384      !!----------------------------------------------------------------------
385      !!                  ***  ROUTINE dom_stiff  ***
386      !!                     
387      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
388      !!
389      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
390      !!                Save the maximum in the vertical direction
391      !!                (this number is only relevant in s-coordinates)
392      !!
393      !!                Haney, R. L., 1991: On the pressure gradient force
394      !!                over steep topography in sigma coordinate ocean models.
395      !!                J. Phys. Oceanogr., 21, 610???619.
396      !!----------------------------------------------------------------------
397      INTEGER  ::   ji, jj, jk 
398      REAL(wp) ::   zrxmax
399      REAL(wp), DIMENSION(4) :: zr1
400      !!----------------------------------------------------------------------
401      rx1(:,:) = 0.e0
402      zrxmax   = 0.e0
403      zr1(:)   = 0.e0
404     
405      DO ji = 2, jpim1
406         DO jj = 2, jpjm1
407            DO jk = 1, jpkm1
408               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji-1,jj  ,jk  )  & 
409                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1)) &
410                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji-1,jj  ,jk  )  &
411                    &                         -gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1) + rsmall) )
412               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw(ji+1,jj  ,jk  )-gdepw(ji  ,jj  ,jk  )  &
413                    &                         +gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
414                    &                        /(gdepw(ji+1,jj  ,jk  )+gdepw(ji  ,jj  ,jk  )  &
415                    &                         -gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
416               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw(ji  ,jj+1,jk  )-gdepw(ji  ,jj  ,jk  )  &
417                    &                         +gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
418                    &                        /(gdepw(ji  ,jj+1,jk  )+gdepw(ji  ,jj  ,jk  )  &
419                    &                         -gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
420               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji  ,jj-1,jk  )  &
421                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1)) &
422                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji  ,jj-1,jk  )  &
423                    &                         -gdepw(ji,  jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1) + rsmall) )
424               zrxmax = MAXVAL(zr1(1:4))
425               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
426            END DO
427         END DO
428      END DO
429
430      CALL lbc_lnk( rx1, 'T', 1. )
431
432      zrxmax = MAXVAL(rx1)
433
434      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
435
436      IF(lwp) THEN
437         WRITE(numout,*)
438         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
439         WRITE(numout,*) '~~~~~~~~~'
440      ENDIF
441
442   END SUBROUTINE dom_stiff
443
444
445
446   !!======================================================================
447END MODULE domain
Note: See TracBrowser for help on using the repository browser.