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

source: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

Last change on this file was 5750, checked in by rfurner, 9 years ago

changes to add restart directory names

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