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

source: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 6069

Last change on this file since 6069 was 6069, checked in by timgraham, 8 years ago

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