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

source: branches/2013/dev_r3940_CNRS4_IOCRS/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 4015

Last change on this file since 4015 was 4015, checked in by cetlod, 11 years ago

2013/dev_r3940_CNRS4_IOCRS: 1st step, add new routines for outputs coarsening

  • Property svn:keywords set to Id
File size: 18.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   !!            3.6  !  2013     ( J. Simeon, C. Calone, G. Madec, C. Ethe ) Online coarsening of outputs
15   !!----------------------------------------------------------------------
16   
17   !!----------------------------------------------------------------------
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   !!----------------------------------------------------------------------
22   USE oce             ! ocean variables
23   USE dom_oce         ! domain: ocean
24   USE sbc_oce         ! surface boundary condition: ocean
25   USE phycst          ! physical constants
26   USE closea          ! closed seas
27   USE in_out_manager  ! I/O manager
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
35   USE domvvl          ! variable volume
36   USE c1d             ! 1D vertical configuration
37   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
38   USE timing          ! Timing
39   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   dom_init   ! called by opa.F90
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48   !!-------------------------------------------------------------------------
49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
50   !! $Id$
51   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
52   !!-------------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE dom_init
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE dom_init  ***
58      !!                   
59      !! ** Purpose :   Domain initialization. Call the routines that are
60      !!              required to create the arrays which define the space
61      !!              and time domain of the ocean model.
62      !!
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
69      !!              - 1D configuration, move Coriolis, u and v at T-point
70      !!----------------------------------------------------------------------
71      INTEGER ::   jk          ! dummy loop argument
72      INTEGER ::   iconf = 0   ! local integers
73      !!----------------------------------------------------------------------
74      !
75      IF( nn_timing == 1 )   CALL timing_start('dom_init')
76      !
77      IF(lwp) THEN
78         WRITE(numout,*)
79         WRITE(numout,*) 'dom_init : domain initialization'
80         WRITE(numout,*) '~~~~~~~~'
81      ENDIF
82      !
83                             CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
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 , ln_crs
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      !
207      IF(lwp) THEN
208         WRITE(numout,*)
209         WRITE(numout,*) '   Namelist namdom : space & time domain'
210         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
211         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
212         WRITE(numout,*) '      min number of ocean level (<0)       '
213         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
214         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
215         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
216         WRITE(numout,*) '           = 0   no file created           '
217         WRITE(numout,*) '           = 1   mesh_mask                 '
218         WRITE(numout,*) '           = 2   mesh and mask             '
219         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
220         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
221         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
222         WRITE(numout,*) '      time-splitting: nb of sub time-step   nn_baro   = ', nn_baro
223         WRITE(numout,*) '      acceleration of converge              nn_acc    = ', nn_acc
224         WRITE(numout,*) '        nn_acc=1: surface tracer rdt        rn_rdtmin = ', rn_rdtmin
225         WRITE(numout,*) '                  bottom  tracer rdt        rdtmax    = ', rn_rdtmax
226         WRITE(numout,*) '                  depth of transition       rn_rdth   = ', rn_rdth
227         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
228         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
229      ENDIF
230
231      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
232      e3zps_min = rn_e3zps_min
233      e3zps_rat = rn_e3zps_rat
234      nmsh      = nn_msh
235      nacc      = nn_acc
236      atfp      = rn_atfp
237      rdt       = rn_rdt
238      rdtmin    = rn_rdtmin
239      rdtmax    = rn_rdtmin
240      rdth      = rn_rdth
241
242      REWIND( numnam )              ! Namelist cross land advection
243      READ  ( numnam, namcla )
244      IF(lwp) THEN
245         WRITE(numout,*)
246         WRITE(numout,*) '   Namelist namcla'
247         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
248      ENDIF
249
250#if defined key_netcdf4
251      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
252      REWIND( numnam )                    ! Namelist namnc4 : netcdf4 chunking parameters
253      READ  ( numnam, namnc4 )
254      IF(lwp) THEN                        ! control print
255         WRITE(numout,*)
256         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
257         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
258         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
259         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
260         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
261      ENDIF
262
263      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
264      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
265      snc4set%ni   = nn_nchunks_i
266      snc4set%nj   = nn_nchunks_j
267      snc4set%nk   = nn_nchunks_k
268      snc4set%luse = ln_nc4zip
269#else
270      snc4set%luse = .FALSE.        ! No NetCDF 4 case
271#endif
272      !
273   END SUBROUTINE dom_nam
274
275
276   SUBROUTINE dom_ctl
277      !!----------------------------------------------------------------------
278      !!                     ***  ROUTINE dom_ctl  ***
279      !!
280      !! ** Purpose :   Domain control.
281      !!
282      !! ** Method  :   compute and print extrema of masked scale factors
283      !!----------------------------------------------------------------------
284      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
285      INTEGER, DIMENSION(2) ::   iloc   !
286      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
287      !!----------------------------------------------------------------------
288      !
289      IF(lk_mpp) THEN
290         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
291         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
292         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
293         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
294      ELSE
295         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )   
296         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )   
297         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )   
298         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )   
299
300         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
301         iimi1 = iloc(1) + nimpp - 1
302         ijmi1 = iloc(2) + njmpp - 1
303         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
304         iimi2 = iloc(1) + nimpp - 1
305         ijmi2 = iloc(2) + njmpp - 1
306         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
307         iima1 = iloc(1) + nimpp - 1
308         ijma1 = iloc(2) + njmpp - 1
309         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
310         iima2 = iloc(1) + nimpp - 1
311         ijma2 = iloc(2) + njmpp - 1
312      ENDIF
313      IF(lwp) THEN
314         WRITE(numout,*)
315         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
316         WRITE(numout,*) '~~~~~~~'
317         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
318         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
319         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
320         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
321      ENDIF
322      !
323   END SUBROUTINE dom_ctl
324
325   SUBROUTINE dom_stiff
326      !!----------------------------------------------------------------------
327      !!                  ***  ROUTINE dom_stiff  ***
328      !!                     
329      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
330      !!
331      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
332      !!                Save the maximum in the vertical direction
333      !!                (this number is only relevant in s-coordinates)
334      !!
335      !!                Haney, R. L., 1991: On the pressure gradient force
336      !!                over steep topography in sigma coordinate ocean models.
337      !!                J. Phys. Oceanogr., 21, 610???619.
338      !!----------------------------------------------------------------------
339      INTEGER  ::   ji, jj, jk 
340      REAL(wp) ::   zrxmax
341      REAL(wp), DIMENSION(4) :: zr1
342      !!----------------------------------------------------------------------
343      rx1(:,:) = 0.e0
344      zrxmax   = 0.e0
345      zr1(:)   = 0.e0
346     
347      DO ji = 2, jpim1
348         DO jj = 2, jpjm1
349            DO jk = 1, jpkm1
350               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji-1,jj  ,jk  )  & 
351                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1)) &
352                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji-1,jj  ,jk  )  &
353                    &                         -gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1) + rsmall) )
354               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw(ji+1,jj  ,jk  )-gdepw(ji  ,jj  ,jk  )  &
355                    &                         +gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
356                    &                        /(gdepw(ji+1,jj  ,jk  )+gdepw(ji  ,jj  ,jk  )  &
357                    &                         -gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
358               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw(ji  ,jj+1,jk  )-gdepw(ji  ,jj  ,jk  )  &
359                    &                         +gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
360                    &                        /(gdepw(ji  ,jj+1,jk  )+gdepw(ji  ,jj  ,jk  )  &
361                    &                         -gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
362               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji  ,jj-1,jk  )  &
363                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1)) &
364                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji  ,jj-1,jk  )  &
365                    &                         -gdepw(ji,  jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1) + rsmall) )
366               zrxmax = MAXVAL(zr1(1:4))
367               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
368            END DO
369         END DO
370      END DO
371
372      CALL lbc_lnk( rx1, 'T', 1. )
373
374      zrxmax = MAXVAL(rx1)
375
376      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
377
378      IF(lwp) THEN
379         WRITE(numout,*)
380         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
381         WRITE(numout,*) '~~~~~~~~~'
382      ENDIF
383
384   END SUBROUTINE dom_stiff
385
386
387
388   !!======================================================================
389END MODULE domain
Note: See TracBrowser for help on using the repository browser.