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

source: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 4375

Last change on this file since 4375 was 4375, checked in by hliu, 10 years ago

updated gravity filters
  • Property svn:keywords set to Id
File size: 20.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   !!----------------------------------------------------------------------
15   
16   !!----------------------------------------------------------------------
17   !!   dom_init       : initialize the space and time domain
18   !!   dom_nam        : read and contral domain namelists
19   !!   dom_ctl        : control print for the ocean domain
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean variables
22   USE dom_oce         ! domain: ocean
23   USE sbc_oce         ! surface boundary condition: ocean
24   USE phycst          ! physical constants
25   USE closea          ! closed seas
26   USE in_out_manager  ! I/O manager
27   USE lib_mpp         ! distributed memory computing library
28
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   USE timing          ! Timing
38   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   dom_init   ! called by opa.F90
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47   !!-------------------------------------------------------------------------
48   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
49   !! $Id$
50   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
51   !!-------------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE dom_init
55      !!----------------------------------------------------------------------
56      !!                  ***  ROUTINE dom_init  ***
57      !!                   
58      !! ** Purpose :   Domain initialization. Call the routines that are
59      !!              required to create the arrays which define the space
60      !!              and time domain of the ocean model.
61      !!
62      !! ** Method  : - dom_msk: compute the masks from the bathymetry file
63      !!              - dom_hgr: compute or read the horizontal grid-point position
64      !!                         and scale factors, and the coriolis factor
65      !!              - dom_zgr: define the vertical coordinate and the bathymetry
66      !!              - dom_stp: defined the model time step
67      !!              - dom_wri: create the meshmask file if nmsh=1
68      !!              - 1D configuration, move Coriolis, u and v at T-point
69      !!----------------------------------------------------------------------
70      INTEGER ::   jk          ! dummy loop argument
71      INTEGER ::   iconf = 0   ! local integers
72      !!----------------------------------------------------------------------
73      !
74      IF( nn_timing == 1 )   CALL timing_start('dom_init')
75      !
76      IF(lwp) THEN
77         WRITE(numout,*)
78         WRITE(numout,*) 'dom_init : domain initialization'
79         WRITE(numout,*) '~~~~~~~~'
80      ENDIF
81      !
82                             CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
83                             CALL dom_wad      ! read namelist ( namwad ), wetting/drying
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      IF( lk_vvl         )   CALL dom_vvl      ! Vertical variable mesh
90      !
91      IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
92      !
93      hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points
94      hv(:,:) = 0._wp
95      DO jk = 1, jpk
96         hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
97         hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk)
98      END DO
99      !                                        ! Inverse of the local depth
100      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
101      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
102
103                             CALL dom_stp      ! time step
104      IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file
105      IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control
106      !
107      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
108      !
109   END SUBROUTINE dom_init
110
111
112   SUBROUTINE dom_nam
113      !!----------------------------------------------------------------------
114      !!                     ***  ROUTINE dom_nam  ***
115      !!                   
116      !! ** Purpose :   read domaine namelists and print the variables.
117      !!
118      !! ** input   : - namrun namelist
119      !!              - namdom namelist
120      !!              - namcla namelist
121      !!              - namnc4 namelist   ! "key_netcdf4" only
122      !!----------------------------------------------------------------------
123      USE ioipsl
124      NAMELIST/namrun/ nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
125         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
126         &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz
127      NAMELIST/namdom/ nn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin,   &
128         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,            &
129         &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea
130      NAMELIST/namcla/ nn_cla
131#if defined key_netcdf4
132      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
133#endif
134      !!----------------------------------------------------------------------
135
136      REWIND( numnam )              ! Namelist namrun : parameters of the run
137      READ  ( numnam, namrun )
138      !
139      IF(lwp) THEN                  ! control print
140         WRITE(numout,*)
141         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
142         WRITE(numout,*) '~~~~~~~ '
143         WRITE(numout,*) '   Namelist namrun'
144         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
145         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
146         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
147         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
148         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
149         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
150         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
151         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
152         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
153         WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
154         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
155         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
156         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
157         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
158         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
159      ENDIF
160
161      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
162      cexper = cn_exp
163      nrstdt = nn_rstctl
164      nit000 = nn_it000
165      nitend = nn_itend
166      ndate0 = nn_date0
167      nleapy = nn_leapy
168      ninist = nn_istate
169      nstock = nn_stock
170      nwrite = nn_write
171
172
173      !                             ! control of output frequency
174      IF ( nstock == 0 .OR. nstock > nitend ) THEN
175         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
176         CALL ctl_warn( ctmp1 )
177         nstock = nitend
178      ENDIF
179      IF ( nwrite == 0 ) THEN
180         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
181         CALL ctl_warn( ctmp1 )
182         nwrite = nitend
183      ENDIF
184
185#if defined key_agrif
186      IF( Agrif_Root() ) THEN
187#endif
188      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
189      CASE (  1 ) 
190         CALL ioconf_calendar('gregorian')
191         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
192      CASE (  0 )
193         CALL ioconf_calendar('noleap')
194         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
195      CASE ( 30 )
196         CALL ioconf_calendar('360d')
197         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
198      END SELECT
199#if defined key_agrif
200      ENDIF
201#endif
202
203      REWIND( numnam )              ! Namelist namdom : space & time domain (bathymetry, mesh, timestep)
204      READ  ( numnam, namdom )
205
206      IF(lwp) THEN
207         WRITE(numout,*)
208         WRITE(numout,*) '   Namelist namdom : space & time domain'
209         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
210         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
211         WRITE(numout,*) '      min number of ocean level (<0)       '
212         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
213         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
214         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
215         WRITE(numout,*) '           = 0   no file created           '
216         WRITE(numout,*) '           = 1   mesh_mask                 '
217         WRITE(numout,*) '           = 2   mesh and mask             '
218         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
219         WRITE(numout,*) '      ocean time step                      rn_rdt    = ', rn_rdt
220         WRITE(numout,*) '      asselin time filter parameter        rn_atfp   = ', rn_atfp
221         WRITE(numout,*) '      time-splitting: nb of sub time-step  nn_baro   = ', nn_baro
222         WRITE(numout,*) '      acceleration of converge             nn_acc    = ', nn_acc
223         WRITE(numout,*) '        nn_acc=1: surface tracer rdt       rn_rdtmin = ', rn_rdtmin
224         WRITE(numout,*) '                  bottom  tracer rdt       rdtmax    = ', rn_rdtmax
225         WRITE(numout,*) '                  depth of transition      rn_rdth   = ', rn_rdth
226         WRITE(numout,*) '      suppression of closed seas (=0)      nn_closea = ', nn_closea
227      ENDIF
228
229      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
230      e3zps_min = rn_e3zps_min
231      e3zps_rat = rn_e3zps_rat
232      nmsh      = nn_msh
233      nacc      = nn_acc
234      atfp      = rn_atfp
235      rdt       = rn_rdt
236      rdtmin    = rn_rdtmin
237      rdtmax    = rn_rdtmin
238      rdth      = rn_rdth
239
240      REWIND( numnam )              ! Namelist cross land advection
241      READ  ( numnam, namcla )
242      IF(lwp) THEN
243         WRITE(numout,*)
244         WRITE(numout,*) '   Namelist namcla'
245         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
246      ENDIF
247
248#if defined key_netcdf4
249      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
250      REWIND( numnam )                    ! Namelist namnc4 : netcdf4 chunking parameters
251      READ  ( numnam, namnc4 )
252      IF(lwp) THEN                        ! control print
253         WRITE(numout,*)
254         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
255         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
256         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
257         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
258         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
259      ENDIF
260
261      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
262      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
263      snc4set%ni   = nn_nchunks_i
264      snc4set%nj   = nn_nchunks_j
265      snc4set%nk   = nn_nchunks_k
266      snc4set%luse = ln_nc4zip
267#else
268      snc4set%luse = .FALSE.        ! No NetCDF 4 case
269#endif
270      !
271   END SUBROUTINE dom_nam
272
273   SUBROUTINE dom_wad
274      !!----------------------------------------------------------------------
275      !!                     ***  ROUTINE dom_wad  ***
276      !!                   
277      !! ** Purpose :   read wetting/drying namelists and print the variables.
278      !!
279      !! ** input   : - namwad namelist
280      !!----------------------------------------------------------------------
281      USE ioipsl
282      NAMELIST/namwad/ ln_wad , rn_wadmin, rn_wadmine, rn_landele, nn_waditr
283      !!----------------------------------------------------------------------
284
285      REWIND( numnam )              ! Namelist namrun : parameters of the run
286      READ  ( numnam, namwad )
287      !
288      IF(lwp) THEN                  ! control print
289         WRITE(numout,*)
290         WRITE(numout,*) 'dom_wad  : set up wetting/drying '
291         WRITE(numout,*) '~~~~~~~ '
292         WRITE(numout,*) '   Namelist namwad'
293         WRITE(numout,*) '      turn on(off) wetting/drying     ln_wad        = ', ln_wad
294         WRITE(numout,*) '      minimal dry cell depth          rn_ wadmin    = ', rn_wadmin 
295         WRITE(numout,*) '      tollerance of minimual depth    rn_wadmine    = ', rn_wadmine
296         WRITE(numout,*) ' threshold elevation for land masking rn_landele    = ', rn_landele
297         WRITE(numout,*) ' Maximum number of iteration for W/D limiter        = ', nn_waditr
298      ENDIF
299
300   END SUBROUTINE dom_wad
301
302
303   SUBROUTINE dom_ctl
304      !!----------------------------------------------------------------------
305      !!                     ***  ROUTINE dom_ctl  ***
306      !!
307      !! ** Purpose :   Domain control.
308      !!
309      !! ** Method  :   compute and print extrema of masked scale factors
310      !!----------------------------------------------------------------------
311      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
312      INTEGER, DIMENSION(2) ::   iloc   !
313      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
314      !!----------------------------------------------------------------------
315      !
316      IF(lk_mpp) THEN
317         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
318         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
319         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
320         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
321      ELSE
322         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )   
323         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )   
324         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )   
325         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )   
326
327         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
328         iimi1 = iloc(1) + nimpp - 1
329         ijmi1 = iloc(2) + njmpp - 1
330         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
331         iimi2 = iloc(1) + nimpp - 1
332         ijmi2 = iloc(2) + njmpp - 1
333         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
334         iima1 = iloc(1) + nimpp - 1
335         ijma1 = iloc(2) + njmpp - 1
336         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
337         iima2 = iloc(1) + nimpp - 1
338         ijma2 = iloc(2) + njmpp - 1
339      ENDIF
340      IF(lwp) THEN
341         WRITE(numout,*)
342         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
343         WRITE(numout,*) '~~~~~~~'
344         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
345         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
346         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
347         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
348      ENDIF
349      !
350   END SUBROUTINE dom_ctl
351
352   SUBROUTINE dom_stiff
353      !!----------------------------------------------------------------------
354      !!                  ***  ROUTINE dom_stiff  ***
355      !!                     
356      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
357      !!
358      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
359      !!                Save the maximum in the vertical direction
360      !!                (this number is only relevant in s-coordinates)
361      !!
362      !!                Haney, R. L., 1991: On the pressure gradient force
363      !!                over steep topography in sigma coordinate ocean models.
364      !!                J. Phys. Oceanogr., 21, 610???619.
365      !!----------------------------------------------------------------------
366      INTEGER  ::   ji, jj, jk 
367      REAL(wp) ::   zrxmax
368      REAL(wp), DIMENSION(4) :: zr1
369      !!----------------------------------------------------------------------
370      rx1(:,:) = 0.e0
371      zrxmax   = 0.e0
372      zr1(:)   = 0.e0
373     
374      DO ji = 2, jpim1
375         DO jj = 2, jpjm1
376            DO jk = 1, jpkm1
377               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji-1,jj  ,jk  )  & 
378                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1)) &
379                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji-1,jj  ,jk  )  &
380                    &                         -gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1) + rsmall) )
381               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw(ji+1,jj  ,jk  )-gdepw(ji  ,jj  ,jk  )  &
382                    &                         +gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
383                    &                        /(gdepw(ji+1,jj  ,jk  )+gdepw(ji  ,jj  ,jk  )  &
384                    &                         -gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
385               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw(ji  ,jj+1,jk  )-gdepw(ji  ,jj  ,jk  )  &
386                    &                         +gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
387                    &                        /(gdepw(ji  ,jj+1,jk  )+gdepw(ji  ,jj  ,jk  )  &
388                    &                         -gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
389               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji  ,jj-1,jk  )  &
390                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1)) &
391                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji  ,jj-1,jk  )  &
392                    &                         -gdepw(ji,  jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1) + rsmall) )
393               zrxmax = MAXVAL(zr1(1:4))
394               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
395            END DO
396         END DO
397      END DO
398
399      CALL lbc_lnk( rx1, 'T', 1. )
400
401      zrxmax = MAXVAL(rx1)
402
403      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
404
405      IF(lwp) THEN
406         WRITE(numout,*)
407         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
408         WRITE(numout,*) '~~~~~~~~~'
409      ENDIF
410
411   END SUBROUTINE dom_stiff
412
413
414
415   !!======================================================================
416END MODULE domain
Note: See TracBrowser for help on using the repository browser.