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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 25.1 KB
RevLine 
[3]1MODULE domain
2   !!==============================================================================
3   !!                       ***  MODULE domain   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
[1438]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
[2528]13   !!            3.3  !  2010-11  (G. Madec)  initialisation in C1D configuration
[4152]14   !!            3.6  !  2013     ( J. Simeon, C. Calone, G. Madec, C. Ethe ) Online coarsening of outputs
[3]15   !!----------------------------------------------------------------------
[1438]16   
17   !!----------------------------------------------------------------------
[3]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   !!----------------------------------------------------------------------
[2528]22   USE oce             ! ocean variables
23   USE dom_oce         ! domain: ocean
[888]24   USE sbc_oce         ! surface boundary condition: ocean
[719]25   USE phycst          ! physical constants
[1601]26   USE closea          ! closed seas
[719]27   USE in_out_manager  ! I/O manager
[3]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
[592]35   USE domvvl          ! variable volume
[2528]36   USE c1d             ! 1D vertical configuration
37   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
[3294]38   USE timing          ! Timing
[3680]39   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
[3]40
41   IMPLICIT NONE
42   PRIVATE
43
[1438]44   PUBLIC   dom_init   ! called by opa.F90
[3]45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
[1438]48   !!-------------------------------------------------------------------------
[2528]49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]50   !! $Id$
[2528]51   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
[1438]52   !!-------------------------------------------------------------------------
[3]53CONTAINS
54
55   SUBROUTINE dom_init
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE dom_init  ***
58      !!                   
59      !! ** Purpose :   Domain initialization. Call the routines that are
[1601]60      !!              required to create the arrays which define the space
61      !!              and time domain of the ocean model.
[3]62      !!
[1601]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
[2528]69      !!              - 1D configuration, move Coriolis, u and v at T-point
[3]70      !!----------------------------------------------------------------------
[3764]71      INTEGER ::   jk          ! dummy loop argument
72      INTEGER ::   iconf = 0   ! local integers
[3]73      !!----------------------------------------------------------------------
[1601]74      !
[3764]75      IF( nn_timing == 1 )   CALL timing_start('dom_init')
[3294]76      !
[3]77      IF(lwp) THEN
78         WRITE(numout,*)
79         WRITE(numout,*) 'dom_init : domain initialization'
80         WRITE(numout,*) '~~~~~~~~'
[11101]81         IF(lflush) CALL flush(numout)
[3]82      ENDIF
[1601]83      !
84                             CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
85                             CALL dom_clo      ! Closed seas and lake
86                             CALL dom_hgr      ! Horizontal mesh
87                             CALL dom_zgr      ! Vertical mesh and bathymetry
88                             CALL dom_msk      ! Masks
[3680]89      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
[4490]90      !
91      ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at T-points
92      hu_0(:,:) = 0.0_wp                       ! Reference ocean depth at U-points
93      hv_0(:,:) = 0.0_wp                       ! Reference ocean depth at V-points
94      DO jk = 1, jpk
95         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
96         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
97         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
98      END DO
99      !
[4292]100      IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh
[1601]101      !
[3764]102      IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
[2528]103      !
[4370]104      !
105      hu(:,:) = 0._wp                          ! Ocean depth at U-points
106      hv(:,:) = 0._wp                          ! Ocean depth at V-points
107      ht(:,:) = 0._wp                          ! Ocean depth at T-points
108      DO jk = 1, jpkm1
[4292]109         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
110         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
[4370]111         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
[3]112      END DO
[1601]113      !                                        ! Inverse of the local depth
[4990]114      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
115      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
[216]116
[1601]117                             CALL dom_stp      ! time step
118      IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file
119      IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control
[1438]120      !
[3764]121      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
[3294]122      !
[3]123   END SUBROUTINE dom_init
124
125
126   SUBROUTINE dom_nam
127      !!----------------------------------------------------------------------
128      !!                     ***  ROUTINE dom_nam  ***
129      !!                   
130      !! ** Purpose :   read domaine namelists and print the variables.
131      !!
132      !! ** input   : - namrun namelist
133      !!              - namdom namelist
134      !!              - namcla namelist
[2528]135      !!              - namnc4 namelist   ! "key_netcdf4" only
[3]136      !!----------------------------------------------------------------------
137      USE ioipsl
[5341]138      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
[6491]139         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , ln_rstdate, nn_rstctl,   &
[1601]140         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
[5363]141         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler
[4245]142      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   &
143         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  &
[4370]144         &             rn_rdtmax, rn_rdth     , nn_closea , ln_crs,    &
[4147]145         &             jphgr_msh, &
146         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
147         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
148         &             ppa2, ppkth2, ppacr2
[1601]149      NAMELIST/namcla/ nn_cla
[2528]150#if defined key_netcdf4
151      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
152#endif
[4147]153      INTEGER  ::   ios                 ! Local integer output status for namelist read
[3]154      !!----------------------------------------------------------------------
155
[4147]156      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
157      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
158901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
159
160      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
161      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
162902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
[11101]163      IF(lwm .AND. nprint > 2) WRITE ( numond, namrun )
[1601]164      !
165      IF(lwp) THEN                  ! control print
[3]166         WRITE(numout,*)
167         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
168         WRITE(numout,*) '~~~~~~~ '
[1601]169         WRITE(numout,*) '   Namelist namrun'
170         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
171         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
[4147]172         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
[5341]173         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
[4147]174         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
[5341]175         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
[6491]176         WRITE(numout,*) '      restart logical                 ln_rstart  = ' , ln_rstart
177         WRITE(numout,*) '      datestamping of restarts        ln_rstdate  = ', ln_rstdate
[4370]178         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
[1604]179         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
[1601]180         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
181         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
182         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
183         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
184         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
[5341]185         IF( ln_rst_list ) THEN
186            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
187         ELSE
188            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
189         ENDIF
[1601]190         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
191         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
192         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
[5363]193         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
[1601]194         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
195         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
[11101]196         IF(lflush) CALL flush(numout)
[3]197      ENDIF
198
[1601]199      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
200      cexper = cn_exp
201      nrstdt = nn_rstctl
202      nit000 = nn_it000
203      nitend = nn_itend
204      ndate0 = nn_date0
205      nleapy = nn_leapy
206      ninist = nn_istate
207      nstock = nn_stock
[5341]208      nstocklist = nn_stocklist
[1601]209      nwrite = nn_write
[4370]210      neuler = nn_euler
[5341]211      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
[4370]212         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
213         CALL ctl_warn( ctmp1 )
214         neuler = 0
215      ENDIF
[3]216
[1601]217      !                             ! control of output frequency
[1335]218      IF ( nstock == 0 .OR. nstock > nitend ) THEN
[1601]219         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
[783]220         CALL ctl_warn( ctmp1 )
[1335]221         nstock = nitend
[3]222      ENDIF
223      IF ( nwrite == 0 ) THEN
[1601]224         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
[783]225         CALL ctl_warn( ctmp1 )
226         nwrite = nitend
[3]227      ENDIF
228
[2528]229#if defined key_agrif
[1601]230      IF( Agrif_Root() ) THEN
[2528]231#endif
232      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
233      CASE (  1 ) 
234         CALL ioconf_calendar('gregorian')
235         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
236      CASE (  0 )
237         CALL ioconf_calendar('noleap')
238         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
239      CASE ( 30 )
240         CALL ioconf_calendar('360d')
241         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
242      END SELECT
243#if defined key_agrif
[1601]244      ENDIF
[2528]245#endif
[3]246
[4147]247      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
248      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
249903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
[4152]250 
251      !
[4147]252      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
253      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
254904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
[11101]255      IF(lwm .AND. nprint > 2) WRITE ( numond, namdom )
[4147]256
[3]257      IF(lwp) THEN
[72]258         WRITE(numout,*)
[1601]259         WRITE(numout,*) '   Namelist namdom : space & time domain'
260         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
[4245]261         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
[2528]262         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
263         WRITE(numout,*) '      min number of ocean level (<0)       '
[1601]264         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
265         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
266         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
[2528]267         WRITE(numout,*) '           = 0   no file created           '
268         WRITE(numout,*) '           = 1   mesh_mask                 '
269         WRITE(numout,*) '           = 2   mesh and mask             '
270         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
[4152]271         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
272         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
273         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
274         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
275         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
276         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
277         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
278         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
[4147]279         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
280         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
281         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
282         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
283         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
284         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
285         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
286         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
287         WRITE(numout,*) '                                        ppa0            = ', ppa0
288         WRITE(numout,*) '                                        ppa1            = ', ppa1
289         WRITE(numout,*) '                                        ppkth           = ', ppkth
290         WRITE(numout,*) '                                        ppacr           = ', ppacr
291         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
292         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
293         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
294         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
295         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
296         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
[11101]297         IF(lflush) CALL flush(numout)
[223]298      ENDIF
299
[1601]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
[4147]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 )
[11101]318      IF(lwm .AND. nprint > 2) WRITE( numond, namcla )
[4147]319
[3]320      IF(lwp) THEN
[72]321         WRITE(numout,*)
[1601]322         WRITE(numout,*) '   Namelist namcla'
323         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
[11101]324         IF(lflush) CALL flush(numout)
[3]325      ENDIF
[4147]326      IF ( nn_cla .EQ. 1 ) THEN
327         IF  ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2
328            CONTINUE
329         ELSE
330            CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' )
331         ENDIF
332      ENDIF
[3]333
[2528]334#if defined key_netcdf4
335      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
[4147]336      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
337      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
338907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
339
340      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
341      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
342908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
[11101]343      IF(lwm .AND. nprint > 2) WRITE( numond, namnc4 )
[4147]344
[2528]345      IF(lwp) THEN                        ! control print
346         WRITE(numout,*)
347         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
348         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
349         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
350         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
351         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
[11101]352         IF(lflush) CALL flush(numout)
[2528]353      ENDIF
[1601]354
[2528]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
[1438]364      !
[3]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
[1601]377      INTEGER, DIMENSION(2) ::   iloc   !
[3]378      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
379      !!----------------------------------------------------------------------
[1601]380      !
381      IF(lk_mpp) THEN
[4990]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 )
[181]386      ELSE
[4990]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 )   
[32]391
[4990]392         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
[181]393         iimi1 = iloc(1) + nimpp - 1
394         ijmi1 = iloc(2) + njmpp - 1
[4990]395         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
[181]396         iimi2 = iloc(1) + nimpp - 1
397         ijmi2 = iloc(2) + njmpp - 1
[4990]398         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
[181]399         iima1 = iloc(1) + nimpp - 1
400         ijma1 = iloc(2) + njmpp - 1
[4990]401         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
[181]402         iima2 = iloc(1) + nimpp - 1
403         ijma2 = iloc(2) + njmpp - 1
[32]404      ENDIF
[3]405      IF(lwp) THEN
[1601]406         WRITE(numout,*)
407         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
408         WRITE(numout,*) '~~~~~~~'
[181]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
[11101]413         IF(lflush) CALL flush(numout)
[3]414      ENDIF
[1438]415      !
[3]416   END SUBROUTINE dom_ctl
417
[3680]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.e0
437      zrxmax   = 0.e0
438      zr1(:)   = 0.e0
439     
440      DO ji = 2, jpim1
441         DO jj = 2, jpjm1
442            DO jk = 1, jpkm1
[4292]443               zr1(1) = umask(ji-1,jj  ,jk) *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) )
447               zr1(2) = umask(ji  ,jj  ,jk) *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) )
451               zr1(3) = vmask(ji  ,jj  ,jk) *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) )
455               zr1(4) = vmask(ji  ,jj-1,jk) *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) )
[3680]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
465      CALL lbc_lnk( rx1, 'T', 1. )
466
467      zrxmax = MAXVAL(rx1)
468
469      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
470
471      IF(lwp) THEN
472         WRITE(numout,*)
473         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
474         WRITE(numout,*) '~~~~~~~~~'
[11101]475         IF(lflush) CALL flush(numout)
[3680]476      ENDIF
477
478   END SUBROUTINE dom_stiff
479
480
481
[3]482   !!======================================================================
483END MODULE domain
Note: See TracBrowser for help on using the repository browser.