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/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 5066

Last change on this file since 5066 was 5066, checked in by rfurner, 9 years ago

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