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

source: branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 3600

Last change on this file since 3600 was 3600, checked in by rfurner, 12 years ago

Changes from branch dev_r3435_UKMO7_SCOORDS revision 3435 to 3507 merged in

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