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

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 3764

Last change on this file since 3764 was 3764, checked in by smasson, 11 years ago

dev_MERGE_2012: report bugfixes done in the trunk from r3555 to r3763 into dev_MERGE_2012
  • Property svn:keywords set to Id
File size: 18.7 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_clo      ! Closed seas and lake
84                             CALL dom_hgr      ! Horizontal mesh
85                             CALL dom_zgr      ! Vertical mesh and bathymetry
86                             CALL dom_msk      ! Masks
87      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
88      IF( lk_vvl         )   CALL dom_vvl      ! Vertical variable mesh
89      !
90      IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
91      !
92      hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points
93      hv(:,:) = 0._wp
94      DO jk = 1, jpk
95         hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
96         hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk)
97      END DO
98      !                                        ! Inverse of the local depth
99      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
100      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
101
102                             CALL dom_stp      ! time step
103      IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file
104      IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control
105      !
106      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
107      !
108   END SUBROUTINE dom_init
109
110
111   SUBROUTINE dom_nam
112      !!----------------------------------------------------------------------
113      !!                     ***  ROUTINE dom_nam  ***
114      !!                   
115      !! ** Purpose :   read domaine namelists and print the variables.
116      !!
117      !! ** input   : - namrun namelist
118      !!              - namdom namelist
119      !!              - namcla namelist
120      !!              - namnc4 namelist   ! "key_netcdf4" only
121      !!----------------------------------------------------------------------
122      USE ioipsl
123      NAMELIST/namrun/ nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
124         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
125         &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz
126      NAMELIST/namdom/ nn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin,   &
127         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,            &
128         &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea
129      NAMELIST/namcla/ nn_cla
130#if defined key_netcdf4
131      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
132#endif
133      !!----------------------------------------------------------------------
134
135      REWIND( numnam )              ! Namelist namrun : parameters of the run
136      READ  ( numnam, namrun )
137      !
138      IF(lwp) THEN                  ! control print
139         WRITE(numout,*)
140         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
141         WRITE(numout,*) '~~~~~~~ '
142         WRITE(numout,*) '   Namelist namrun'
143         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
144         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
145         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
146         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
147         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
148         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
149         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
150         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
151         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
152         WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
153         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
154         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
155         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
156         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
157         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
158      ENDIF
159
160      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
161      cexper = cn_exp
162      nrstdt = nn_rstctl
163      nit000 = nn_it000
164      nitend = nn_itend
165      ndate0 = nn_date0
166      nleapy = nn_leapy
167      ninist = nn_istate
168      nstock = nn_stock
169      nwrite = nn_write
170
171
172      !                             ! control of output frequency
173      IF ( nstock == 0 .OR. nstock > nitend ) THEN
174         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
175         CALL ctl_warn( ctmp1 )
176         nstock = nitend
177      ENDIF
178      IF ( nwrite == 0 ) THEN
179         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
180         CALL ctl_warn( ctmp1 )
181         nwrite = nitend
182      ENDIF
183
184#if defined key_agrif
185      IF( Agrif_Root() ) THEN
186#endif
187      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
188      CASE (  1 ) 
189         CALL ioconf_calendar('gregorian')
190         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
191      CASE (  0 )
192         CALL ioconf_calendar('noleap')
193         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
194      CASE ( 30 )
195         CALL ioconf_calendar('360d')
196         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
197      END SELECT
198#if defined key_agrif
199      ENDIF
200#endif
201
202      REWIND( numnam )              ! Namelist namdom : space & time domain (bathymetry, mesh, timestep)
203      READ  ( numnam, namdom )
204
205      IF(lwp) THEN
206         WRITE(numout,*)
207         WRITE(numout,*) '   Namelist namdom : space & time domain'
208         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
209         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
210         WRITE(numout,*) '      min number of ocean level (<0)       '
211         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
212         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
213         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
214         WRITE(numout,*) '           = 0   no file created           '
215         WRITE(numout,*) '           = 1   mesh_mask                 '
216         WRITE(numout,*) '           = 2   mesh and mask             '
217         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
218         WRITE(numout,*) '      ocean time step                      rn_rdt    = ', rn_rdt
219         WRITE(numout,*) '      asselin time filter parameter        rn_atfp   = ', rn_atfp
220         WRITE(numout,*) '      time-splitting: nb of sub time-step  nn_baro   = ', nn_baro
221         WRITE(numout,*) '      acceleration of converge             nn_acc    = ', nn_acc
222         WRITE(numout,*) '        nn_acc=1: surface tracer rdt       rn_rdtmin = ', rn_rdtmin
223         WRITE(numout,*) '                  bottom  tracer rdt       rdtmax    = ', rn_rdtmax
224         WRITE(numout,*) '                  depth of transition      rn_rdth   = ', rn_rdth
225         WRITE(numout,*) '      suppression of closed seas (=0)      nn_closea = ', nn_closea
226      ENDIF
227
228      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
229      e3zps_min = rn_e3zps_min
230      e3zps_rat = rn_e3zps_rat
231      nmsh      = nn_msh
232      nacc      = nn_acc
233      atfp      = rn_atfp
234      rdt       = rn_rdt
235      rdtmin    = rn_rdtmin
236      rdtmax    = rn_rdtmin
237      rdth      = rn_rdth
238
239      REWIND( numnam )              ! Namelist cross land advection
240      READ  ( numnam, namcla )
241      IF(lwp) THEN
242         WRITE(numout,*)
243         WRITE(numout,*) '   Namelist namcla'
244         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
245      ENDIF
246
247#if defined key_netcdf4
248      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
249      REWIND( numnam )                    ! Namelist namnc4 : netcdf4 chunking parameters
250      READ  ( numnam, namnc4 )
251      IF(lwp) THEN                        ! control print
252         WRITE(numout,*)
253         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
254         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
255         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
256         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
257         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
258      ENDIF
259
260      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
261      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
262      snc4set%ni   = nn_nchunks_i
263      snc4set%nj   = nn_nchunks_j
264      snc4set%nk   = nn_nchunks_k
265      snc4set%luse = ln_nc4zip
266#else
267      snc4set%luse = .FALSE.        ! No NetCDF 4 case
268#endif
269      !
270   END SUBROUTINE dom_nam
271
272
273   SUBROUTINE dom_ctl
274      !!----------------------------------------------------------------------
275      !!                     ***  ROUTINE dom_ctl  ***
276      !!
277      !! ** Purpose :   Domain control.
278      !!
279      !! ** Method  :   compute and print extrema of masked scale factors
280      !!----------------------------------------------------------------------
281      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
282      INTEGER, DIMENSION(2) ::   iloc   !
283      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
284      !!----------------------------------------------------------------------
285      !
286      IF(lk_mpp) THEN
287         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
288         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
289         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
290         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
291      ELSE
292         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )   
293         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )   
294         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )   
295         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )   
296
297         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
298         iimi1 = iloc(1) + nimpp - 1
299         ijmi1 = iloc(2) + njmpp - 1
300         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
301         iimi2 = iloc(1) + nimpp - 1
302         ijmi2 = iloc(2) + njmpp - 1
303         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
304         iima1 = iloc(1) + nimpp - 1
305         ijma1 = iloc(2) + njmpp - 1
306         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
307         iima2 = iloc(1) + nimpp - 1
308         ijma2 = iloc(2) + njmpp - 1
309      ENDIF
310      IF(lwp) THEN
311         WRITE(numout,*)
312         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
313         WRITE(numout,*) '~~~~~~~'
314         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
315         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
316         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
317         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
318      ENDIF
319      !
320   END SUBROUTINE dom_ctl
321
322   SUBROUTINE dom_stiff
323      !!----------------------------------------------------------------------
324      !!                  ***  ROUTINE dom_stiff  ***
325      !!                     
326      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
327      !!
328      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
329      !!                Save the maximum in the vertical direction
330      !!                (this number is only relevant in s-coordinates)
331      !!
332      !!                Haney, R. L., 1991: On the pressure gradient force
333      !!                over steep topography in sigma coordinate ocean models.
334      !!                J. Phys. Oceanogr., 21, 610???619.
335      !!----------------------------------------------------------------------
336      INTEGER  ::   ji, jj, jk 
337      REAL(wp) ::   zrxmax
338      REAL(wp), DIMENSION(4) :: zr1
339      !!----------------------------------------------------------------------
340      rx1(:,:) = 0.e0
341      zrxmax   = 0.e0
342      zr1(:)   = 0.e0
343     
344      DO ji = 2, jpim1
345         DO jj = 2, jpjm1
346            DO jk = 1, jpkm1
347               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji-1,jj  ,jk  )  & 
348                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1)) &
349                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji-1,jj  ,jk  )  &
350                    &                         -gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1) + rsmall) )
351               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw(ji+1,jj  ,jk  )-gdepw(ji  ,jj  ,jk  )  &
352                    &                         +gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
353                    &                        /(gdepw(ji+1,jj  ,jk  )+gdepw(ji  ,jj  ,jk  )  &
354                    &                         -gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
355               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw(ji  ,jj+1,jk  )-gdepw(ji  ,jj  ,jk  )  &
356                    &                         +gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
357                    &                        /(gdepw(ji  ,jj+1,jk  )+gdepw(ji  ,jj  ,jk  )  &
358                    &                         -gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
359               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji  ,jj-1,jk  )  &
360                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1)) &
361                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji  ,jj-1,jk  )  &
362                    &                         -gdepw(ji,  jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1) + rsmall) )
363               zrxmax = MAXVAL(zr1(1:4))
364               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
365            END DO
366         END DO
367      END DO
368
369      CALL lbc_lnk( rx1, 'T', 1. )
370
371      zrxmax = MAXVAL(rx1)
372
373      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
374
375      IF(lwp) THEN
376         WRITE(numout,*)
377         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
378         WRITE(numout,*) '~~~~~~~~~'
379      ENDIF
380
381   END SUBROUTINE dom_stiff
382
383
384
385   !!======================================================================
386END MODULE domain
Note: See TracBrowser for help on using the repository browser.