source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

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