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/dev_r10171_test_crs_AMM7/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_r10171_test_crs_AMM7/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 10207

Last change on this file since 10207 was 10207, checked in by cmao, 6 years ago

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