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

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 5779

Last change on this file since 5779 was 5779, checked in by mathiot, 9 years ago

ISF coupling branch: correct some compilation issues, remove code related to MISOMIP/ISOMIP+ and polishing
  • Property svn:keywords set to Id
File size: 24.9 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   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   dom_init   ! called by opa.F90
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48   !!-------------------------------------------------------------------------
49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
50   !! $Id$
51   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
52   !!-------------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE dom_init
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE dom_init  ***
58      !!                   
59      !! ** Purpose :   Domain initialization. Call the routines that are
60      !!              required to create the arrays which define the space
61      !!              and time domain of the ocean model.
62      !!
63      !! ** Method  : - dom_msk: compute the masks from the bathymetry file
64      !!              - dom_hgr: compute or read the horizontal grid-point position
65      !!                         and scale factors, and the coriolis factor
66      !!              - dom_zgr: define the vertical coordinate and the bathymetry
67      !!              - dom_stp: defined the model time step
68      !!              - dom_wri: create the meshmask file if nmsh=1
69      !!              - 1D configuration, move Coriolis, u and v at T-point
70      !!----------------------------------------------------------------------
71      INTEGER ::   jk          ! dummy loop argument
72      INTEGER ::   iconf = 0   ! local integers
73      !!----------------------------------------------------------------------
74      !
75      IF( nn_timing == 1 )   CALL timing_start('dom_init')
76      !
77      IF(lwp) THEN
78         WRITE(numout,*)
79         WRITE(numout,*) 'dom_init : domain initialization'
80         WRITE(numout,*) '~~~~~~~~'
81      ENDIF
82      !
83                             CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
84                             CALL dom_clo      ! Closed seas and lake
85                             CALL dom_hgr      ! Horizontal mesh
86                             CALL dom_zgr      ! Vertical mesh and bathymetry
87                             CALL dom_msk      ! Masks
88      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
89      !
90      ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at T-points
91      hu_0(:,:) = 0.0_wp                       ! Reference ocean depth at U-points
92      hv_0(:,:) = 0.0_wp                       ! Reference ocean depth at V-points
93      DO jk = 1, jpk
94         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
95         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
96         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
97      END DO
98      !
99      IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh
100      !
101      IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
102      !
103      !
104      hu(:,:) = 0._wp                          ! Ocean depth at U-points
105      hv(:,:) = 0._wp                          ! Ocean depth at V-points
106      ht(:,:) = 0._wp                          ! Ocean depth at T-points
107      DO jk = 1, jpkm1
108         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
109         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
110         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
111      END DO
112      !                                        ! Inverse of the local depth
113      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
114      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
115
116                             CALL dom_stp      ! time step
117      !
118      IF( nmsh /= 0 .AND. .NOT. ln_iscpl )                         CALL dom_wri      ! Create a domain file
119      IF( nmsh /= 0 .AND.       ln_iscpl .AND. .NOT. ln_rstart )   CALL dom_wri      ! Create a domain file
120      !
121      IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control
122      !
123      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
124      !
125   END SUBROUTINE dom_init
126
127
128   SUBROUTINE dom_nam
129      !!----------------------------------------------------------------------
130      !!                     ***  ROUTINE dom_nam  ***
131      !!                   
132      !! ** Purpose :   read domaine namelists and print the variables.
133      !!
134      !! ** input   : - namrun namelist
135      !!              - namdom namelist
136      !!              - namcla namelist
137      !!              - namnc4 namelist   ! "key_netcdf4" only
138      !!----------------------------------------------------------------------
139      USE ioipsl
140      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
141         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
142         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
143         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler, ln_iscpl
144      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   &
145         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  &
146         &             rn_rdtmax, rn_rdth     , nn_closea , ln_crs,    &
147         &             jphgr_msh, &
148         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
149         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
150         &             ppa2, ppkth2, ppacr2
151      NAMELIST/namcla/ nn_cla
152#if defined key_netcdf4
153      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
154#endif
155      INTEGER  ::   ios                 ! Local integer output status for namelist read
156      !!----------------------------------------------------------------------
157
158      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
159      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
160901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
161
162      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
163      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
164902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
165      IF(lwm) WRITE ( numond, namrun )
166      !
167      IF(lwp) THEN                  ! control print
168         WRITE(numout,*)
169         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
170         WRITE(numout,*) '~~~~~~~ '
171         WRITE(numout,*) '   Namelist namrun'
172         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
173         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
174         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
175         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
176         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
177         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
178         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
179         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
180         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
181         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
182         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
183         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
184         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
185         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
186         IF( ln_rst_list ) THEN
187            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
188         ELSE
189            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
190         ENDIF
191         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
192         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
193         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
194         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
195         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
196         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
197         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl   = ', ln_iscpl
198      ENDIF
199
200      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
201      cexper = cn_exp
202      nrstdt = nn_rstctl
203      nit000 = nn_it000
204      nitend = nn_itend
205      ndate0 = nn_date0
206      nleapy = nn_leapy
207      ninist = nn_istate
208      nstock = nn_stock
209      nstocklist = nn_stocklist
210      nwrite = nn_write
211      neuler = nn_euler
212      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
213         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
214         CALL ctl_warn( ctmp1 )
215         neuler = 0
216      ENDIF
217
218      !                             ! control of output frequency
219      IF ( nstock == 0 .OR. nstock > nitend ) THEN
220         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
221         CALL ctl_warn( ctmp1 )
222         nstock = nitend
223      ENDIF
224      IF ( nwrite == 0 ) THEN
225         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
226         CALL ctl_warn( ctmp1 )
227         nwrite = nitend
228      ENDIF
229
230#if defined key_agrif
231      IF( Agrif_Root() ) THEN
232#endif
233      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
234      CASE (  1 ) 
235         CALL ioconf_calendar('gregorian')
236         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
237      CASE (  0 )
238         CALL ioconf_calendar('noleap')
239         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
240      CASE ( 30 )
241         CALL ioconf_calendar('360d')
242         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
243      END SELECT
244#if defined key_agrif
245      ENDIF
246#endif
247
248      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
249      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
250903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
251 
252      !
253      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
254      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
255904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
256      IF(lwm) WRITE ( numond, namdom )
257
258      IF(lwp) THEN
259         WRITE(numout,*)
260         WRITE(numout,*) '   Namelist namdom : space & time domain'
261         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
262         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
263         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
264         WRITE(numout,*) '      min number of ocean level (<0)       '
265         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
266         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
267         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
268         WRITE(numout,*) '           = 0   no file created           '
269         WRITE(numout,*) '           = 1   mesh_mask                 '
270         WRITE(numout,*) '           = 2   mesh and mask             '
271         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
272         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
273         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
274         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
275         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
276         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
277         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
278         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
279         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
280         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
281         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
282         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
283         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
284         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
285         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
286         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
287         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
288         WRITE(numout,*) '                                        ppa0            = ', ppa0
289         WRITE(numout,*) '                                        ppa1            = ', ppa1
290         WRITE(numout,*) '                                        ppkth           = ', ppkth
291         WRITE(numout,*) '                                        ppacr           = ', ppacr
292         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
293         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
294         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
295         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
296         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
297         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
298      ENDIF
299
300      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
301      e3zps_min = rn_e3zps_min
302      e3zps_rat = rn_e3zps_rat
303      nmsh      = nn_msh
304      nacc      = nn_acc
305      atfp      = rn_atfp
306      rdt       = rn_rdt
307      rdtmin    = rn_rdtmin
308      rdtmax    = rn_rdtmin
309      rdth      = rn_rdth
310
311      REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection
312      READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)
313905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )
314
315      REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection
316      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
317906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
318      IF(lwm) WRITE( numond, namcla )
319
320      IF(lwp) THEN
321         WRITE(numout,*)
322         WRITE(numout,*) '   Namelist namcla'
323         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
324      ENDIF
325      IF ( nn_cla .EQ. 1 ) THEN
326         IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2
327            CONTINUE
328         ELSE
329            CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' )
330         ENDIF
331      ENDIF
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   SUBROUTINE dom_stiff
416      !!----------------------------------------------------------------------
417      !!                  ***  ROUTINE dom_stiff  ***
418      !!                     
419      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
420      !!
421      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
422      !!                Save the maximum in the vertical direction
423      !!                (this number is only relevant in s-coordinates)
424      !!
425      !!                Haney, R. L., 1991: On the pressure gradient force
426      !!                over steep topography in sigma coordinate ocean models.
427      !!                J. Phys. Oceanogr., 21, 610???619.
428      !!----------------------------------------------------------------------
429      INTEGER  ::   ji, jj, jk 
430      REAL(wp) ::   zrxmax
431      REAL(wp), DIMENSION(4) :: zr1
432      !!----------------------------------------------------------------------
433      rx1(:,:) = 0.e0
434      zrxmax   = 0.e0
435      zr1(:)   = 0.e0
436     
437      DO ji = 2, jpim1
438         DO jj = 2, jpjm1
439            DO jk = 1, jpkm1
440               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji-1,jj  ,jk  )  & 
441                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1)) &
442                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji-1,jj  ,jk  )  &
443                    &                         -gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji-1,jj  ,jk+1) + rsmall) )
444               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw_0(ji+1,jj  ,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
445                    &                         +gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
446                    &                        /(gdepw_0(ji+1,jj  ,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
447                    &                         -gdepw_0(ji+1,jj  ,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
448               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw_0(ji  ,jj+1,jk  )-gdepw_0(ji  ,jj  ,jk  )  &
449                    &                         +gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1)) &
450                    &                        /(gdepw_0(ji  ,jj+1,jk  )+gdepw_0(ji  ,jj  ,jk  )  &
451                    &                         -gdepw_0(ji  ,jj+1,jk+1)-gdepw_0(ji  ,jj  ,jk+1) + rsmall) )
452               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw_0(ji  ,jj  ,jk  )-gdepw_0(ji  ,jj-1,jk  )  &
453                    &                         +gdepw_0(ji  ,jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1)) &
454                    &                        /(gdepw_0(ji  ,jj  ,jk  )+gdepw_0(ji  ,jj-1,jk  )  &
455                    &                         -gdepw_0(ji,  jj  ,jk+1)-gdepw_0(ji  ,jj-1,jk+1) + rsmall) )
456               zrxmax = MAXVAL(zr1(1:4))
457               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
458            END DO
459         END DO
460      END DO
461
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
477
478   !!======================================================================
479END MODULE domain
Note: See TracBrowser for help on using the repository browser.