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

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 6748

Last change on this file since 6748 was 6748, checked in by mocavero, 8 years ago

GYRE hybrid parallelization
  • Property svn:keywords set to Id
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   !!            3.7  !  2015-11  (G. Madec, A. Coward)  time varying zgr by default
16   !!----------------------------------------------------------------------
17   
18   !!----------------------------------------------------------------------
19   !!   dom_init       : initialize the space and time domain
20   !!   dom_nam        : read and contral domain namelists
21   !!   dom_ctl        : control print for the ocean domain
22   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean variables
25   USE dom_oce         ! domain: ocean
26   USE sbc_oce         ! surface boundary condition: ocean
27   USE phycst          ! physical constants
28   USE closea          ! closed seas
29   USE domhgr          ! domain: set the horizontal mesh
30   USE domzgr          ! domain: set the vertical mesh
31   USE domstp          ! domain: set the time-step
32   USE dommsk          ! domain: set the mask system
33   USE domwri          ! domain: write the meshmask file
34   USE domvvl          ! variable volume
35   USE c1d             ! 1D vertical configuration
36   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
37   !
38   USE in_out_manager  ! I/O manager
39   USE wrk_nemo        ! Memory Allocation
40   USE lib_mpp         ! distributed memory computing library
41   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
42   USE timing          ! Timing
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   dom_init   ! called by opa.F90
48
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, jj, ji          ! dummy loop indices
73      INTEGER ::   iconf = 0   ! local integers
74      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_hu_0, z1_hv_0
75      !!----------------------------------------------------------------------
76      !
77      IF( nn_timing == 1 )   CALL timing_start('dom_init')
78      !
79      IF(lwp) THEN
80         WRITE(numout,*)
81         WRITE(numout,*) 'dom_init : domain initialization'
82         WRITE(numout,*) '~~~~~~~~'
83      ENDIF
84      !
85      !                       !==  Reference coordinate system  ==!
86      !
87                     CALL dom_nam               ! read namelist ( namrun, namdom )
88                     CALL dom_clo               ! Closed seas and lake
89                     CALL dom_hgr               ! Horizontal mesh
90                     CALL dom_zgr               ! Vertical mesh and bathymetry
91                     CALL dom_msk               ! Masks
92      IF( ln_sco )   CALL dom_stiff             ! Maximum stiffness ratio/hydrostatic consistency
93      !
94!$OMP PARALLEL WORKSHARE
95      ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1)   ! Reference ocean thickness
96      hu_0(:,:) = e3u_0(:,:,1) * umask(:,:,1)
97      hv_0(:,:) = e3v_0(:,:,1) * vmask(:,:,1)
98!$OMP END PARALLEL WORKSHARE
99      DO jk = 2, jpk
100!$OMP PARALLEL DO schedule(static) private(jj,ji)
101         DO jj =1, jpj
102            DO ji=1, jpi
103                ht_0(ji,jj) = ht_0(ji,jj) + e3t_0(ji,jj,jk) * tmask(ji,jj,jk)
104                hu_0(ji,jj) = hu_0(ji,jj) + e3u_0(ji,jj,jk) * umask(ji,jj,jk)
105                hv_0(ji,jj) = hv_0(ji,jj) + e3v_0(ji,jj,jk) * vmask(ji,jj,jk)
106            END DO
107         END DO
108      END DO
109      !
110      !              !==  time varying part of coordinate system  ==!
111      !
112      IF( ln_linssh ) THEN          ! Fix in time : set to the reference one for all
113         !       before        !          now          !       after         !
114         ;  gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
115         ;  gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
116         ;                     ;   gde3w_n = gde3w_0   !        ---          !
117         !                                                                 
118         ;    e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
119         ;    e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
120         ;    e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
121         ;                     ;     e3f_n =   e3f_0   !        ---          !
122         ;    e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
123         ;   e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
124         ;   e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
125         !
126         CALL wrk_alloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
127         !
128!$OMP PARALLEL WORKSHARE
129         z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
130         z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
131!$OMP END PARALLEL WORKSHARE
132         !
133         !        before       !          now          !       after         !
134         ;                     ;      ht_n =    ht_0   !                     ! water column thickness
135         ;     hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   !
136         ;     hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   !
137         ;  r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness
138         ;  r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   !
139         !
140         CALL wrk_dealloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
141         !
142      ELSE                         ! time varying : initialize before/now/after variables
143         !
144         CALL dom_vvl_init 
145         !
146      ENDIF
147      !
148      IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
149      !
150                             CALL dom_stp       ! time step
151      IF( nmsh /= 0 .AND. .NOT. ln_iscpl )                         CALL dom_wri      ! Create a domain file
152      IF( nmsh /= 0 .AND.       ln_iscpl .AND. .NOT. ln_rstart )   CALL dom_wri      ! Create a domain file
153      IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control
154      !
155      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
156      !
157   END SUBROUTINE dom_init
158
159
160   SUBROUTINE dom_nam
161      !!----------------------------------------------------------------------
162      !!                     ***  ROUTINE dom_nam  ***
163      !!                   
164      !! ** Purpose :   read domaine namelists and print the variables.
165      !!
166      !! ** input   : - namrun namelist
167      !!              - namdom namelist
168      !!              - namnc4 namelist   ! "key_netcdf4" only
169      !!----------------------------------------------------------------------
170      USE ioipsl
171      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 &
172                       nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     &
173         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
174         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     &
175         &             ln_cfmeta, ln_iscpl
176      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, &
177         &             rn_atfp , rn_rdt   , nn_closea   , ln_crs      , jphgr_msh ,                  &
178         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         &
179         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  &
180         &             ppa2, ppkth2, ppacr2
181#if defined key_netcdf4
182      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
183#endif
184      INTEGER  ::   ios                 ! Local integer output status for namelist read
185      !!----------------------------------------------------------------------
186
187      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
188      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
189901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
190
191      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
192      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
193902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
194      IF(lwm) WRITE ( numond, namrun )
195      !
196      IF(lwp) THEN                  ! control print
197         WRITE(numout,*)
198         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
199         WRITE(numout,*) '~~~~~~~ '
200         WRITE(numout,*) '   Namelist namrun'
201         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
202         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
203         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
204         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
205         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
206         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
207         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
208         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
209         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
210         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
211         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
212         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
213         WRITE(numout,*) '      initial time of day in hhmm     nn_time0   = ', nn_time0
214         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
215         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
216         IF( ln_rst_list ) THEN
217            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
218         ELSE
219            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
220         ENDIF
221         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
222         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
223         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
224         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
225         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
226         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl   = ', ln_iscpl
227      ENDIF
228
229      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
230      cexper = cn_exp
231      nrstdt = nn_rstctl
232      nit000 = nn_it000
233      nitend = nn_itend
234      ndate0 = nn_date0
235      nleapy = nn_leapy
236      ninist = nn_istate
237      nstock = nn_stock
238      nstocklist = nn_stocklist
239      nwrite = nn_write
240      neuler = nn_euler
241      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
242         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
243         CALL ctl_warn( ctmp1 )
244         neuler = 0
245      ENDIF
246
247      !                             ! control of output frequency
248      IF ( nstock == 0 .OR. nstock > nitend ) THEN
249         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
250         CALL ctl_warn( ctmp1 )
251         nstock = nitend
252      ENDIF
253      IF ( nwrite == 0 ) THEN
254         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
255         CALL ctl_warn( ctmp1 )
256         nwrite = nitend
257      ENDIF
258
259#if defined key_agrif
260      IF( Agrif_Root() ) THEN
261#endif
262      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
263      CASE (  1 ) 
264         CALL ioconf_calendar('gregorian')
265         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
266      CASE (  0 )
267         CALL ioconf_calendar('noleap')
268         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
269      CASE ( 30 )
270         CALL ioconf_calendar('360d')
271         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
272      END SELECT
273#if defined key_agrif
274      ENDIF
275#endif
276
277      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
278      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
279903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
280 
281      !
282      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
283      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
284904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
285      IF(lwm) WRITE ( numond, namdom )
286      !
287      IF(lwp) THEN
288         WRITE(numout,*)
289         WRITE(numout,*) '   Namelist namdom : space & time domain'
290         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
291         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
292         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
293         WRITE(numout,*) '      min number of ocean level (<0)       '
294         WRITE(numout,*) '      treshold to open the isf cavity   rn_isfhmin   = ', rn_isfhmin, ' (m)'
295         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
296         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
297         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
298         WRITE(numout,*) '           = 0   no file created           '
299         WRITE(numout,*) '           = 1   mesh_mask                 '
300         WRITE(numout,*) '           = 2   mesh and mask             '
301         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
302         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
303         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
304         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
305         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
306         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
307         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
308         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
309         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
310         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
311         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
312         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
313         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
314         WRITE(numout,*) '                                        ppa0            = ', ppa0
315         WRITE(numout,*) '                                        ppa1            = ', ppa1
316         WRITE(numout,*) '                                        ppkth           = ', ppkth
317         WRITE(numout,*) '                                        ppacr           = ', ppacr
318         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
319         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
320         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
321         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
322         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
323         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
324      ENDIF
325      !
326      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
327      e3zps_min = rn_e3zps_min
328      e3zps_rat = rn_e3zps_rat
329      nmsh      = nn_msh
330      atfp      = rn_atfp
331      rdt       = rn_rdt
332
333#if defined key_netcdf4
334      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
335      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
336      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
337907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
338
339      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
340      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
341908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
342      IF(lwm) WRITE( numond, namnc4 )
343
344      IF(lwp) THEN                        ! control print
345         WRITE(numout,*)
346         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
347         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
348         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
349         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
350         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
351      ENDIF
352
353      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
354      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
355      snc4set%ni   = nn_nchunks_i
356      snc4set%nj   = nn_nchunks_j
357      snc4set%nk   = nn_nchunks_k
358      snc4set%luse = ln_nc4zip
359#else
360      snc4set%luse = .FALSE.        ! No NetCDF 4 case
361#endif
362      !
363   END SUBROUTINE dom_nam
364
365
366   SUBROUTINE dom_ctl
367      !!----------------------------------------------------------------------
368      !!                     ***  ROUTINE dom_ctl  ***
369      !!
370      !! ** Purpose :   Domain control.
371      !!
372      !! ** Method  :   compute and print extrema of masked scale factors
373      !!----------------------------------------------------------------------
374      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
375      INTEGER, DIMENSION(2) ::   iloc   !
376      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
377      !!----------------------------------------------------------------------
378      !
379      IF(lk_mpp) THEN
380         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
381         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
382         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
383         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
384      ELSE
385         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
386         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
387         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
388         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
389
390         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
391         iimi1 = iloc(1) + nimpp - 1
392         ijmi1 = iloc(2) + njmpp - 1
393         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
394         iimi2 = iloc(1) + nimpp - 1
395         ijmi2 = iloc(2) + njmpp - 1
396         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
397         iima1 = iloc(1) + nimpp - 1
398         ijma1 = iloc(2) + njmpp - 1
399         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
400         iima2 = iloc(1) + nimpp - 1
401         ijma2 = iloc(2) + njmpp - 1
402      ENDIF
403      IF(lwp) THEN
404         WRITE(numout,*)
405         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
406         WRITE(numout,*) '~~~~~~~'
407         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
408         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
409         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
410         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
411      ENDIF
412      !
413   END SUBROUTINE dom_ctl
414
415
416   SUBROUTINE dom_stiff
417      !!----------------------------------------------------------------------
418      !!                  ***  ROUTINE dom_stiff  ***
419      !!                     
420      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
421      !!
422      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
423      !!                Save the maximum in the vertical direction
424      !!                (this number is only relevant in s-coordinates)
425      !!
426      !!                Haney, R. L., 1991: On the pressure gradient force
427      !!                over steep topography in sigma coordinate ocean models.
428      !!                J. Phys. Oceanogr., 21, 610???619.
429      !!----------------------------------------------------------------------
430      INTEGER  ::   ji, jj, jk 
431      REAL(wp) ::   zrxmax
432      REAL(wp), DIMENSION(4) ::   zr1
433      !!----------------------------------------------------------------------
434      rx1(:,:) = 0._wp
435      zrxmax   = 0._wp
436      zr1(:)   = 0._wp
437      !
438      DO ji = 2, jpim1
439         DO jj = 2, jpjm1
440            DO jk = 1, jpkm1
441               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
442                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
443                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
444                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
445               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
446                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
447                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
448                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
449               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
450                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
451                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
452                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
453               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
454                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
455                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
456                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
457               zrxmax = MAXVAL( zr1(1:4) )
458               rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax )
459            END DO
460         END DO
461      END DO
462      CALL lbc_lnk( rx1, 'T', 1. )
463      !
464      zrxmax = MAXVAL( rx1 )
465      !
466      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
467      !
468      IF(lwp) THEN
469         WRITE(numout,*)
470         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
471         WRITE(numout,*) '~~~~~~~~~'
472      ENDIF
473      !
474   END SUBROUTINE dom_stiff
475
476   !!======================================================================
477END MODULE domain
Note: See TracBrowser for help on using the repository browser.