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 @ 3993

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

Configuration setting/Step3 bugfixes,doc, and redistribute variables in namcfg and namdom see ticket:#1074

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