source: branches/UKMO/dev_r6393_CO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 7019

Last change on this file since 7019 was 7019, checked in by deazer, 5 years ago

Cleared svn keywords

File size: 25.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   !!            3.7  !  2015-11  (G. Madec, A. Coward)  time varying zgr by default
16   !!----------------------------------------------------------------------
17   
18   !!----------------------------------------------------------------------
19   !!   dom_init       : initialize the space and time domain
20   !!   dom_nam        : read and contral domain namelists
21   !!   dom_ctl        : control print for the ocean domain
22   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean variables
25   USE dom_oce         ! domain: ocean
26   USE sbc_oce         ! surface boundary condition: ocean
27   USE phycst          ! physical constants
28   USE closea          ! closed seas
29   USE domhgr          ! domain: set the horizontal mesh
30   USE domzgr          ! domain: set the vertical mesh
31   USE domstp          ! domain: set the time-step
32   USE dommsk          ! domain: set the mask system
33   USE domwri          ! domain: write the meshmask file
34   USE domvvl          ! variable volume
35   USE c1d             ! 1D vertical configuration
36   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
37   !
38   USE in_out_manager  ! I/O manager
39   USE wrk_nemo        ! Memory Allocation
40   USE lib_mpp         ! distributed memory computing library
41   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
42   USE timing          ! Timing
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   dom_init   ! called by opa.F90
48
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 indices
73      INTEGER ::   iconf = 0   ! local integers
74      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_hu_0, z1_hv_0
75      !!----------------------------------------------------------------------
76      !
77      IF( nn_timing == 1 )   CALL timing_start('dom_init')
78      !
79      IF(lwp) THEN
80         WRITE(numout,*)
81         WRITE(numout,*) 'dom_init : domain initialization'
82         WRITE(numout,*) '~~~~~~~~'
83      ENDIF
84      !
85      !                       !==  Reference coordinate system  ==!
86      !
87                     CALL dom_nam               ! read namelist ( namrun, namdom )
88                     CALL dom_clo               ! Closed seas and lake
89                     CALL dom_hgr               ! Horizontal mesh
90                     CALL dom_zgr               ! Vertical mesh and bathymetry
91                     CALL dom_msk               ! Masks
92      IF( ln_sco )   CALL dom_stiff             ! Maximum stiffness ratio/hydrostatic consistency
93      !
94      ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1)   ! Reference ocean thickness
95      hu_0(:,:) = e3u_0(:,:,1) * umask(:,:,1)
96      hv_0(:,:) = e3v_0(:,:,1) * vmask(:,:,1)
97      DO jk = 2, jpk
98         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
99         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
100         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
101      END DO
102      !
103      !              !==  time varying part of coordinate system  ==!
104      !
105      IF( ln_linssh ) THEN          ! Fix in time : set to the reference one for all
106         !       before        !          now          !       after         !
107         ;  gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
108         ;  gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
109         ;                     ;   gde3w_n = gde3w_0   !        ---          !
110         !                                                                 
111         ;    e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
112         ;    e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
113         ;    e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
114         ;                     ;     e3f_n =   e3f_0   !        ---          !
115         ;    e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
116         ;   e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
117         ;   e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
118         !
119         CALL wrk_alloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
120         !
121         z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
122         z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
123         !
124         !        before       !          now          !       after         !
125         ;                     ;      ht_n =    ht_0   !                     ! water column thickness
126         ;     hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   !
127         ;     hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   !
128         ;  r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness
129         ;  r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   !
130         !
131         CALL wrk_dealloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
132         !
133      ELSE                         ! time varying : initialize before/now/after variables
134         !
135         CALL dom_vvl_init 
136         !
137      ENDIF
138      !
139      IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
140      !
141                             CALL dom_stp       ! time step
142      IF( nmsh /= 0 .AND. .NOT. ln_iscpl )                         CALL dom_wri      ! Create a domain file
143      IF( nmsh /= 0 .AND.       ln_iscpl .AND. .NOT. ln_rstart )   CALL dom_wri      ! Create a domain file
144      IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control
145      !
146      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
147      !
148   END SUBROUTINE dom_init
149
150
151   SUBROUTINE dom_nam
152      !!----------------------------------------------------------------------
153      !!                     ***  ROUTINE dom_nam  ***
154      !!                   
155      !! ** Purpose :   read domaine namelists and print the variables.
156      !!
157      !! ** input   : - namrun namelist
158      !!              - namdom namelist
159      !!              - namnc4 namelist   ! "key_netcdf4" only
160      !!----------------------------------------------------------------------
161      USE ioipsl
162      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 &
163                       nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     &
164         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
165         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     &
166         &             ln_cfmeta, ln_iscpl
167      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, &
168         &             rn_atfp , rn_rdt   , nn_closea   , ln_crs      , jphgr_msh ,                  &
169         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         &
170         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  &
171         &             ppa2, ppkth2, ppacr2
172#if defined key_netcdf4
173      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
174#endif
175      INTEGER  ::   ios                 ! Local integer output status for namelist read
176      !!----------------------------------------------------------------------
177
178      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
179      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
180901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
181
182      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
183      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
184902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
185      IF(lwm) WRITE ( numond, namrun )
186      !
187      IF(lwp) THEN                  ! control print
188         WRITE(numout,*)
189         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
190         WRITE(numout,*) '~~~~~~~ '
191         WRITE(numout,*) '   Namelist namrun'
192         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
193         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
194         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
195         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
196         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
197         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
198         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
199         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
200         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
201         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
202         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
203         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
204         WRITE(numout,*) '      initial time of day in hhmm     nn_time0   = ', nn_time0
205         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
206         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
207         IF( ln_rst_list ) THEN
208            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
209         ELSE
210            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
211         ENDIF
212         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
213         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
214         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
215         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
216         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
217         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl   = ', ln_iscpl
218      ENDIF
219
220      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
221      cexper = cn_exp
222      nrstdt = nn_rstctl
223      nit000 = nn_it000
224      nitend = nn_itend
225      ndate0 = nn_date0
226      nleapy = nn_leapy
227      ninist = nn_istate
228      nstock = nn_stock
229      nstocklist = nn_stocklist
230      nwrite = nn_write
231      neuler = nn_euler
232      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
233         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
234         CALL ctl_warn( ctmp1 )
235         neuler = 0
236      ENDIF
237
238      !                             ! control of output frequency
239      IF ( nstock == 0 .OR. nstock > nitend ) THEN
240         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
241         CALL ctl_warn( ctmp1 )
242         nstock = nitend
243      ENDIF
244      IF ( nwrite == 0 ) THEN
245         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
246         CALL ctl_warn( ctmp1 )
247         nwrite = nitend
248      ENDIF
249
250#if defined key_agrif
251      IF( Agrif_Root() ) THEN
252#endif
253      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
254      CASE (  1 ) 
255         CALL ioconf_calendar('gregorian')
256         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
257      CASE (  0 )
258         CALL ioconf_calendar('noleap')
259         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
260      CASE ( 30 )
261         CALL ioconf_calendar('360d')
262         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
263      END SELECT
264#if defined key_agrif
265      ENDIF
266#endif
267
268      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
269      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
270903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
271 
272      !
273      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
274      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
275904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
276      IF(lwm) WRITE ( numond, namdom )
277      !
278      IF(lwp) THEN
279         WRITE(numout,*)
280         WRITE(numout,*) '   Namelist namdom : space & time domain'
281         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
282         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
283         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
284         WRITE(numout,*) '      min number of ocean level (<0)       '
285         WRITE(numout,*) '      treshold to open the isf cavity   rn_isfhmin   = ', rn_isfhmin, ' (m)'
286         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
287         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
288         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
289         WRITE(numout,*) '           = 0   no file created           '
290         WRITE(numout,*) '           = 1   mesh_mask                 '
291         WRITE(numout,*) '           = 2   mesh and mask             '
292         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
293         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
294         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
295         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
296         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
297         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
298         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
299         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
300         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
301         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
302         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
303         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
304         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
305         WRITE(numout,*) '                                        ppa0            = ', ppa0
306         WRITE(numout,*) '                                        ppa1            = ', ppa1
307         WRITE(numout,*) '                                        ppkth           = ', ppkth
308         WRITE(numout,*) '                                        ppacr           = ', ppacr
309         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
310         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
311         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
312         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
313         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
314         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
315      ENDIF
316      !
317      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
318      e3zps_min = rn_e3zps_min
319      e3zps_rat = rn_e3zps_rat
320      nmsh      = nn_msh
321      atfp      = rn_atfp
322      rdt       = rn_rdt
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
407   SUBROUTINE dom_stiff
408      !!----------------------------------------------------------------------
409      !!                  ***  ROUTINE dom_stiff  ***
410      !!                     
411      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
412      !!
413      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
414      !!                Save the maximum in the vertical direction
415      !!                (this number is only relevant in s-coordinates)
416      !!
417      !!                Haney, R. L., 1991: On the pressure gradient force
418      !!                over steep topography in sigma coordinate ocean models.
419      !!                J. Phys. Oceanogr., 21, 610???619.
420      !!----------------------------------------------------------------------
421      INTEGER  ::   ji, jj, jk 
422      REAL(wp) ::   zrxmax
423      REAL(wp), DIMENSION(4) ::   zr1
424      !!----------------------------------------------------------------------
425      rx1(:,:) = 0._wp
426      zrxmax   = 0._wp
427      zr1(:)   = 0._wp
428      !
429      DO ji = 2, jpim1
430         DO jj = 2, jpjm1
431            DO jk = 1, jpkm1
432               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
433                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
434                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
435                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
436               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
437                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
438                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
439                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
440               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
441                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
442                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
443                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
444               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
445                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
446                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
447                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
448               zrxmax = MAXVAL( zr1(1:4) )
449               rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax )
450            END DO
451         END DO
452      END DO
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   !!======================================================================
468END MODULE domain
Note: See TracBrowser for help on using the repository browser.