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

source: branches/UKMO/dev_r5107_restart_date/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 5494

Last change on this file since 5494 was 5494, checked in by dancopsey, 9 years ago

Upgraded NEMO3,5 branch dev/hadco/vn3.5_beta_restart_date to this revision.

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