source: branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 9353

Last change on this file since 9353 was 9353, checked in by jchanut, 2 years ago

Tilde coordinate update: first draft

  • Property svn:keywords set to Id
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   !!----------------------------------------------------------------------
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 trc_oce         ! shared ocean-passive tracers variables
26   USE phycst          ! physical constants
27   USE closea          ! closed seas
28   USE in_out_manager  ! I/O manager
29   USE lib_mpp         ! distributed memory computing library
30
31   USE domhgr          ! domain: set the horizontal mesh
32   USE domzgr          ! domain: set the vertical mesh
33   USE domstp          ! domain: set the time-step
34   USE dommsk          ! domain: set the mask system
35   USE domwri          ! domain: write the meshmask file
36   USE domvvl          ! variable volume
37   USE c1d             ! 1D vertical configuration
38   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
39   USE timing          ! Timing
40   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   dom_init   ! called by opa.F90
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
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 ::   jj, jk      ! dummy loop argument
73      INTEGER ::   iconf = 0   ! local integers
74      !!----------------------------------------------------------------------
75      !
76      IF( nn_timing == 1 )   CALL timing_start('dom_init')
77      !
78      IF(lwp) THEN
79         WRITE(numout,*)
80         WRITE(numout,*) 'dom_init : domain initialization'
81         WRITE(numout,*) '~~~~~~~~'
82      ENDIF
83      !
84                             CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
85                             CALL dom_clo      ! Closed seas and lake
86                             CALL dom_hgr      ! Horizontal mesh
87                             CALL dom_zgr      ! Vertical mesh and bathymetry
88                             CALL dom_msk      ! Masks
89      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
90      !
91      ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at T-points
92      hu_0(:,:) = 0.0_wp                       ! Reference ocean depth at U-points
93      hv_0(:,:) = 0.0_wp                       ! Reference ocean depth at V-points
94      hf_0(:,:) = 0.0_wp                       ! Reference ocean depth at F-points
95      DO jk = 1, jpk
96         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
97         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
98         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
99      END DO
100      !   
101      DO jj=1, jpjm1
102         DO jk = 1, jpk
103            hf_0(:,jj) = hf_0(:,jj) + e3f_0(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk)
104         END DO
105      END DO
106      CALL lbc_lnk( hf_0, 'F', 1._wp )
107      !
108      IF( lk_c1d )           CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
109      !
110      IF( .NOT.lk_offline ) THEN
111        !
112        IF( lk_vvl )         CALL dom_vvl_init ! Vertical variable mesh
113        !
114        hu(:,:) = 0._wp                          ! Ocean depth at U-points
115        hv(:,:) = 0._wp                          ! Ocean depth at V-points
116        ht(:,:) = 0._wp                          ! Ocean depth at T-points
117        DO jk = 1, jpkm1
118           hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
119           hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
120           ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
121        END DO
122        !                                        ! Inverse of the local depth
123        hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
124        hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
125        !
126      ENDIF
127
128                             CALL dom_stp      ! time step
129      IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file
130      IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control
131      !
132      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
133      !
134   END SUBROUTINE dom_init
135
136
137   SUBROUTINE dom_nam
138      !!----------------------------------------------------------------------
139      !!                     ***  ROUTINE dom_nam  ***
140      !!                   
141      !! ** Purpose :   read domaine namelists and print the variables.
142      !!
143      !! ** input   : - namrun namelist
144      !!              - namdom namelist
145      !!              - namcla namelist
146      !!              - namnc4 namelist   ! "key_netcdf4" only
147      !!----------------------------------------------------------------------
148      USE ioipsl
149      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
150         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
151         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
152         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler
153      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   &
154         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  &
155         &             rn_rdtmax, rn_rdth     , nn_closea , ln_crs,    &
156         &             jphgr_msh, &
157         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
158         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
159         &             ppa2, ppkth2, ppacr2
160      NAMELIST/namcla/ nn_cla
161#if defined key_netcdf4
162      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
163#endif
164      INTEGER  ::   ios                 ! Local integer output status for namelist read
165      !!----------------------------------------------------------------------
166
167      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
168      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
169901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
170
171      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
172      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
173902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
174      IF(lwm) WRITE ( numond, namrun )
175      !
176      IF(lwp) THEN                  ! control print
177         WRITE(numout,*)
178         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
179         WRITE(numout,*) '~~~~~~~ '
180         WRITE(numout,*) '   Namelist namrun'
181         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
182         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
183         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
184         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
185         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
186         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
187         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
188         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
189         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
190         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
191         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
192         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
193         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
194         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
195         IF( ln_rst_list ) THEN
196            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
197         ELSE
198            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
199         ENDIF
200         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
201         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
202         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
203         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
204         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
205         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
206      ENDIF
207
208      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
209      cexper = cn_exp
210      nrstdt = nn_rstctl
211      nit000 = nn_it000
212      nitend = nn_itend
213      ndate0 = nn_date0
214      nleapy = nn_leapy
215      ninist = nn_istate
216      nstock = nn_stock
217      nstocklist = nn_stocklist
218      nwrite = nn_write
219      neuler = nn_euler
220      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
221         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
222         CALL ctl_warn( ctmp1 )
223         neuler = 0
224      ENDIF
225
226      !                             ! control of output frequency
227      IF ( nstock == 0 .OR. nstock > nitend ) THEN
228         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
229         CALL ctl_warn( ctmp1 )
230         nstock = nitend
231      ENDIF
232      IF ( nwrite == 0 ) THEN
233         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
234         CALL ctl_warn( ctmp1 )
235         nwrite = nitend
236      ENDIF
237
238#if defined key_agrif
239      IF( Agrif_Root() ) THEN
240#endif
241      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
242      CASE (  1 ) 
243         CALL ioconf_calendar('gregorian')
244         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
245      CASE (  0 )
246         CALL ioconf_calendar('noleap')
247         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
248      CASE ( 30 )
249         CALL ioconf_calendar('360d')
250         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
251      END SELECT
252#if defined key_agrif
253      ENDIF
254#endif
255
256      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
257      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
258903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
259 
260      !
261      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
262      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
263904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
264      IF(lwm) WRITE ( numond, namdom )
265
266      IF(lwp) THEN
267         WRITE(numout,*)
268         WRITE(numout,*) '   Namelist namdom : space & time domain'
269         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
270         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
271         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
272         WRITE(numout,*) '      min number of ocean level (<0)       '
273         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
274         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
275         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
276         WRITE(numout,*) '           = 0   no file created           '
277         WRITE(numout,*) '           = 1   mesh_mask                 '
278         WRITE(numout,*) '           = 2   mesh and mask             '
279         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
280         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
281         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
282         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
283         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
284         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
285         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
286         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
287         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
288         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
289         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
290         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
291         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
292         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
293         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
294         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
295         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
296         WRITE(numout,*) '                                        ppa0            = ', ppa0
297         WRITE(numout,*) '                                        ppa1            = ', ppa1
298         WRITE(numout,*) '                                        ppkth           = ', ppkth
299         WRITE(numout,*) '                                        ppacr           = ', ppacr
300         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
301         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
302         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
303         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
304         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
305         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
306      ENDIF
307
308      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
309      e3zps_min = rn_e3zps_min
310      e3zps_rat = rn_e3zps_rat
311      nmsh      = nn_msh
312      nacc      = nn_acc
313      atfp      = rn_atfp
314      rdt       = rn_rdt
315      rdtmin    = rn_rdtmin
316      rdtmax    = rn_rdtmin
317      rdth      = rn_rdth
318
319      REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection
320      READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)
321905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )
322
323      REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection
324      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
325906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
326      IF(lwm) WRITE( numond, namcla )
327
328      IF(lwp) THEN
329         WRITE(numout,*)
330         WRITE(numout,*) '   Namelist namcla'
331         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
332      ENDIF
333      IF ( nn_cla .EQ. 1 ) THEN
334         IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2
335            CONTINUE
336         ELSE
337            CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' )
338         ENDIF
339      ENDIF
340
341#if defined key_netcdf4
342      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
343      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
344      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
345907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
346
347      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
348      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
349908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
350      IF(lwm) WRITE( numond, namnc4 )
351
352      IF(lwp) THEN                        ! control print
353         WRITE(numout,*)
354         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
355         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
356         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
357         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
358         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
359      ENDIF
360
361      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
362      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
363      snc4set%ni   = nn_nchunks_i
364      snc4set%nj   = nn_nchunks_j
365      snc4set%nk   = nn_nchunks_k
366      snc4set%luse = ln_nc4zip
367#else
368      snc4set%luse = .FALSE.        ! No NetCDF 4 case
369#endif
370      !
371   END SUBROUTINE dom_nam
372
373
374   SUBROUTINE dom_ctl
375      !!----------------------------------------------------------------------
376      !!                     ***  ROUTINE dom_ctl  ***
377      !!
378      !! ** Purpose :   Domain control.
379      !!
380      !! ** Method  :   compute and print extrema of masked scale factors
381      !!----------------------------------------------------------------------
382      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
383      INTEGER, DIMENSION(2) ::   iloc   !
384      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
385      !!----------------------------------------------------------------------
386      !
387      IF(lk_mpp) THEN
388         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
389         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
390         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
391         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
392      ELSE
393         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
394         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
395         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
396         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
397
398         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
399         iimi1 = iloc(1) + nimpp - 1
400         ijmi1 = iloc(2) + njmpp - 1
401         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
402         iimi2 = iloc(1) + nimpp - 1
403         ijmi2 = iloc(2) + njmpp - 1
404         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
405         iima1 = iloc(1) + nimpp - 1
406         ijma1 = iloc(2) + njmpp - 1
407         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
408         iima2 = iloc(1) + nimpp - 1
409         ijma2 = iloc(2) + njmpp - 1
410      ENDIF
411      IF(lwp) THEN
412         WRITE(numout,*)
413         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
414         WRITE(numout,*) '~~~~~~~'
415         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
416         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
417         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
418         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
419      ENDIF
420      !
421   END SUBROUTINE dom_ctl
422
423   SUBROUTINE dom_stiff
424      !!----------------------------------------------------------------------
425      !!                  ***  ROUTINE dom_stiff  ***
426      !!                     
427      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
428      !!
429      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
430      !!                Save the maximum in the vertical direction
431      !!                (this number is only relevant in s-coordinates)
432      !!
433      !!                Haney, R. L., 1991: On the pressure gradient force
434      !!                over steep topography in sigma coordinate ocean models.
435      !!                J. Phys. Oceanogr., 21, 610???619.
436      !!----------------------------------------------------------------------
437      INTEGER  ::   ji, jj, jk 
438      REAL(wp) ::   zrxmax
439      REAL(wp), DIMENSION(4) :: zr1
440      !!----------------------------------------------------------------------
441      rx1(:,:) = 0.e0
442      zrxmax   = 0.e0
443      zr1(:)   = 0.e0
444     
445      DO ji = 2, jpim1
446         DO jj = 2, jpjm1
447            DO jk = 1, jpkm1
448               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji-1,jj  ,jk  )  & 
449                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1)) &
450                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji-1,jj  ,jk  )  &
451                    &                         -gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1) + rsmall) )
452               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw_0(ji+1,jj  ,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
453                    &                         +gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
454                    &                        /(gdepw_0(ji+1,jj  ,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
455                    &                         -gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
456               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw_0(ji  ,jj+1,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
457                    &                         +gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
458                    &                        /(gdepw_0(ji  ,jj+1,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
459                    &                         -gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
460               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji  ,jj-1,jk  )  &
461                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1)) &
462                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji  ,jj-1,jk  )  &
463                    &                         -gdepw_0(ji,  jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1) + rsmall) )
464               zrxmax = MAXVAL(zr1(1:4))
465               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
466            END DO
467         END DO
468      END DO
469
470      CALL lbc_lnk( rx1, 'T', 1. )
471
472      zrxmax = MAXVAL(rx1)
473
474      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
475
476      IF(lwp) THEN
477         WRITE(numout,*)
478         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
479         WRITE(numout,*) '~~~~~~~~~'
480      ENDIF
481
482   END SUBROUTINE dom_stiff
483
484
485
486   !!======================================================================
487END MODULE domain
Note: See TracBrowser for help on using the repository browser.