New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domain.F90 in branches/UKMO/2015_CO6_CO5_zenv_wr_direct_dwl_temp/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

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

Last change on this file since 5389 was 5389, checked in by deazer, 9 years ago

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