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

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

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 5845

Last change on this file since 5845 was 5845, checked in by gm, 9 years ago

#1613: vvl by default: suppression of domzgr_substitute.h90

  • Property svn:keywords set to Id
File size: 25.2 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   !!             -   ! 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 argument
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) * tmask(:,:,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( lk_vvl ) THEN                         ! time varying : initialize before/now/after variables
106         CALL dom_vvl_init 
107         !
108      ELSE                                      ! Fix in time : set to the reference one for all
109         !    before         !          now          !       after         !
110         gdept_b = gdept_0   ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
111         gdepw_b = gdepw_0   ;   gdepw_n = gdepw_0   !        ---          !
112                                 gde3w_n = gde3w_0   !        ---          !
113         !                                                                 
114           e3t_b =   e3t_0   ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
115           e3u_b =   e3u_0   ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
116           e3v_b =   e3v_0   ;     e3v_n =   e3u_0   ;   e3v_a =  e3v_0    !
117                             ;     e3f_n =   e3f_0   !        ---          !
118           e3w_b =   e3w_0   ;     e3w_n =   e3w_0   !        ---          !
119          e3uw_b =  e3uw_0   ;    e3uw_n =  e3uw_0   !        ---          !
120          e3vw_b =  e3vw_0   ;    e3vw_n =  e3vw_0   !        ---          !
121         !
122         !                                            !
123         CALL wrk_alloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
124         !
125         z1_hu_0(:,:) = 1._wp / ( hu_0(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)    ! _i mask due to ISF
126         z1_hv_0(:,:) = 1._wp / ( hv_0(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
127         !
128         !        before       !         now         !       after         !
129         ;                     ;     ht_n =    hu_0  ;    ht_a =    hu_0   ! water column thickness
130         ;     hu_b =    hu_0  ;     hu_n =    hu_0  ;    hu_a =    hu_0   !
131         ;     hv_b =    hv_0  ;     hv_n =    hv_0  ;    hv_a =    hv_0   !
132         ;  r1_hu_b = z1_hu_0  ;  r1_hu_n = z1_hu_0  ; r1_hu_a = z1_hu_0   ! inverse of water column thickness
133         ;  r1_hv_b = z1_hv_0  ;  r1_hv_n = z1_hv_0  ; r1_hv_a = z1_hv_0   !
134         !
135         CALL wrk_dealloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
136      ENDIF
137      !
138      IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
139      !
140                             CALL dom_stp       ! time step
141      IF( nmsh /= 0      )   CALL dom_wri       ! Create a domain file
142      IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control
143      !
144      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
145      !
146   END SUBROUTINE dom_init
147
148
149   SUBROUTINE dom_nam
150      !!----------------------------------------------------------------------
151      !!                     ***  ROUTINE dom_nam  ***
152      !!                   
153      !! ** Purpose :   read domaine namelists and print the variables.
154      !!
155      !! ** input   : - namrun namelist
156      !!              - namdom namelist
157      !!              - namnc4 namelist   ! "key_netcdf4" only
158      !!----------------------------------------------------------------------
159      USE ioipsl
160      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
161         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
162         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
163         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler
164      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   &
165         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  &
166         &             rn_rdtmax, rn_rdth     , nn_closea , ln_crs,    &
167         &             jphgr_msh, &
168         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
169         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
170         &             ppa2, ppkth2, ppacr2
171#if defined key_netcdf4
172      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
173#endif
174      INTEGER  ::   ios                 ! Local integer output status for namelist read
175      !!----------------------------------------------------------------------
176
177      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
178      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
179901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
180
181      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
182      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
183902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
184      IF(lwm) WRITE ( numond, namrun )
185      !
186      IF(lwp) THEN                  ! control print
187         WRITE(numout,*)
188         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
189         WRITE(numout,*) '~~~~~~~ '
190         WRITE(numout,*) '   Namelist namrun'
191         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
192         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
193         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
194         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
195         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
196         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
197         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
198         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
199         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
200         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
201         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
202         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
203         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
204         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
205         IF( ln_rst_list ) THEN
206            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
207         ELSE
208            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
209         ENDIF
210         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
211         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
212         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
213         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
214         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
215         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
216      ENDIF
217
218      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
219      cexper = cn_exp
220      nrstdt = nn_rstctl
221      nit000 = nn_it000
222      nitend = nn_itend
223      ndate0 = nn_date0
224      nleapy = nn_leapy
225      ninist = nn_istate
226      nstock = nn_stock
227      nstocklist = nn_stocklist
228      nwrite = nn_write
229      neuler = nn_euler
230      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
231         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
232         CALL ctl_warn( ctmp1 )
233         neuler = 0
234      ENDIF
235
236      !                             ! control of output frequency
237      IF ( nstock == 0 .OR. nstock > nitend ) THEN
238         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
239         CALL ctl_warn( ctmp1 )
240         nstock = nitend
241      ENDIF
242      IF ( nwrite == 0 ) THEN
243         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
244         CALL ctl_warn( ctmp1 )
245         nwrite = nitend
246      ENDIF
247
248#if defined key_agrif
249      IF( Agrif_Root() ) THEN
250#endif
251      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
252      CASE (  1 ) 
253         CALL ioconf_calendar('gregorian')
254         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
255      CASE (  0 )
256         CALL ioconf_calendar('noleap')
257         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
258      CASE ( 30 )
259         CALL ioconf_calendar('360d')
260         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
261      END SELECT
262#if defined key_agrif
263      ENDIF
264#endif
265
266      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
267      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
268903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
269 
270      !
271      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
272      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
273904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
274      IF(lwm) WRITE ( numond, namdom )
275      !
276      IF(lwp) THEN
277         WRITE(numout,*)
278         WRITE(numout,*) '   Namelist namdom : space & time domain'
279         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
280         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
281         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
282         WRITE(numout,*) '      min number of ocean level (<0)       '
283         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
284         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
285         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
286         WRITE(numout,*) '           = 0   no file created           '
287         WRITE(numout,*) '           = 1   mesh_mask                 '
288         WRITE(numout,*) '           = 2   mesh and mask             '
289         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
290         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
291         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
292         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
293         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
294         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
295         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
296         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
297         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
298         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
299         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
300         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
301         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
302         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
303         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
304         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
305         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
306         WRITE(numout,*) '                                        ppa0            = ', ppa0
307         WRITE(numout,*) '                                        ppa1            = ', ppa1
308         WRITE(numout,*) '                                        ppkth           = ', ppkth
309         WRITE(numout,*) '                                        ppacr           = ', ppacr
310         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
311         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
312         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
313         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
314         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
315         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
316      ENDIF
317      !
318      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
319      e3zps_min = rn_e3zps_min
320      e3zps_rat = rn_e3zps_rat
321      nmsh      = nn_msh
322      nacc      = nn_acc
323      atfp      = rn_atfp
324      rdt       = rn_rdt
325      rdtmin    = rn_rdtmin
326      rdtmax    = rn_rdtmin
327      rdth      = rn_rdth
328
329#if defined key_netcdf4
330      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
331      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
332      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
333907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
334
335      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
336      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
337908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
338      IF(lwm) WRITE( numond, namnc4 )
339
340      IF(lwp) THEN                        ! control print
341         WRITE(numout,*)
342         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
343         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
344         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
345         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
346         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
347      ENDIF
348
349      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
350      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
351      snc4set%ni   = nn_nchunks_i
352      snc4set%nj   = nn_nchunks_j
353      snc4set%nk   = nn_nchunks_k
354      snc4set%luse = ln_nc4zip
355#else
356      snc4set%luse = .FALSE.        ! No NetCDF 4 case
357#endif
358      !
359   END SUBROUTINE dom_nam
360
361
362   SUBROUTINE dom_ctl
363      !!----------------------------------------------------------------------
364      !!                     ***  ROUTINE dom_ctl  ***
365      !!
366      !! ** Purpose :   Domain control.
367      !!
368      !! ** Method  :   compute and print extrema of masked scale factors
369      !!----------------------------------------------------------------------
370      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
371      INTEGER, DIMENSION(2) ::   iloc   !
372      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
373      !!----------------------------------------------------------------------
374      !
375      IF(lk_mpp) THEN
376         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
377         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
378         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
379         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
380      ELSE
381         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
382         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
383         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
384         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
385
386         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
387         iimi1 = iloc(1) + nimpp - 1
388         ijmi1 = iloc(2) + njmpp - 1
389         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
390         iimi2 = iloc(1) + nimpp - 1
391         ijmi2 = iloc(2) + njmpp - 1
392         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
393         iima1 = iloc(1) + nimpp - 1
394         ijma1 = iloc(2) + njmpp - 1
395         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
396         iima2 = iloc(1) + nimpp - 1
397         ijma2 = iloc(2) + njmpp - 1
398      ENDIF
399      IF(lwp) THEN
400         WRITE(numout,*)
401         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
402         WRITE(numout,*) '~~~~~~~'
403         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
404         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
405         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
406         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
407      ENDIF
408      !
409   END SUBROUTINE dom_ctl
410
411
412   SUBROUTINE dom_stiff
413      !!----------------------------------------------------------------------
414      !!                  ***  ROUTINE dom_stiff  ***
415      !!                     
416      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
417      !!
418      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
419      !!                Save the maximum in the vertical direction
420      !!                (this number is only relevant in s-coordinates)
421      !!
422      !!                Haney, R. L., 1991: On the pressure gradient force
423      !!                over steep topography in sigma coordinate ocean models.
424      !!                J. Phys. Oceanogr., 21, 610???619.
425      !!----------------------------------------------------------------------
426      INTEGER  ::   ji, jj, jk 
427      REAL(wp) ::   zrxmax
428      REAL(wp), DIMENSION(4) :: zr1
429      !!----------------------------------------------------------------------
430      rx1(:,:) = 0._wp
431      zrxmax   = 0._wp
432      zr1(:)   = 0._wp
433      !
434      DO ji = 2, jpim1
435         DO jj = 2, jpjm1
436            DO jk = 1, jpkm1
437               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
438                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
439                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
440                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
441               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
442                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
443                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
444                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
445               zr1(3) =ABS(  (  gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
446                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
447                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
448                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
449               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
450                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
451                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
452                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
453               zrxmax = MAXVAL( zr1(1:4) )
454               rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax )
455            END DO
456         END DO
457      END DO
458      CALL lbc_lnk( rx1, 'T', 1. )
459      !
460      zrxmax = MAXVAL( rx1 )
461      !
462      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
463      !
464      IF(lwp) THEN
465         WRITE(numout,*)
466         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
467         WRITE(numout,*) '~~~~~~~~~'
468      ENDIF
469      !
470   END SUBROUTINE dom_stiff
471
472   !!======================================================================
473END MODULE domain
Note: See TracBrowser for help on using the repository browser.