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

source: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 4780

Last change on this file since 4780 was 4780, checked in by edblockley, 10 years ago

Second commit in UKMO11 development branch.

This change allows the user to replace the nn_stock frequency-based restart dump writing functionality with a list-based version (nn_stocklist).
This is conterolled using the logical ln_rst_list which defaults to false.
At present the list is hard-wired to have maximum 10 entries but this could be modified if required.

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