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

Last change on this file since 5747 was 5747, checked in by rfurner, 5 years ago

bug fix to ln_rstdate in domain.F90

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