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

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 6006

Last change on this file since 6006 was 6006, checked in by mathiot, 8 years ago

Merged ice sheet coupling branch
  • Property svn:keywords set to Id
File size: 24.4 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   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
22   !!----------------------------------------------------------------------
23   USE oce             ! ocean variables
24   USE dom_oce         ! domain: ocean
25   USE sbc_oce         ! surface boundary condition: ocean
26   USE phycst          ! physical constants
27   USE closea          ! closed seas
28   USE domhgr          ! domain: set the horizontal mesh
29   USE domzgr          ! domain: set the vertical mesh
30   USE domstp          ! domain: set the time-step
31   USE dommsk          ! domain: set the mask system
32   USE domwri          ! domain: write the meshmask file
33   USE domvvl          ! variable volume
34   USE c1d             ! 1D vertical configuration
35   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
36   !
37   USE in_out_manager  ! I/O manager
38   USE lib_mpp         ! distributed memory computing library
39   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
40   USE timing          ! Timing
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   dom_init   ! called by opa.F90
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49   !!-------------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
53   !!-------------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE dom_init
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE dom_init  ***
59      !!                   
60      !! ** Purpose :   Domain initialization. Call the routines that are
61      !!              required to create the arrays which define the space
62      !!              and time domain of the ocean model.
63      !!
64      !! ** Method  : - dom_msk: compute the masks from the bathymetry file
65      !!              - dom_hgr: compute or read the horizontal grid-point position
66      !!                         and scale factors, and the coriolis factor
67      !!              - dom_zgr: define the vertical coordinate and the bathymetry
68      !!              - dom_stp: defined the model time step
69      !!              - dom_wri: create the meshmask file if nmsh=1
70      !!              - 1D configuration, move Coriolis, u and v at T-point
71      !!----------------------------------------------------------------------
72      INTEGER ::   jk          ! dummy loop argument
73      INTEGER ::   iconf = 0   ! local integers
74      !!----------------------------------------------------------------------
75      !
76      IF( nn_timing == 1 )   CALL timing_start('dom_init')
77      !
78      IF(lwp) THEN
79         WRITE(numout,*)
80         WRITE(numout,*) 'dom_init : domain initialization'
81         WRITE(numout,*) '~~~~~~~~'
82      ENDIF
83      !
84                             CALL dom_nam      ! read namelist ( namrun, namdom )
85                             CALL dom_clo      ! Closed seas and lake
86                             CALL dom_hgr      ! Horizontal mesh
87                             CALL dom_zgr      ! Vertical mesh and bathymetry
88                             CALL dom_msk      ! Masks
89      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
90      !
91      ht_0(:,:) = 0._wp                        ! Reference ocean depth at T-points
92      hu_0(:,:) = 0._wp                        ! Reference ocean depth at U-points
93      hv_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 - ssumask(:,:) ) * ssumask(:,:)
115      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
116
117                             CALL dom_stp      ! time step
118      !
119      IF( nmsh /= 0 .AND. .NOT. ln_iscpl )                         CALL dom_wri      ! Create a domain file
120      IF( nmsh /= 0 .AND.       ln_iscpl .AND. .NOT. ln_rstart )   CALL dom_wri      ! Create a domain file
121      !
122      IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control
123      !
124      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
125      !
126   END SUBROUTINE dom_init
127
128
129   SUBROUTINE dom_nam
130      !!----------------------------------------------------------------------
131      !!                     ***  ROUTINE dom_nam  ***
132      !!                   
133      !! ** Purpose :   read domaine namelists and print the variables.
134      !!
135      !! ** input   : - namrun namelist
136      !!              - namdom namelist
137      !!              - namnc4 namelist   ! "key_netcdf4" only
138      !!----------------------------------------------------------------------
139      USE ioipsl
140      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 &
141                       nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     &
142         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
143         &             nn_stock, nn_write , ln_dimgnnn  , ln_mskland   , ln_clobber, nn_chunksz,     &
144         &             nn_euler, ln_cfmeta, ln_iscpl
145      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, &
146         &             nn_acc   , rn_atfp , rn_rdt      , rn_rdtmin   ,                              &
147         &             rn_rdtmax, rn_rdth , nn_closea   , ln_crs      ,                              &
148         &             jphgr_msh,                                                                    &
149         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         &
150         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  &
151         &             ppa2, ppkth2, ppacr2
152#if defined key_netcdf4
153      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
154#endif
155      INTEGER  ::   ios                 ! Local integer output status for namelist read
156      !!----------------------------------------------------------------------
157
158      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
159      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
160901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
161
162      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
163      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
164902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
165      IF(lwm) WRITE ( numond, namrun )
166      !
167      IF(lwp) THEN                  ! control print
168         WRITE(numout,*)
169         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
170         WRITE(numout,*) '~~~~~~~ '
171         WRITE(numout,*) '   Namelist namrun'
172         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
173         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
174         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
175         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
176         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
177         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
178         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
179         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
180         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
181         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
182         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
183         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
184         WRITE(numout,*) '      initial time of day in hhmm     nn_time0   = ', nn_time0
185         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
186         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
187         IF( ln_rst_list ) THEN
188            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
189         ELSE
190            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
191         ENDIF
192         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
193         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
194         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
195         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
196         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
197         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
198         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl   = ', ln_iscpl
199      ENDIF
200
201      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
202      cexper = cn_exp
203      nrstdt = nn_rstctl
204      nit000 = nn_it000
205      nitend = nn_itend
206      ndate0 = nn_date0
207      nleapy = nn_leapy
208      ninist = nn_istate
209      nstock = nn_stock
210      nstocklist = nn_stocklist
211      nwrite = nn_write
212      neuler = nn_euler
213      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
214         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
215         CALL ctl_warn( ctmp1 )
216         neuler = 0
217      ENDIF
218
219      !                             ! control of output frequency
220      IF ( nstock == 0 .OR. nstock > nitend ) THEN
221         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
222         CALL ctl_warn( ctmp1 )
223         nstock = nitend
224      ENDIF
225      IF ( nwrite == 0 ) THEN
226         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
227         CALL ctl_warn( ctmp1 )
228         nwrite = nitend
229      ENDIF
230
231#if defined key_agrif
232      IF( Agrif_Root() ) THEN
233#endif
234      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
235      CASE (  1 ) 
236         CALL ioconf_calendar('gregorian')
237         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
238      CASE (  0 )
239         CALL ioconf_calendar('noleap')
240         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
241      CASE ( 30 )
242         CALL ioconf_calendar('360d')
243         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
244      END SELECT
245#if defined key_agrif
246      ENDIF
247#endif
248
249      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
250      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
251903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
252 
253      !
254      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
255      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
256904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
257      IF(lwm) WRITE ( numond, namdom )
258      !
259      IF(lwp) THEN
260         WRITE(numout,*)
261         WRITE(numout,*) '   Namelist namdom : space & time domain'
262         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
263         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
264         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
265         WRITE(numout,*) '      min number of ocean level (<0)       '
266         WRITE(numout,*) '      treshold to open the isf cavity   rn_isfhmin   = ', rn_isfhmin, ' (m)'
267         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
268         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
269         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
270         WRITE(numout,*) '           = 0   no file created           '
271         WRITE(numout,*) '           = 1   mesh_mask                 '
272         WRITE(numout,*) '           = 2   mesh and mask             '
273         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
274         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
275         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
276         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
277         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
278         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
279         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
280         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
281         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
282         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
283         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
284         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
285         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
286         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
287         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
288         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
289         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
290         WRITE(numout,*) '                                        ppa0            = ', ppa0
291         WRITE(numout,*) '                                        ppa1            = ', ppa1
292         WRITE(numout,*) '                                        ppkth           = ', ppkth
293         WRITE(numout,*) '                                        ppacr           = ', ppacr
294         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
295         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
296         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
297         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
298         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
299         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
300      ENDIF
301      !
302      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
303      e3zps_min = rn_e3zps_min
304      e3zps_rat = rn_e3zps_rat
305      nmsh      = nn_msh
306      nacc      = nn_acc
307      atfp      = rn_atfp
308      rdt       = rn_rdt
309      rdtmin    = rn_rdtmin
310      rdtmax    = rn_rdtmin
311      rdth      = rn_rdth
312
313#if defined key_netcdf4
314      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
315      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
316      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
317907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
318
319      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
320      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
321908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
322      IF(lwm) WRITE( numond, namnc4 )
323
324      IF(lwp) THEN                        ! control print
325         WRITE(numout,*)
326         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
327         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
328         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
329         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
330         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
331      ENDIF
332
333      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
334      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
335      snc4set%ni   = nn_nchunks_i
336      snc4set%nj   = nn_nchunks_j
337      snc4set%nk   = nn_nchunks_k
338      snc4set%luse = ln_nc4zip
339#else
340      snc4set%luse = .FALSE.        ! No NetCDF 4 case
341#endif
342      !
343   END SUBROUTINE dom_nam
344
345
346   SUBROUTINE dom_ctl
347      !!----------------------------------------------------------------------
348      !!                     ***  ROUTINE dom_ctl  ***
349      !!
350      !! ** Purpose :   Domain control.
351      !!
352      !! ** Method  :   compute and print extrema of masked scale factors
353      !!----------------------------------------------------------------------
354      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
355      INTEGER, DIMENSION(2) ::   iloc   !
356      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
357      !!----------------------------------------------------------------------
358      !
359      IF(lk_mpp) THEN
360         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
361         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
362         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
363         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
364      ELSE
365         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
366         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
367         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
368         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
369
370         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
371         iimi1 = iloc(1) + nimpp - 1
372         ijmi1 = iloc(2) + njmpp - 1
373         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
374         iimi2 = iloc(1) + nimpp - 1
375         ijmi2 = iloc(2) + njmpp - 1
376         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
377         iima1 = iloc(1) + nimpp - 1
378         ijma1 = iloc(2) + njmpp - 1
379         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
380         iima2 = iloc(1) + nimpp - 1
381         ijma2 = iloc(2) + njmpp - 1
382      ENDIF
383      IF(lwp) THEN
384         WRITE(numout,*)
385         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
386         WRITE(numout,*) '~~~~~~~'
387         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
388         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
389         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
390         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
391      ENDIF
392      !
393   END SUBROUTINE dom_ctl
394
395
396   SUBROUTINE dom_stiff
397      !!----------------------------------------------------------------------
398      !!                  ***  ROUTINE dom_stiff  ***
399      !!                     
400      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
401      !!
402      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
403      !!                Save the maximum in the vertical direction
404      !!                (this number is only relevant in s-coordinates)
405      !!
406      !!                Haney, R. L., 1991: On the pressure gradient force
407      !!                over steep topography in sigma coordinate ocean models.
408      !!                J. Phys. Oceanogr., 21, 610???619.
409      !!----------------------------------------------------------------------
410      INTEGER  ::   ji, jj, jk 
411      REAL(wp) ::   zrxmax
412      REAL(wp), DIMENSION(4) :: zr1
413      !!----------------------------------------------------------------------
414      rx1(:,:) = 0._wp
415      zrxmax   = 0._wp
416      zr1(:)   = 0._wp
417      !
418      DO ji = 2, jpim1
419         DO jj = 2, jpjm1
420            DO jk = 1, jpkm1
421               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji-1,jj  ,jk  )  & 
422                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1)) &
423                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji-1,jj  ,jk  )  &
424                    &                         -gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1) + rsmall) )
425               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw_0(ji+1,jj  ,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
426                    &                         +gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
427                    &                        /(gdepw_0(ji+1,jj  ,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
428                    &                         -gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
429               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw_0(ji  ,jj+1,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
430                    &                         +gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
431                    &                        /(gdepw_0(ji  ,jj+1,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
432                    &                         -gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
433               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji  ,jj-1,jk  )  &
434                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1)) &
435                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji  ,jj-1,jk  )  &
436                    &                         -gdepw_0(ji,  jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1) + rsmall) )
437               zrxmax = MAXVAL(zr1(1:4))
438               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
439            END DO
440         END DO
441      END DO
442      CALL lbc_lnk( rx1, 'T', 1. )
443      !
444      zrxmax = MAXVAL( rx1 )
445      !
446      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
447      !
448      IF(lwp) THEN
449         WRITE(numout,*)
450         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
451         WRITE(numout,*) '~~~~~~~~~'
452      ENDIF
453      !
454   END SUBROUTINE dom_stiff
455
456   !!======================================================================
457END MODULE domain
Note: See TracBrowser for help on using the repository browser.