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

source: branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 3970

Last change on this file since 3970 was 3970, checked in by cbricaud, 11 years ago

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