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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 4294

Last change on this file since 4294 was 4292, checked in by cetlod, 10 years ago

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