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

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef_tools/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90 @ 6827

Last change on this file since 6827 was 6827, checked in by flavoni, 8 years ago

#1692 - branch SIMPLIF_2_usrdef: commit routines Fortran to create domain_cfg.nc file, to have domain (input) files for new simplification branch. To be moved in TOOLS directory
  • Property svn:keywords set to Id
File size: 31.7 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   !!            3.7  !  2015-11  (G. Madec, A. Coward)  time varying zgr by default
16   !!----------------------------------------------------------------------
17   
18   !!----------------------------------------------------------------------
19   !!   dom_init       : initialize the space and time domain
20   !!   dom_nam        : read and contral domain namelists
21   !!   dom_ctl        : control print for the ocean domain
22   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean variables
25   USE dom_oce         ! domain: ocean
26   USE phycst          ! physical constants
27   USE closea          ! closed seas
28   USE domhgr          ! domain: set the horizontal mesh
29   USE domzgr          ! domain: set the vertical mesh
30   USE domstp          ! domain: set the time-step
31   USE dommsk          ! domain: set the mask system
32   USE domwri          ! domain: write the meshmask file
33   USE domvvl          ! variable volume
34   !
35   USE in_out_manager  ! I/O manager
36   USE iom             !
37   USE wrk_nemo        ! Memory Allocation
38   USE lib_mpp         ! distributed memory computing library
39   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
40   USE timing          ! Timing
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   dom_init   ! called by opa.F90
46
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 indices
71      INTEGER ::   iconf = 0   ! local integers
72      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_hu_0, z1_hv_0
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      !                       !==  Reference coordinate system  ==!
84      !
85                     CALL dom_nam               ! read namelist ( namrun, namdom )
86                     CALL dom_clo               ! Closed seas and lake
87                     CALL dom_hgr               ! Horizontal mesh
88                     CALL dom_zgr               ! Vertical mesh and bathymetry
89                     CALL dom_msk               ! Masks
90      IF( ln_sco )   CALL dom_stiff             ! Maximum stiffness ratio/hydrostatic consistency
91      !
92      ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1)   ! Reference ocean thickness
93      hu_0(:,:) = e3u_0(:,:,1) * umask(:,:,1)
94      hv_0(:,:) = e3v_0(:,:,1) * vmask(:,:,1)
95      DO jk = 2, jpk
96         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
97         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
98         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
99      END DO
100      !
101      !              !==  time varying part of coordinate system  ==!
102      !
103      IF( ln_linssh ) THEN          ! Fix in time : set to the reference one for all
104         !       before        !          now          !       after         !
105         ;  gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
106         ;  gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
107         ;                     ;   gde3w_n = gde3w_0   !        ---          !
108         !                                                                 
109         ;    e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
110         ;    e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
111         ;    e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
112         ;                     ;     e3f_n =   e3f_0   !        ---          !
113         ;    e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
114         ;   e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
115         ;   e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
116         !
117         CALL wrk_alloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
118         !
119         z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
120         z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
121         !
122         !        before       !          now          !       after         !
123         ;                     ;      ht_n =    ht_0   !                     ! water column thickness
124         ;     hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   !
125         ;     hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   !
126         ;  r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness
127         ;  r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   !
128         !
129         CALL wrk_dealloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
130         !
131      ELSE                         ! time varying : initialize before/now/after variables
132         !
133         CALL dom_vvl_init 
134         !
135      ENDIF
136      !
137      CALL cfg_write         ! create the configuration file
138      !
139      IF( nn_timing == 1 )   CALL timing_stop('dom_init')
140      !
141   END SUBROUTINE dom_init
142
143
144   SUBROUTINE dom_nam
145      !!----------------------------------------------------------------------
146      !!                     ***  ROUTINE dom_nam  ***
147      !!                   
148      !! ** Purpose :   read domaine namelists and print the variables.
149      !!
150      !! ** input   : - namrun namelist
151      !!              - namdom namelist
152      !!              - namnc4 namelist   ! "key_netcdf4" only
153      !!----------------------------------------------------------------------
154      USE ioipsl
155      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 &
156                       nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     &
157         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
158         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     &
159         &             ln_cfmeta, ln_iscpl
160      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, &
161         &             rn_atfp , rn_rdt   , nn_closea   , ln_crs      , jphgr_msh ,                  &
162         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         &
163         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  &
164         &             ppa2, ppkth2, ppacr2
165#if defined key_netcdf4
166      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
167#endif
168      INTEGER  ::   ios                 ! Local integer output status for namelist read
169      !!----------------------------------------------------------------------
170
171      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
172      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
173901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
174
175      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
176      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
177902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
178      IF(lwm) WRITE ( numond, namrun )
179      !
180      IF(lwp) THEN                  ! control print
181         WRITE(numout,*)
182         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
183         WRITE(numout,*) '~~~~~~~ '
184         WRITE(numout,*) '   Namelist namrun'
185         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
186         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
187         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in
188         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir
189         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out
190         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir
191         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
192         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler
193         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
194         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
195         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
196         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
197         WRITE(numout,*) '      initial time of day in hhmm     nn_time0   = ', nn_time0
198         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
199         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
200         IF( ln_rst_list ) THEN
201            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist
202         ELSE
203            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
204         ENDIF
205         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
206         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
207         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
208         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
209         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
210         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl   = ', ln_iscpl
211      ENDIF
212
213      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
214      cexper = cn_exp
215      nrstdt = nn_rstctl
216      nit000 = nn_it000
217      nitend = nn_itend
218      ndate0 = nn_date0
219      nleapy = nn_leapy
220      ninist = nn_istate
221      nstock = nn_stock
222      nstocklist = nn_stocklist
223      nwrite = nn_write
224      neuler = nn_euler
225      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN
226         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 '
227         CALL ctl_warn( ctmp1 )
228         neuler = 0
229      ENDIF
230
231      !                             ! control of output frequency
232      IF ( nstock == 0 .OR. nstock > nitend ) THEN
233         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
234         CALL ctl_warn( ctmp1 )
235         nstock = nitend
236      ENDIF
237      IF ( nwrite == 0 ) THEN
238         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
239         CALL ctl_warn( ctmp1 )
240         nwrite = nitend
241      ENDIF
242
243#if defined key_agrif
244      IF( Agrif_Root() ) THEN
245#endif
246      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
247      CASE (  1 ) 
248         CALL ioconf_calendar('gregorian')
249         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
250      CASE (  0 )
251         CALL ioconf_calendar('noleap')
252         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
253      CASE ( 30 )
254         CALL ioconf_calendar('360d')
255         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
256      END SELECT
257#if defined key_agrif
258      ENDIF
259#endif
260
261      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
262      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
263903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
264 
265      !
266      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
267      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
268904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
269      IF(lwm) WRITE ( numond, namdom )
270      !
271      IF(lwp) THEN
272         WRITE(numout,*)
273         WRITE(numout,*) '   Namelist namdom : space & time domain'
274         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
275         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
276         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
277         WRITE(numout,*) '      min number of ocean level (<0)       '
278         WRITE(numout,*) '      treshold to open the isf cavity   rn_isfhmin   = ', rn_isfhmin, ' (m)'
279         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
280         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
281         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
282         WRITE(numout,*) '           = 0   no file created           '
283         WRITE(numout,*) '           = 1   mesh_mask                 '
284         WRITE(numout,*) '           = 2   mesh and mask             '
285         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask'
286         WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt
287         WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp
288         WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea
289         WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs
290         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
291         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
292         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
293         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
294         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
295         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
296         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
297         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
298         WRITE(numout,*) '                                        ppa0            = ', ppa0
299         WRITE(numout,*) '                                        ppa1            = ', ppa1
300         WRITE(numout,*) '                                        ppkth           = ', ppkth
301         WRITE(numout,*) '                                        ppacr           = ', ppacr
302         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
303         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
304         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
305         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
306         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
307         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
308      ENDIF
309      !
310      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
311      e3zps_min = rn_e3zps_min
312      e3zps_rat = rn_e3zps_rat
313      nmsh      = nn_msh
314      atfp      = rn_atfp
315      rdt       = rn_rdt
316
317#if defined key_netcdf4
318      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
319      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
320      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
321907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
322
323      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
324      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
325908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
326      IF(lwm) WRITE( numond, namnc4 )
327
328      IF(lwp) THEN                        ! control print
329         WRITE(numout,*)
330         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
331         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
332         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
333         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
334         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
335      ENDIF
336
337      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
338      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
339      snc4set%ni   = nn_nchunks_i
340      snc4set%nj   = nn_nchunks_j
341      snc4set%nk   = nn_nchunks_k
342      snc4set%luse = ln_nc4zip
343#else
344      snc4set%luse = .FALSE.        ! No NetCDF 4 case
345#endif
346      !
347   END SUBROUTINE dom_nam
348
349
350   SUBROUTINE dom_ctl
351      !!----------------------------------------------------------------------
352      !!                     ***  ROUTINE dom_ctl  ***
353      !!
354      !! ** Purpose :   Domain control.
355      !!
356      !! ** Method  :   compute and print extrema of masked scale factors
357      !!----------------------------------------------------------------------
358      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
359      INTEGER, DIMENSION(2) ::   iloc   !
360      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
361      !!----------------------------------------------------------------------
362      !
363      IF(lk_mpp) THEN
364         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
365         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
366         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
367         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
368      ELSE
369         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
370         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
371         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )   
372         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )   
373
374         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
375         iimi1 = iloc(1) + nimpp - 1
376         ijmi1 = iloc(2) + njmpp - 1
377         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
378         iimi2 = iloc(1) + nimpp - 1
379         ijmi2 = iloc(2) + njmpp - 1
380         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
381         iima1 = iloc(1) + nimpp - 1
382         ijma1 = iloc(2) + njmpp - 1
383         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
384         iima2 = iloc(1) + nimpp - 1
385         ijma2 = iloc(2) + njmpp - 1
386      ENDIF
387      IF(lwp) THEN
388         WRITE(numout,*)
389         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
390         WRITE(numout,*) '~~~~~~~'
391         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
392         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
393         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
394         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
395      ENDIF
396      !
397   END SUBROUTINE dom_ctl
398
399
400   SUBROUTINE dom_stiff
401      !!----------------------------------------------------------------------
402      !!                  ***  ROUTINE dom_stiff  ***
403      !!                     
404      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
405      !!
406      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
407      !!                Save the maximum in the vertical direction
408      !!                (this number is only relevant in s-coordinates)
409      !!
410      !!                Haney, R. L., 1991: On the pressure gradient force
411      !!                over steep topography in sigma coordinate ocean models.
412      !!                J. Phys. Oceanogr., 21, 610???619.
413      !!----------------------------------------------------------------------
414      INTEGER  ::   ji, jj, jk 
415      REAL(wp) ::   zrxmax
416      REAL(wp), DIMENSION(4) ::   zr1
417      !!----------------------------------------------------------------------
418      rx1(:,:) = 0._wp
419      zrxmax   = 0._wp
420      zr1(:)   = 0._wp
421      !
422      DO ji = 2, jpim1
423         DO jj = 2, jpjm1
424            DO jk = 1, jpkm1
425               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
426                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
427                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
428                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
429               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
430                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
431                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
432                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
433               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
434                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
435                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
436                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
437               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
438                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
439                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
440                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
441               zrxmax = MAXVAL( zr1(1:4) )
442               rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax )
443            END DO
444         END DO
445      END DO
446      CALL lbc_lnk( rx1, 'T', 1. )
447      !
448      zrxmax = MAXVAL( rx1 )
449      !
450      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
451      !
452      IF(lwp) THEN
453         WRITE(numout,*)
454         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
455         WRITE(numout,*) '~~~~~~~~~'
456      ENDIF
457      !
458   END SUBROUTINE dom_stiff
459
460
461   SUBROUTINE cfg_write
462      !!----------------------------------------------------------------------
463      !!                  ***  ROUTINE cfg_write  ***
464      !!                   
465      !! ** Purpose :   Create the "domain_cfg" file, a NetCDF file which
466      !!              contains all the ocean domain informations required to
467      !!              define an ocean configuration.
468      !!
469      !! ** Method  :   Write in a file all the arrays required to set up an
470      !!              ocean configuration.
471      !!
472      !! ** output file :   domain_cfg.nc : domain size, characteristics,
473      !horizontal mesh,
474      !!                              Coriolis parameter, depth and vertical
475      !scale factors
476      !!----------------------------------------------------------------------
477      INTEGER           ::   ji, jj, jk   ! dummy loop indices
478      INTEGER           ::   izco, izps, isco, icav
479      INTEGER           ::   inum     ! temprary units for 'domain_cfg.nc' file
480      CHARACTER(len=21) ::   clnam    ! filename (mesh and mask informations)
481      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! workspace
482      !!----------------------------------------------------------------------
483      !
484      IF(lwp) WRITE(numout,*)
485      IF(lwp) WRITE(numout,*) 'cfg_write : create the "domain_cfg.nc" file containing all required configuration information'
486      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
487      !
488      !                       ! ============================= !
489      !                       !  create 'domain_cfg.nc' file  !
490      !                       ! ============================= !
491      !         
492      clnam = 'domain_cfg'  ! filename (configuration information)
493      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
494     
495      !                             !==  global domain size  ==!
496      !
497      CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
498      CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
499      CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpk   , wp), ktype = jp_i4 )
500      !
501      !                             !==  domain characteristics  ==!
502      !
503      !                                   ! lateral boundary of the global
504      !                                   domain
505      CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
506      !
507      !                                   ! type of vertical coordinate
508      IF( ln_zco    ) THEN   ;   izco = 1   ;   ELSE   ;   izco = 0   ;   ENDIF
509      IF( ln_zps    ) THEN   ;   izps = 1   ;   ELSE   ;   izps = 0   ;   ENDIF
510      IF( ln_sco    ) THEN   ;   isco = 1   ;   ELSE   ;   isco = 0   ;   ENDIF
511      CALL iom_rstput( 0, 0, inum, 'ln_zco'   , REAL( izco, wp), ktype = jp_i4 )
512      CALL iom_rstput( 0, 0, inum, 'ln_zps'   , REAL( izps, wp), ktype = jp_i4 )
513      CALL iom_rstput( 0, 0, inum, 'ln_sco'   , REAL( isco, wp), ktype = jp_i4 )
514      !
515      !                                   ! ocean cavities under iceshelves
516      IF( ln_isfcav ) THEN   ;   icav = 1   ;   ELSE   ;   icav = 0   ;   ENDIF
517      CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
518      !
519      !                             !==  horizontal mesh  !
520      !
521      CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 )   ! latitude
522      CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
523      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
524      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
525      !                               
526      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )   ! longitude
527      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
528      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
529      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
530      !                               
531      CALL iom_rstput( 0, 0, inum, 'e1t'  , e1t  , ktype = jp_r8 )   ! i-scale factors (e1.)
532      CALL iom_rstput( 0, 0, inum, 'e1u'  , e1u  , ktype = jp_r8 )
533      CALL iom_rstput( 0, 0, inum, 'e1v'  , e1v  , ktype = jp_r8 )
534      CALL iom_rstput( 0, 0, inum, 'e1f'  , e1f  , ktype = jp_r8 )
535      !
536      CALL iom_rstput( 0, 0, inum, 'e2t'  , e2t  , ktype = jp_r8 )   ! j-scale factors (e2.)
537      CALL iom_rstput( 0, 0, inum, 'e2u'  , e2u  , ktype = jp_r8 )
538      CALL iom_rstput( 0, 0, inum, 'e2v'  , e2v  , ktype = jp_r8 )
539      CALL iom_rstput( 0, 0, inum, 'e2f'  , e2f  , ktype = jp_r8 )
540      !
541      CALL iom_rstput( 0, 0, inum, 'ff_f' , ff_f , ktype = jp_r8 )   ! coriolis factor
542      CALL iom_rstput( 0, 0, inum, 'ff_t' , ff_t , ktype = jp_r8 )
543      !
544      !                             !==  vertical mesh - 3D mask  ==!
545      !                                                     
546      CALL iom_rstput( 0, 0, inum, 'gdept_1d', gdept_1d, ktype = jp_r8 )   !  reference 1D-coordinate
547      CALL iom_rstput( 0, 0, inum, 'gdepw_1d', gdepw_1d, ktype = jp_r8 )
548      CALL iom_rstput( 0, 0, inum, 'e3t_1d'  , e3t_1d  , ktype = jp_r8 )
549      CALL iom_rstput( 0, 0, inum, 'e3w_1d'  , e3w_1d  , ktype = jp_r8 )
550      !                                                     
551      CALL iom_rstput( 0, 0, inum, 'gdept_0' , gdept_0 , ktype = jp_r8 )   !  depth (t- & w-points)
552      CALL iom_rstput( 0, 0, inum, 'gdepw_0' , gdepw_0 , ktype = jp_r8 )
553      !
554      CALL iom_rstput( 0, 0, inum, 'e3t_0'   , e3t_0   , ktype = jp_r8 )   !  vertical scale factors (e
555      CALL iom_rstput( 0, 0, inum, 'e3u_0'   , e3u_0   , ktype = jp_r8 )
556      CALL iom_rstput( 0, 0, inum, 'e3v_0'   , e3v_0   , ktype = jp_r8 )
557      CALL iom_rstput( 0, 0, inum, 'e3f_0'   , e3f_0   , ktype = jp_r8 )
558      CALL iom_rstput( 0, 0, inum, 'e3w_0'   , e3w_0   , ktype = jp_r8 )
559      CALL iom_rstput( 0, 0, inum, 'e3uw_0'  , e3uw_0  , ktype = jp_r8 )
560      CALL iom_rstput( 0, 0, inum, 'e3vw_0'  , e3vw_0  , ktype = jp_r8 )
561      !                                         
562      !                             !==  ocean top and bottom level  ==!
563      !
564      CALL iom_rstput( 0, 0, inum, 'bottom level' , REAL( mbkt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points
565      CALL iom_rstput( 0, 0, inum, 'top    level' , REAL( mikt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points (ISF)
566      !
567      IF( ln_sco ) THEN             ! s-coordinate: store grid stiffness ratio (Not required anyway)
568         CALL dom_stiff
569         !SF  CALL dom_stiff( z2d ) !commented because at compilation error:  The
570                                   !number of actual arguments cannot be greater than the number of dummy
571                                   !arguments.   [DOM_STIFF]
572                                   !CALL dom_stiff( z2d )
573                                   ! --------------^ compilation aborted for
574                                   !/workgpfs/rech/gzi/rgzi011/commit-simplif2-TOOLS/NEMOGCM/CONFIG/TEST/BLD/ppsrc/nemo/domain.f90
575         CALL iom_rstput( 0, 0, inum, 'stiffness', z2d )        !    ! Max. grid stiffness ratio
576      ENDIF
577      !
578      !                                ! ============================
579      !                                !        close the files
580      !                                ! ============================
581      CALL iom_close( inum )
582      !
583   END SUBROUTINE cfg_write
584
585
586
587   !!======================================================================
588END MODULE domain
Note: See TracBrowser for help on using the repository browser.