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 @ 7512

Last change on this file since 7512 was 7508, checked in by mocavero, 7 years ago

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