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

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 3865

Last change on this file since 3865 was 3865, checked in by acc, 11 years ago

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