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

source: branches/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 4044

Last change on this file since 4044 was 4044, checked in by clevy, 11 years ago

Configuration setting/add SETTE compatibility, see ticket:#1074

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