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

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

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

Last change on this file was 11738, checked in by marc, 5 years ago

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