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

Last change on this file since 5997 was 5997, checked in by timgraham, 8 years ago

Merged branches/2014/dev_r4650_UKMO7_STARTHOUR into branch

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