source: branches/UKMO/2015_CO6_CO5_zenv_wr_direct_dwl_temp/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 5418

Last change on this file since 5418 was 5418, checked in by deazer, 6 years ago

Removed SVN KEYWORDS ready for adding code changes before fcm merges

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