New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domain.F90 in branches/UKMO/dev_r7573_xios_write/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_r7573_xios_write/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 8310

Last change on this file since 8310 was 8310, checked in by andmirek, 6 years ago

#1882 works with AGRIF and few small fixes/changes
File size: 25.4 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   USE iom_def, ONLY:lwxios, wxioso ! write restart flag and output type
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   dom_init   ! called by opa.F90
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49   !!-------------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
53   !!-------------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE dom_init
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE dom_init  ***
59      !!                   
60      !! ** Purpose :   Domain initialization. Call the routines that are
61      !!              required to create the arrays which define the space
62      !!              and time domain of the ocean model.
63      !!
64      !! ** Method  : - dom_msk: compute the masks from the bathymetry file
65      !!              - dom_hgr: compute or read the horizontal grid-point position
66      !!                         and scale factors, and the coriolis factor
67      !!              - dom_zgr: define the vertical coordinate and the bathymetry
68      !!              - dom_stp: defined the model time step
69      !!              - dom_wri: create the meshmask file if nmsh=1
70      !!              - 1D configuration, move Coriolis, u and v at T-point
71      !!----------------------------------------------------------------------
72      INTEGER ::   jk          ! dummy loop argument
73      INTEGER ::   iconf = 0   ! local integers
74      !!----------------------------------------------------------------------
75      !
76      IF( nn_timing == 1 )   CALL timing_start('dom_init')
77      !
78      IF(lwp) THEN
79         WRITE(numout,*)
80         WRITE(numout,*) 'dom_init : domain initialization'
81         WRITE(numout,*) '~~~~~~~~'
82      ENDIF
83      !
84                             CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
85                             CALL dom_clo      ! Closed seas and lake
86                             CALL dom_hgr      ! Horizontal mesh
87                             CALL dom_zgr      ! Vertical mesh and bathymetry
88                             CALL dom_msk      ! Masks
89      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
90      !
91      ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at T-points
92      hu_0(:,:) = 0.0_wp                       ! Reference ocean depth at U-points
93      hv_0(:,:) = 0.0_wp                       ! Reference ocean depth at V-points
94      DO jk = 1, jpk
95         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
96         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
97         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
98      END DO
99      !
100      IF( lk_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_i(:,:) ) * umask_i(:,:)
115      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
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      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
139         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , ln_rstdate, nn_rstctl,   &
140         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
141         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler, &
142         &             nn_wxios
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      nn_wxios = 0
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,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
175         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
176         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
177         WRITE(numout,*) '      restart logical                 ln_rstart  = ' , ln_rstart
178         WRITE(numout,*) '      datestamping of restarts        ln_rstdate  = ', ln_rstdate
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         IF( ln_rst_list ) THEN
187            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
188         ELSE
189            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
190         ENDIF
191         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
192         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
193         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
194         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
195         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
196         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
197         IF( TRIM(Agrif_CFixed()) == '0' ) &
198     &      WRITE(numout,*) '      Write restart using XIOS        nn_wxios   = ', nn_wxios
199      ENDIF
200
201      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
202      cexper = cn_exp
203      nrstdt = nn_rstctl
204      nit000 = nn_it000
205      nitend = nn_itend
206      ndate0 = nn_date0
207      nleapy = nn_leapy
208      ninist = nn_istate
209      nstock = nn_stock
210      nstocklist = nn_stocklist
211      nwrite = nn_write
212      neuler = nn_euler
213      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
214         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
215         CALL ctl_warn( ctmp1 )
216         neuler = 0
217      ENDIF
218
219      !                             ! control of output frequency
220      IF ( nstock == 0 .OR. nstock > nitend ) THEN
221         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
222         CALL ctl_warn( ctmp1 )
223         nstock = nitend
224      ENDIF
225      IF ( nwrite == 0 ) THEN
226         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
227         CALL ctl_warn( ctmp1 )
228         nwrite = nitend
229      ENDIF
230
231#if defined key_agrif
232      IF( Agrif_Root() ) THEN
233#endif
234      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
235      CASE (  1 ) 
236         CALL ioconf_calendar('gregorian')
237         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
238      CASE (  0 )
239         CALL ioconf_calendar('noleap')
240         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
241      CASE ( 30 )
242         CALL ioconf_calendar('360d')
243         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
244      END SELECT
245#if defined key_agrif
246      ENDIF
247#endif
248
249      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
250      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
251903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
252 
253      !
254      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
255      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
256904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
257      IF(lwm) WRITE ( numond, namdom )
258
259      IF(lwp) THEN
260         WRITE(numout,*)
261         WRITE(numout,*) '   Namelist namdom : space & time domain'
262         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
263         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
264         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
265         WRITE(numout,*) '      min number of ocean level (<0)       '
266         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
267         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
268         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
269         WRITE(numout,*) '           = 0   no file created           '
270         WRITE(numout,*) '           = 1   mesh_mask                 '
271         WRITE(numout,*) '           = 2   mesh and mask             '
272         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
273         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
274         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
275         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
276         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
277         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
278         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
279         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
280         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
281         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
282         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
283         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
284         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
285         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
286         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
287         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
288         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
289         WRITE(numout,*) '                                        ppa0            = ', ppa0
290         WRITE(numout,*) '                                        ppa1            = ', ppa1
291         WRITE(numout,*) '                                        ppkth           = ', ppkth
292         WRITE(numout,*) '                                        ppacr           = ', ppacr
293         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
294         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
295         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
296         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
297         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
298         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
299      ENDIF
300
301      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
302      e3zps_min = rn_e3zps_min
303      e3zps_rat = rn_e3zps_rat
304      nmsh      = nn_msh
305      nacc      = nn_acc
306      atfp      = rn_atfp
307      rdt       = rn_rdt
308      rdtmin    = rn_rdtmin
309      rdtmax    = rn_rdtmin
310      rdth      = rn_rdth
311
312      IF( TRIM(Agrif_CFixed()) == '0' ) THEN
313!set output file type for XIOS based on NEMO namelist
314         if (nn_wxios > 0) lwxios = .TRUE. 
315         wxioso = nn_wxios
316      ELSE
317         IF(lwp) THEN
318            write(numout,*)
319            write(numout,*) "AGRIF: nn_wxios will be ingored. See setting for NEMO" 
320            write(numout,*)
321         ENDIF
322      ENDIF
323
324      REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection
325      READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)
326905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )
327
328      REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection
329      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
330906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
331      IF(lwm) WRITE( numond, namcla )
332
333      IF(lwp) THEN
334         WRITE(numout,*)
335         WRITE(numout,*) '   Namelist namcla'
336         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
337      ENDIF
338      IF ( nn_cla .EQ. 1 ) THEN
339         IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2
340            CONTINUE
341         ELSE
342            CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' )
343         ENDIF
344      ENDIF
345
346#if defined key_netcdf4
347      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
348      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
349      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
350907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
351
352      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
353      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
354908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
355      IF(lwm) WRITE( numond, namnc4 )
356
357      IF(lwp) THEN                        ! control print
358         WRITE(numout,*)
359         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
360         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
361         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
362         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
363         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
364      ENDIF
365
366      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
367      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
368      snc4set%ni   = nn_nchunks_i
369      snc4set%nj   = nn_nchunks_j
370      snc4set%nk   = nn_nchunks_k
371      snc4set%luse = ln_nc4zip
372#else
373      snc4set%luse = .FALSE.        ! No NetCDF 4 case
374#endif
375      !
376   END SUBROUTINE dom_nam
377
378
379   SUBROUTINE dom_ctl
380      !!----------------------------------------------------------------------
381      !!                     ***  ROUTINE dom_ctl  ***
382      !!
383      !! ** Purpose :   Domain control.
384      !!
385      !! ** Method  :   compute and print extrema of masked scale factors
386      !!----------------------------------------------------------------------
387      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
388      INTEGER, DIMENSION(2) ::   iloc   !
389      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
390      !!----------------------------------------------------------------------
391      !
392      IF(lk_mpp) THEN
393         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
394         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
395         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
396         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
397      ELSE
398         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
399         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
400         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
401         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
402
403         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
404         iimi1 = iloc(1) + nimpp - 1
405         ijmi1 = iloc(2) + njmpp - 1
406         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
407         iimi2 = iloc(1) + nimpp - 1
408         ijmi2 = iloc(2) + njmpp - 1
409         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
410         iima1 = iloc(1) + nimpp - 1
411         ijma1 = iloc(2) + njmpp - 1
412         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
413         iima2 = iloc(1) + nimpp - 1
414         ijma2 = iloc(2) + njmpp - 1
415      ENDIF
416      IF(lwp) THEN
417         WRITE(numout,*)
418         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
419         WRITE(numout,*) '~~~~~~~'
420         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
421         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
422         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
423         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
424      ENDIF
425      !
426   END SUBROUTINE dom_ctl
427
428   SUBROUTINE dom_stiff
429      !!----------------------------------------------------------------------
430      !!                  ***  ROUTINE dom_stiff  ***
431      !!                     
432      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
433      !!
434      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
435      !!                Save the maximum in the vertical direction
436      !!                (this number is only relevant in s-coordinates)
437      !!
438      !!                Haney, R. L., 1991: On the pressure gradient force
439      !!                over steep topography in sigma coordinate ocean models.
440      !!                J. Phys. Oceanogr., 21, 610???619.
441      !!----------------------------------------------------------------------
442      INTEGER  ::   ji, jj, jk 
443      REAL(wp) ::   zrxmax
444      REAL(wp), DIMENSION(4) :: zr1
445      !!----------------------------------------------------------------------
446      rx1(:,:) = 0.e0
447      zrxmax   = 0.e0
448      zr1(:)   = 0.e0
449     
450      DO ji = 2, jpim1
451         DO jj = 2, jpjm1
452            DO jk = 1, jpkm1
453               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji-1,jj  ,jk  )  & 
454                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1)) &
455                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji-1,jj  ,jk  )  &
456                    &                         -gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1) + rsmall) )
457               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw_0(ji+1,jj  ,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
458                    &                         +gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
459                    &                        /(gdepw_0(ji+1,jj  ,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
460                    &                         -gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
461               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw_0(ji  ,jj+1,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
462                    &                         +gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
463                    &                        /(gdepw_0(ji  ,jj+1,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
464                    &                         -gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
465               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji  ,jj-1,jk  )  &
466                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1)) &
467                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji  ,jj-1,jk  )  &
468                    &                         -gdepw_0(ji,  jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1) + rsmall) )
469               zrxmax = MAXVAL(zr1(1:4))
470               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
471            END DO
472         END DO
473      END DO
474
475      CALL lbc_lnk( rx1, 'T', 1. )
476
477      zrxmax = MAXVAL(rx1)
478
479      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
480
481      IF(lwp) THEN
482         WRITE(numout,*)
483         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
484         WRITE(numout,*) '~~~~~~~~~'
485      ENDIF
486
487   END SUBROUTINE dom_stiff
488
489
490
491   !!======================================================================
492END MODULE domain
Note: See TracBrowser for help on using the repository browser.