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

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

Time split cleaning / AMM12 restart / BDY and Agrif. See tickets #1228, #1227, #1225 and #1133

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