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

source: branches/UKMO/test_moci_test_suite/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 8269

Last change on this file since 8269 was 8269, checked in by andmirek, 7 years ago

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