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
RevLine 
[3]1MODULE domain
2   !!==============================================================================
3   !!                       ***  MODULE domain   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
[1438]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
[2528]13   !!            3.3  !  2010-11  (G. Madec)  initialisation in C1D configuration
[4015]14   !!            3.6  !  2013     ( J. Simeon, C. Calone, G. Madec, C. Ethe ) Online coarsening of outputs
[3]15   !!----------------------------------------------------------------------
[1438]16   
17   !!----------------------------------------------------------------------
[3]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   !!----------------------------------------------------------------------
[2528]22   USE oce             ! ocean variables
23   USE dom_oce         ! domain: ocean
[888]24   USE sbc_oce         ! surface boundary condition: ocean
[719]25   USE phycst          ! physical constants
[1601]26   USE closea          ! closed seas
[719]27   USE in_out_manager  ! I/O manager
[3]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
[592]35   USE domvvl          ! variable volume
[2528]36   USE c1d             ! 1D vertical configuration
37   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
[3294]38   USE timing          ! Timing
[3680]39   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
[3]40
41   IMPLICIT NONE
42   PRIVATE
43
[1438]44   PUBLIC   dom_init   ! called by opa.F90
[3]45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
[1438]48   !!-------------------------------------------------------------------------
[2528]49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]50   !! $Id$
[2528]51   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
[1438]52   !!-------------------------------------------------------------------------
[3]53CONTAINS
54
55   SUBROUTINE dom_init
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE dom_init  ***
58      !!                   
59      !! ** Purpose :   Domain initialization. Call the routines that are
[1601]60      !!              required to create the arrays which define the space
61      !!              and time domain of the ocean model.
[3]62      !!
[1601]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
[2528]69      !!              - 1D configuration, move Coriolis, u and v at T-point
[3]70      !!----------------------------------------------------------------------
[3764]71      INTEGER ::   jk          ! dummy loop argument
72      INTEGER ::   iconf = 0   ! local integers
[3]73      !!----------------------------------------------------------------------
[1601]74      !
[3764]75      IF( nn_timing == 1 )   CALL timing_start('dom_init')
[3294]76      !
[3]77      IF(lwp) THEN
78         WRITE(numout,*)
79         WRITE(numout,*) 'dom_init : domain initialization'
80         WRITE(numout,*) '~~~~~~~~'
81      ENDIF
[1601]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
[3680]88      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
[1601]89      IF( lk_vvl         )   CALL dom_vvl      ! Vertical variable mesh
90      !
[3764]91      IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
[2528]92      !
[3764]93      hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points
94      hv(:,:) = 0._wp
[3]95      DO jk = 1, jpk
96         hu(:,:) = hu(:,:) + fse3u(:,:,jk) * umask(:,:,jk)
97         hv(:,:) = hv(:,:) + fse3v(:,:,jk) * vmask(:,:,jk)
98      END DO
[1601]99      !                                        ! Inverse of the local depth
[3764]100      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
101      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
[216]102
[1601]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
[1438]106      !
[3764]107      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
[3294]108      !
[3]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
[2528]121      !!              - namnc4 namelist   ! "key_netcdf4" only
[3]122      !!----------------------------------------------------------------------
123      USE ioipsl
[1601]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
[2528]127      NAMELIST/namdom/ nn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin,   &
128         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,            &
[4015]129         &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea , ln_crs
[1601]130      NAMELIST/namcla/ nn_cla
[2528]131#if defined key_netcdf4
132      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
133#endif
[3]134      !!----------------------------------------------------------------------
135
[1601]136      REWIND( numnam )              ! Namelist namrun : parameters of the run
137      READ  ( numnam, namrun )
138      !
139      IF(lwp) THEN                  ! control print
[3]140         WRITE(numout,*)
141         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
142         WRITE(numout,*) '~~~~~~~ '
[1601]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
[1604]147         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
[1601]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
[3]159      ENDIF
160
[1601]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
[3]171
[1601]172
173      !                             ! control of output frequency
[1335]174      IF ( nstock == 0 .OR. nstock > nitend ) THEN
[1601]175         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
[783]176         CALL ctl_warn( ctmp1 )
[1335]177         nstock = nitend
[3]178      ENDIF
179      IF ( nwrite == 0 ) THEN
[1601]180         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
[783]181         CALL ctl_warn( ctmp1 )
182         nwrite = nitend
[3]183      ENDIF
184
[2528]185#if defined key_agrif
[1601]186      IF( Agrif_Root() ) THEN
[2528]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
[1601]200      ENDIF
[2528]201#endif
[3]202
[2528]203      REWIND( numnam )              ! Namelist namdom : space & time domain (bathymetry, mesh, timestep)
[3]204      READ  ( numnam, namdom )
[4015]205 
206      !
[3]207      IF(lwp) THEN
[72]208         WRITE(numout,*)
[1601]209         WRITE(numout,*) '   Namelist namdom : space & time domain'
210         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
[2528]211         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
212         WRITE(numout,*) '      min number of ocean level (<0)       '
[1601]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
[2528]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'
[4015]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
[223]229      ENDIF
230
[1601]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
[2528]242      REWIND( numnam )              ! Namelist cross land advection
[3]243      READ  ( numnam, namcla )
244      IF(lwp) THEN
[72]245         WRITE(numout,*)
[1601]246         WRITE(numout,*) '   Namelist namcla'
247         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
[3]248      ENDIF
249
[2528]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
[1601]262
[2528]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
[1438]272      !
[3]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
[1601]285      INTEGER, DIMENSION(2) ::   iloc   !
[3]286      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
287      !!----------------------------------------------------------------------
[1601]288      !
289      IF(lk_mpp) THEN
[181]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
[3764]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 )   
[32]299
[3764]300         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
[181]301         iimi1 = iloc(1) + nimpp - 1
302         ijmi1 = iloc(2) + njmpp - 1
[3764]303         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
[181]304         iimi2 = iloc(1) + nimpp - 1
305         ijmi2 = iloc(2) + njmpp - 1
[3764]306         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp )
[181]307         iima1 = iloc(1) + nimpp - 1
308         ijma1 = iloc(2) + njmpp - 1
[3764]309         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp )
[181]310         iima2 = iloc(1) + nimpp - 1
311         ijma2 = iloc(2) + njmpp - 1
[32]312      ENDIF
[3]313      IF(lwp) THEN
[1601]314         WRITE(numout,*)
315         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
316         WRITE(numout,*) '~~~~~~~'
[181]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
[3]321      ENDIF
[1438]322      !
[3]323   END SUBROUTINE dom_ctl
324
[3680]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
[3]388   !!======================================================================
389END MODULE domain
Note: See TracBrowser for help on using the repository browser.