source: branches/UKMO/dev_r5518_GO6_package_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 10005

Last change on this file since 10005 was 10005, checked in by timgraham, 2 years ago

Included all of the functional changes from STARTHOUR branch in branches/2014

File size: 24.9 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_i(:,:) ) * umask_i(:,:)
114      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
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 , ln_rstdate, nn_rstctl,   &
139         &             nn_it000, nn_itend  , nn_date0    , nn_time0,      nn_leapy     , nn_istate , nn_stock ,   &
140         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , 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,*) '      datestamping of restarts        ln_rstdate  = ', ln_rstdate
177         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
178         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
179         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
180         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
181         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
182         WRITE(numout,*) '      initial time of day in hhmm     nn_time0   = ', nn_time0
183         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
184         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
185         IF( ln_rst_list ) THEN
186            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
187         ELSE
188            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
189         ENDIF
190         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
191         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
192         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
193         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
194         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
195         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
196      ENDIF
197
198      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
199      cexper = cn_exp
200      nrstdt = nn_rstctl
201      nit000 = nn_it000
202      nitend = nn_itend
203      ndate0 = nn_date0
204      nleapy = nn_leapy
205      ninist = nn_istate
206      nstock = nn_stock
207      nstocklist = nn_stocklist
208      nwrite = nn_write
209      neuler = nn_euler
210      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
211         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
212         CALL ctl_warn( ctmp1 )
213         neuler = 0
214      ENDIF
215
216      !                             ! control of output frequency
217      IF ( nstock == 0 .OR. nstock > nitend ) THEN
218         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
219         CALL ctl_warn( ctmp1 )
220         nstock = nitend
221      ENDIF
222      IF ( nwrite == 0 ) THEN
223         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
224         CALL ctl_warn( ctmp1 )
225         nwrite = nitend
226      ENDIF
227
228#if defined key_agrif
229      IF( Agrif_Root() ) THEN
230#endif
231      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
232      CASE (  1 ) 
233         CALL ioconf_calendar('gregorian')
234         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
235      CASE (  0 )
236         CALL ioconf_calendar('noleap')
237         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
238      CASE ( 30 )
239         CALL ioconf_calendar('360d')
240         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
241      END SELECT
242#if defined key_agrif
243      ENDIF
244#endif
245
246      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
247      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
248903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
249 
250      !
251      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
252      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
253904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
254      IF(lwm) WRITE ( numond, namdom )
255
256      IF(lwp) THEN
257         WRITE(numout,*)
258         WRITE(numout,*) '   Namelist namdom : space & time domain'
259         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
260         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
261         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
262         WRITE(numout,*) '      min number of ocean level (<0)       '
263         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
264         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
265         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
266         WRITE(numout,*) '           = 0   no file created           '
267         WRITE(numout,*) '           = 1   mesh_mask                 '
268         WRITE(numout,*) '           = 2   mesh and mask             '
269         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
270         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
271         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
272         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
273         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
274         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
275         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
276         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
277         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
278         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
279         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
280         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
281         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
282         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
283         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
284         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
285         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
286         WRITE(numout,*) '                                        ppa0            = ', ppa0
287         WRITE(numout,*) '                                        ppa1            = ', ppa1
288         WRITE(numout,*) '                                        ppkth           = ', ppkth
289         WRITE(numout,*) '                                        ppacr           = ', ppacr
290         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
291         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
292         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
293         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
294         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
295         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
296      ENDIF
297
298      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
299      e3zps_min = rn_e3zps_min
300      e3zps_rat = rn_e3zps_rat
301      nmsh      = nn_msh
302      nacc      = nn_acc
303      atfp      = rn_atfp
304      rdt       = rn_rdt
305      rdtmin    = rn_rdtmin
306      rdtmax    = rn_rdtmin
307      rdth      = rn_rdth
308
309      REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection
310      READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)
311905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )
312
313      REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection
314      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
315906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
316      IF(lwm) WRITE( numond, namcla )
317
318      IF(lwp) THEN
319         WRITE(numout,*)
320         WRITE(numout,*) '   Namelist namcla'
321         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
322      ENDIF
323      IF ( nn_cla .EQ. 1 ) THEN
324         IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2
325            CONTINUE
326         ELSE
327            CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' )
328         ENDIF
329      ENDIF
330
331#if defined key_netcdf4
332      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
333      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
334      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
335907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
336
337      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
338      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
339908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
340      IF(lwm) WRITE( numond, namnc4 )
341
342      IF(lwp) THEN                        ! control print
343         WRITE(numout,*)
344         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
345         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
346         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
347         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
348         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
349      ENDIF
350
351      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
352      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
353      snc4set%ni   = nn_nchunks_i
354      snc4set%nj   = nn_nchunks_j
355      snc4set%nk   = nn_nchunks_k
356      snc4set%luse = ln_nc4zip
357#else
358      snc4set%luse = .FALSE.        ! No NetCDF 4 case
359#endif
360      !
361   END SUBROUTINE dom_nam
362
363
364   SUBROUTINE dom_ctl
365      !!----------------------------------------------------------------------
366      !!                     ***  ROUTINE dom_ctl  ***
367      !!
368      !! ** Purpose :   Domain control.
369      !!
370      !! ** Method  :   compute and print extrema of masked scale factors
371      !!----------------------------------------------------------------------
372      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
373      INTEGER, DIMENSION(2) ::   iloc   !
374      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
375      !!----------------------------------------------------------------------
376      !
377      IF(lk_mpp) THEN
378         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
379         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
380         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
381         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
382      ELSE
383         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
384         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
385         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
386         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
387
388         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
389         iimi1 = iloc(1) + nimpp - 1
390         ijmi1 = iloc(2) + njmpp - 1
391         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
392         iimi2 = iloc(1) + nimpp - 1
393         ijmi2 = iloc(2) + njmpp - 1
394         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
395         iima1 = iloc(1) + nimpp - 1
396         ijma1 = iloc(2) + njmpp - 1
397         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
398         iima2 = iloc(1) + nimpp - 1
399         ijma2 = iloc(2) + njmpp - 1
400      ENDIF
401      IF(lwp) THEN
402         WRITE(numout,*)
403         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
404         WRITE(numout,*) '~~~~~~~'
405         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
406         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
407         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
408         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
409      ENDIF
410      !
411   END SUBROUTINE dom_ctl
412
413   SUBROUTINE dom_stiff
414      !!----------------------------------------------------------------------
415      !!                  ***  ROUTINE dom_stiff  ***
416      !!                     
417      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
418      !!
419      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
420      !!                Save the maximum in the vertical direction
421      !!                (this number is only relevant in s-coordinates)
422      !!
423      !!                Haney, R. L., 1991: On the pressure gradient force
424      !!                over steep topography in sigma coordinate ocean models.
425      !!                J. Phys. Oceanogr., 21, 610???619.
426      !!----------------------------------------------------------------------
427      INTEGER  ::   ji, jj, jk 
428      REAL(wp) ::   zrxmax
429      REAL(wp), DIMENSION(4) :: zr1
430      !!----------------------------------------------------------------------
431      rx1(:,:) = 0.e0
432      zrxmax   = 0.e0
433      zr1(:)   = 0.e0
434     
435      DO ji = 2, jpim1
436         DO jj = 2, jpjm1
437            DO jk = 1, jpkm1
438               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji-1,jj  ,jk  )  & 
439                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1)) &
440                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji-1,jj  ,jk  )  &
441                    &                         -gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1) + rsmall) )
442               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw_0(ji+1,jj  ,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
443                    &                         +gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
444                    &                        /(gdepw_0(ji+1,jj  ,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
445                    &                         -gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
446               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw_0(ji  ,jj+1,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
447                    &                         +gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
448                    &                        /(gdepw_0(ji  ,jj+1,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
449                    &                         -gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
450               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji  ,jj-1,jk  )  &
451                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1)) &
452                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji  ,jj-1,jk  )  &
453                    &                         -gdepw_0(ji,  jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1) + rsmall) )
454               zrxmax = MAXVAL(zr1(1:4))
455               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
456            END DO
457         END DO
458      END DO
459
460      CALL lbc_lnk( rx1, 'T', 1. )
461
462      zrxmax = MAXVAL(rx1)
463
464      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
465
466      IF(lwp) THEN
467         WRITE(numout,*)
468         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
469         WRITE(numout,*) '~~~~~~~~~'
470      ENDIF
471
472   END SUBROUTINE dom_stiff
473
474
475
476   !!======================================================================
477END MODULE domain
Note: See TracBrowser for help on using the repository browser.