source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 7522

Last change on this file since 7522 was 7522, checked in by cetlod, 4 years ago

3.6 stable : update the offline routines to be able to run passive tracers offline with linear free surface, see ticket #1775

  • Property svn:keywords set to Id
File size: 39.5 KB
Line 
1MODULE domrea
2   !!==============================================================================
3   !!                       ***  MODULE domrea   ***
4   !! Ocean initialization : domain initialization
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dom_init       : initialize the space and time domain
9   !!   dom_nam        : read and contral domain namelists
10   !!   dom_ctl        : control print for the ocean domain
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE oce             !
14   USE dom_oce         ! ocean space and time domain
15   USE phycst          ! physical constants
16   USE in_out_manager  ! I/O manager
17   USE lib_mpp         ! distributed memory computing library
18
19   USE iom
20   USE domstp          ! domain: set the time-step
21
22   USE lbclnk          ! lateral boundary condition - MPP exchanges
23   USE trc_oce         ! shared ocean/biogeochemical variables
24   USE wrk_nemo 
25   
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Routine accessibility
30   PUBLIC dom_rea       ! called by opa.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
37   !! $Id$
38   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43   SUBROUTINE dom_rea
44      !!----------------------------------------------------------------------
45      !!                  ***  ROUTINE dom_rea  ***
46      !!                   
47      !! ** Purpose :   Domain initialization. Call the routines that are
48      !!      required to create the arrays which define the space and time
49      !!      domain of the ocean model.
50      !!
51      !! ** Method  :
52      !!      - dom_stp: defined the model time step
53      !!      - dom_rea: read the meshmask file if nmsh=1
54      !!
55      !! History :
56      !!        !  90-10  (C. Levy - G. Madec)  Original code
57      !!        !  91-11  (G. Madec)
58      !!        !  92-01  (M. Imbard) insert time step initialization
59      !!        !  96-06  (G. Madec) generalized vertical coordinate
60      !!        !  97-02  (G. Madec) creation of domwri.F
61      !!        !  01-05  (E.Durand - G. Madec) insert closed sea
62      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
63      !!----------------------------------------------------------------------
64      !! * Local declarations
65      INTEGER ::   jk                ! dummy loop argument
66      INTEGER ::   iconf = 0         ! temporary integers
67      !!----------------------------------------------------------------------
68
69      IF(lwp) THEN
70         WRITE(numout,*)
71         WRITE(numout,*) 'dom_init : domain initialization'
72         WRITE(numout,*) '~~~~~~~~'
73      ENDIF
74
75      CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
76      CALL dom_msk      ! Masks
77      CALL dom_hgr      ! Horizontal grid
78      CALL dom_zgr      ! Vertical mesh and bathymetry option
79      !
80      e12t    (:,:) = e1t(:,:) * e2t(:,:)
81      e1e2t   (:,:) = e1t(:,:) * e2t(:,:)
82      e12u    (:,:) = e1u(:,:) * e2u(:,:)
83      e12v    (:,:) = e1v(:,:) * e2v(:,:)
84      e12f    (:,:) = e1f(:,:) * e2f(:,:)
85      r1_e12t (:,:) = 1._wp    / e12t(:,:)
86      r1_e12u (:,:) = 1._wp    / e12u(:,:)
87      r1_e12v (:,:) = 1._wp    / e12v(:,:)
88      r1_e12f (:,:) = 1._wp    / e12f(:,:)
89      re2u_e1u(:,:) = e2u(:,:) / e1u(:,:)
90      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
91      !
92      CALL dom_stp      ! Time step
93      CALL dom_ctl      ! Domain control
94
95   END SUBROUTINE dom_rea
96
97   SUBROUTINE dom_nam
98      !!----------------------------------------------------------------------
99      !!                     ***  ROUTINE dom_nam  ***
100      !!                   
101      !! ** Purpose :   read domaine namelists and print the variables.
102      !!
103      !! ** input   : - namrun namelist
104      !!              - namdom namelist
105      !!              - namcla namelist
106      !!----------------------------------------------------------------------
107      USE ioipsl
108      INTEGER  ::   ios                 ! Local integer output status for namelist read
109      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
110         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
111         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
112         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler
113      NAMELIST/namdom/ nn_bathy , rn_bathy, rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin,   &
114         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,            &
115         &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea , ln_crs, &
116         &             jphgr_msh, &
117         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
118         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
119         &             ppa2, ppkth2, ppacr2
120      NAMELIST/namcla/ nn_cla
121#if defined key_netcdf4
122      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
123#endif
124      !!----------------------------------------------------------------------
125
126      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
127      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
128901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
129
130      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
131      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
132902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
133      IF(lwm) WRITE ( numond, namrun )
134      !
135      IF(lwp) THEN                  ! control print
136         WRITE(numout,*)
137         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
138         WRITE(numout,*) '~~~~~~~ '
139         WRITE(numout,*) '   Namelist namrun' 
140         WRITE(numout,*) '      job number                      nn_no      = ', nn_no
141         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp
142         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart
143         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl
144         WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000
145         WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend
146         WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0
147         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy
148         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate
149         WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock
150         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write
151         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn
152         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland
153         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta
154         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
155         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
156      ENDIF
157      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
158      cexper = cn_exp
159      nrstdt = nn_rstctl
160      nit000 = nn_it000
161      nitend = nn_itend
162      ndate0 = nn_date0
163      nleapy = nn_leapy
164      ninist = nn_istate
165      nstock = nn_stock
166      nstocklist = nn_stocklist
167      nwrite = nn_write
168      !                             ! control of output frequency
169      IF ( nstock == 0 .OR. nstock > nitend ) THEN
170         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
171         CALL ctl_warn( ctmp1 )
172         nstock = nitend
173      ENDIF
174      IF ( nwrite == 0 ) THEN
175         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
176         CALL ctl_warn( ctmp1 )
177         nwrite = nitend
178      ENDIF
179
180      ! parameters correspondting to nit000 - 1 (as we start the step loop with a call to day)
181      ndastp = ndate0 - 1        ! ndate0 read in the namelist in dom_nam, we assume that we start run at 00:00
182      adatrj = ( REAL( nit000-1, wp ) * rdttra(1) ) / rday
183
184#if defined key_agrif
185      IF( Agrif_Root() ) THEN
186#endif
187      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
188      CASE (  1 ) 
189         CALL ioconf_calendar('gregorian')
190         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
191      CASE (  0 )
192         CALL ioconf_calendar('noleap')
193         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
194      CASE ( 30 )
195         CALL ioconf_calendar('360d')
196         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
197      END SELECT
198#if defined key_agrif
199      ENDIF
200#endif
201
202      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
203      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
204903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
205
206      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
207      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
208904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
209      IF(lwm) WRITE ( numond, namdom )
210
211
212      IF(lwp) THEN
213         WRITE(numout,*) 
214         WRITE(numout,*) '   Namelist namdom : space & time domain'
215         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy
216         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy
217         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
218         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
219         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
220         WRITE(numout,*) '      create mesh/mask file(s)          nn_msh       = ', nn_msh
221         WRITE(numout,*) '           = 0   no file created                 '
222         WRITE(numout,*) '           = 1   mesh_mask                       '
223         WRITE(numout,*) '           = 2   mesh and mask                   '
224         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask      '
225         WRITE(numout,*) '      ocean time step                      rn_rdt    = ', rn_rdt
226         WRITE(numout,*) '      asselin time filter parameter        rn_atfp   = ', rn_atfp
227         WRITE(numout,*) '      time-splitting: nb of sub time-step  nn_baro   = ', nn_baro
228         WRITE(numout,*) '      acceleration of converge             nn_acc    = ', nn_acc
229         WRITE(numout,*) '        nn_acc=1: surface tracer rdt       rn_rdtmin = ', rn_rdtmin
230         WRITE(numout,*) '                  bottom  tracer rdt       rdtmax    = ', rn_rdtmax
231         WRITE(numout,*) '                  depth of transition      rn_rdth   = ', rn_rdth
232         WRITE(numout,*) '      suppression of closed seas (=0)      nn_closea = ', nn_closea
233         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh
234         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0
235         WRITE(numout,*) '      latitude  of first raw and column T-point ppgphi0 = ', ppgphi0
236         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_deg        = ', ppe1_deg
237         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_deg        = ', ppe2_deg
238         WRITE(numout,*) '      zonal      grid-spacing (degrees) ppe1_m          = ', ppe1_m
239         WRITE(numout,*) '      meridional grid-spacing (degrees) ppe2_m          = ', ppe2_m
240         WRITE(numout,*) '      ORCA r4, r2 and r05 coefficients  ppsur           = ', ppsur
241         WRITE(numout,*) '                                        ppa0            = ', ppa0
242         WRITE(numout,*) '                                        ppa1            = ', ppa1
243         WRITE(numout,*) '                                        ppkth           = ', ppkth
244         WRITE(numout,*) '                                        ppacr           = ', ppacr
245         WRITE(numout,*) '      Minimum vertical spacing ppdzmin                  = ', ppdzmin
246         WRITE(numout,*) '      Maximum depth pphmax                              = ', pphmax
247         WRITE(numout,*) '      Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh
248         WRITE(numout,*) '      Double tanh function parameters ppa2              = ', ppa2
249         WRITE(numout,*) '                                      ppkth2            = ', ppkth2
250         WRITE(numout,*) '                                      ppacr2            = ', ppacr2
251      ENDIF
252
253      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon)
254      e3zps_min = rn_e3zps_min
255      e3zps_rat = rn_e3zps_rat
256      nmsh      = nn_msh
257      nacc      = nn_acc
258      atfp      = rn_atfp
259      rdt       = rn_rdt
260      rdtmin    = rn_rdtmin
261      rdtmax    = rn_rdtmin
262      rdth      = rn_rdth
263
264      REWIND( numnam_ref )              ! Namelist namcla in reference namelist : Cross land advection
265      READ  ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)
266905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )
267
268      REWIND( numnam_cfg )              ! Namelist namcla in configuration namelist : Cross land advection
269      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )
270906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )
271      IF(lwm) WRITE( numond, namcla )
272
273      IF(lwp) THEN
274         WRITE(numout,*)
275         WRITE(numout,*) '   Namelist namcla'
276         WRITE(numout,*) '      cross land advection                 nn_cla    = ', nn_cla
277      ENDIF
278
279#if defined key_netcdf4
280      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
281      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
282      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
283907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
284
285      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
286      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
287908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
288      IF(lwm) WRITE( numond, namnc4 )
289      IF(lwp) THEN                        ! control print
290         WRITE(numout,*)
291         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
292         WRITE(numout,*) '      number of chunks in i-dimension      nn_nchunks_i   = ', nn_nchunks_i
293         WRITE(numout,*) '      number of chunks in j-dimension      nn_nchunks_j   = ', nn_nchunks_j
294         WRITE(numout,*) '      number of chunks in k-dimension      nn_nchunks_k   = ', nn_nchunks_k
295         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression ln_nc4zip = ', ln_nc4zip
296      ENDIF
297
298      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
299      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
300      snc4set%ni   = nn_nchunks_i
301      snc4set%nj   = nn_nchunks_j
302      snc4set%nk   = nn_nchunks_k
303      snc4set%luse = ln_nc4zip
304#else
305      snc4set%luse = .FALSE.        ! No NetCDF 4 case
306#endif
307      !
308   END SUBROUTINE dom_nam
309
310   SUBROUTINE dom_msk
311      !!---------------------------------------------------------------------
312      !!                 ***  ROUTINE dom_msk  ***
313      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
314      !!      ocean mask informations and defines the interior domain T-mask.
315      !!
316      !! ** Method  :  Read in a file all the arrays generated in routines
317      !!               dommsk:   'mask.nc' file
318      !!              The interior ocean/land mask is computed from tmask
319      !!              setting to zero the duplicated row and lines due to
320      !!              MPP exchange halos, est-west cyclic and north fold
321      !!              boundary conditions.
322      !!
323      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point
324      !!               tpol     : ???
325      !!----------------------------------------------------------------------
326      !
327      INTEGER  ::  inum   ! local integers
328      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
329      INTEGER  ::   iif, iil, ijf, ijl       ! local integers
330      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
331      !
332      !!---------------------------------------------------------------------
333     
334
335
336      IF(lwp) WRITE(numout,*)
337      IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)'
338      IF(lwp) WRITE(numout,*) '~~~~~~~'
339
340      CALL wrk_alloc( jpi, jpj, zmbk )
341      zmbk(:,:) = 0._wp
342
343      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
344      CALL iom_open( 'mask', inum )
345
346         !                                                         ! masks (inum2)
347      CALL iom_get( inum, jpdom_data, 'tmask', tmask )
348      CALL iom_get( inum, jpdom_data, 'umask', umask )
349      CALL iom_get( inum, jpdom_data, 'vmask', vmask )
350      CALL iom_get( inum, jpdom_data, 'fmask', fmask )
351
352      CALL lbc_lnk( tmask, 'T', 1._wp )    ! Lateral boundary conditions
353      CALL lbc_lnk( umask, 'U', 1._wp )     
354      CALL lbc_lnk( vmask, 'V', 1._wp )
355      CALL lbc_lnk( fmask, 'F', 1._wp )
356
357#if defined key_c1d
358      ! set umask and vmask equal tmask in 1D configuration
359      IF(lwp) WRITE(numout,*)
360      IF(lwp) WRITE(numout,*) '**********  1D configuration : set umask and vmask equal tmask ********'
361      IF(lwp) WRITE(numout,*) '**********                                                     ********'
362
363      umask(:,:,:) = tmask(:,:,:)
364      vmask(:,:,:) = tmask(:,:,:)
365#endif
366
367#if defined key_degrad
368      CALL iom_get( inum, jpdom_data, 'facvolt', facvol )
369#endif
370
371      CALL iom_get( inum, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points
372      mbathy (:,:) = INT( zmbk(:,:) )
373      misfdep(:,:) = 1                                               ! ice shelf case not yet done
374     
375      CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points
376
377      !                                     ! ============================
378      !                                     !        close the files
379      !                                     ! ============================
380
381      !
382      ! Interior domain mask (used for global sum)
383      ! --------------------
384      ssmask(:,:)  = tmask(:,:,1)
385      tmask_i(:,:) = tmask(:,:,1)
386      iif = jpreci                        ! thickness of exchange halos in i-axis
387      iil = nlci - jpreci + 1
388      ijf = jprecj                        ! thickness of exchange halos in j-axis
389      ijl = nlcj - jprecj + 1
390      !
391      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns
392      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns)
393      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows
394      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows)
395      !
396      !                                   ! north fold mask
397      tpol(1:jpiglo) = 1._wp
398      !                               
399      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot
400      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot
401      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row
402         IF( mjg(ijl-1) == jpjglo-1 ) THEN
403            DO ji = iif+1, iil-1
404               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji))
405            END DO
406         ENDIF
407      ENDIF 
408      !
409      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at
410      ! least 1 wet u point
411      DO jj = 1, jpjm1
412         DO ji = 1, fs_jpim1   ! vector loop
413            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
414            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
415         END DO
416         DO ji = 1, jpim1      ! NO vector opt.
417            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
418               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
419         END DO
420      END DO
421      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
422      CALL lbc_lnk( vmask_i, 'V', 1._wp )
423      CALL lbc_lnk( fmask_i, 'F', 1._wp )
424
425      ! 3. Ocean/land mask at wu-, wv- and w points
426      !----------------------------------------------
427      wmask (:,:,1) = tmask(:,:,1) ! ????????
428      wumask(:,:,1) = umask(:,:,1) ! ????????
429      wvmask(:,:,1) = vmask(:,:,1) ! ????????
430      DO jk = 2, jpk
431         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
432         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
433         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
434      END DO
435      !
436      CALL wrk_dealloc( jpi, jpj, zmbk )
437      !
438      CALL iom_close( inum )
439      !
440   END SUBROUTINE dom_msk
441
442   SUBROUTINE zgr_bot_level
443      !!----------------------------------------------------------------------
444      !!                    ***  ROUTINE zgr_bot_level  ***
445      !!
446      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
447      !!
448      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
449      !!
450      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
451      !!                                     ocean level at t-, u- & v-points
452      !!                                     (min value = 1 over land)
453      !!----------------------------------------------------------------------
454      !
455      INTEGER ::   ji, jj   ! dummy loop indices
456      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
457      !!----------------------------------------------------------------------
458
459      !
460      IF(lwp) WRITE(numout,*)
461      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
462      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
463      !
464      CALL wrk_alloc( jpi, jpj, zmbk )
465      !
466      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
467      mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf)
468      !                                     ! bottom k-index of W-level = mbkt+1
469      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
470         DO ji = 1, jpim1
471            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
472            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
473         END DO
474      END DO
475      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
476      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
477      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
478      !
479      CALL wrk_dealloc( jpi, jpj, zmbk )
480      !
481   END SUBROUTINE zgr_bot_level
482
483   SUBROUTINE dom_hgr
484      !!----------------------------------------------------------------------
485      !!                  ***  ROUTINE dom_hgr  ***
486      !!                   
487      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
488      !!      ocean horizontal mesh informations
489      !!
490      !! ** Method  :   Read in a file all the arrays generated in routines
491      !!                domhgr:   'mesh_hgr.nc' file
492      !!----------------------------------------------------------------------
493      !!
494      INTEGER ::   ji, jj   ! dummy loop indices
495      INTEGER  ::  inum    ! local integers
496      !!----------------------------------------------------------------------
497
498      IF(lwp) WRITE(numout,*)
499      IF(lwp) WRITE(numout,*) 'dom_grd_hgr : read NetCDF mesh and mask information file(s)'
500      IF(lwp) WRITE(numout,*) '~~~~~~~'
501
502      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
503      CALL iom_open( 'mesh_hgr', inum )
504
505      !                                                         ! horizontal mesh (inum3)
506      CALL iom_get( inum, jpdom_data, 'glamt', glamt )
507      CALL iom_get( inum, jpdom_data, 'glamu', glamu )
508      CALL iom_get( inum, jpdom_data, 'glamv', glamv )
509      CALL iom_get( inum, jpdom_data, 'glamf', glamf )
510
511      CALL iom_get( inum, jpdom_data, 'gphit', gphit )
512      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu )
513      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv )
514      CALL iom_get( inum, jpdom_data, 'gphif', gphif )
515
516      CALL iom_get( inum, jpdom_data, 'e1t', e1t )
517      CALL iom_get( inum, jpdom_data, 'e1u', e1u )
518      CALL iom_get( inum, jpdom_data, 'e1v', e1v )
519     
520      CALL iom_get( inum, jpdom_data, 'e2t', e2t )
521      CALL iom_get( inum, jpdom_data, 'e2u', e2u )
522      CALL iom_get( inum, jpdom_data, 'e2v', e2v )
523
524      CALL iom_get( inum, jpdom_data, 'ff', ff )
525
526
527      ! Control printing : Grid informations (if not restart)
528      ! ----------------
529
530      IF(lwp .AND. .NOT.ln_rstart ) THEN
531         WRITE(numout,*)
532         WRITE(numout,*) '          longitude and e1 scale factors'
533         WRITE(numout,*) '          ------------------------------'
534         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
535            glamv(ji,1), glamf(ji,1),   &
536            e1t(ji,1), e1u(ji,1),   &
537            e1v(ji,1), ji = 1, jpi,10)
538
539         WRITE(numout,*)
540         WRITE(numout,*) '          latitude and e2 scale factors'
541         WRITE(numout,*) '          -----------------------------'
542         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
543            &                     gphiv(1,jj), gphif(1,jj),   &
544            &                     e2t  (1,jj), e2u  (1,jj),   &
545            &                     e2v  (1,jj), jj = 1, jpj, 10 )
546      ENDIF
547
548      !                                     ! ============================
549      !                                     !        close the files
550      !                                     ! ============================
551      CALL iom_close( inum )
552      !
5539300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
554            f19.10, 1x, f19.10, 1x, f19.10 )
555   END SUBROUTINE dom_hgr
556
557
558   SUBROUTINE dom_zgr
559      !!----------------------------------------------------------------------
560      !!                ***  ROUTINE dom_zgr  ***
561      !!                   
562      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
563      !!      ocean horizontal mesh informations and/or set the depth of model levels
564      !!      and the resulting vertical scale factors.
565      !!
566      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
567      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
568      !!              - vertical coordinate (gdep., e3.) depending on the
569      !!                coordinate chosen :
570      !!                   ln_zco=T   z-coordinate 
571      !!                   ln_zps=T   z-coordinate with partial steps
572      !!                   ln_zco=T   s-coordinate
573      !!
574      !! ** Action  :   define gdep., e3., mbathy and bathy
575      !!----------------------------------------------------------------------
576      INTEGER  ::  ioptio = 0   ! temporary integer
577      INTEGER  ::  inum, ios
578      INTEGER  ::  ji, jj, jk, ik
579      REAL(wp) ::  zrefdep
580      !!
581      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
582      REAL(wp), POINTER, DIMENSION(:,:) :: zprt, zprw
583      !!----------------------------------------------------------------------
584
585      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
586      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
587901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
588
589      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
590      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
591902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
592      IF(lwm) WRITE ( numond, namzgr )
593
594      IF(lwp) THEN                     ! Control print
595         WRITE(numout,*)
596         WRITE(numout,*) 'dom_zgr : vertical coordinate'
597         WRITE(numout,*) '~~~~~~~'
598         WRITE(numout,*) '          Namelist namzgr : set vertical coordinate'
599         WRITE(numout,*) '             z-coordinate - full steps      ln_zco    = ', ln_zco
600         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps    = ', ln_zps
601         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
602         WRITE(numout,*) '             ice shelf cavity               ln_isfcav = ', ln_isfcav
603      ENDIF
604
605      ioptio = 0                       ! Check Vertical coordinate options
606      IF( ln_zco ) ioptio = ioptio + 1
607      IF( ln_zps ) ioptio = ioptio + 1
608      IF( ln_sco ) ioptio = ioptio + 1
609      IF( ln_isfcav ) ioptio = 33
610      IF ( ioptio /= 1  )   CALL ctl_stop( ' none or several vertical coordinate options used' )
611      IF ( ioptio == 33 )   CALL ctl_stop( ' isf cavity with off line module not yet done    ' )
612
613      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
614      CALL iom_open( 'mesh_zgr', inum )
615
616      CALL iom_get( inum, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
617      CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', gdepw_1d )
618      IF( ln_zco .OR. ln_zps ) THEN
619         CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors
620         CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
621      ENDIF
622
623!!gm BUG in s-coordinate this does not work!
624      ! deepest/shallowest W level Above/Below ~10m
625      zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness)
626      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
627      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
628!!gm end bug
629
630      IF(lwp) THEN
631         WRITE(numout,*)
632         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
633         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
634         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
635      ENDIF
636
637      DO jk = 1, jpk
638         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' )
639         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' )
640      END DO
641
642      IF( lk_vvl ) THEN
643          CALL iom_get( inum, jpdom_data, 'e3t_0', e3t_0(:,:,:) )
644          CALL iom_get( inum, jpdom_data, 'e3u_0', e3u_0(:,:,:) )
645          CALL iom_get( inum, jpdom_data, 'e3v_0', e3v_0(:,:,:) )
646          CALL iom_get( inum, jpdom_data, 'e3w_0', e3w_0(:,:,:) )
647          CALL iom_get( inum, jpdom_data, 'gdept_0', gdept_0(:,:,:) )
648          CALL iom_get( inum, jpdom_data, 'gdepw_0', gdepw_0(:,:,:) )
649          ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at  T-points
650          DO jk = 1, jpk
651             ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
652          END DO
653      ELSE
654         IF( ln_sco ) THEN                                         ! s-coordinate
655            CALL iom_get( inum, jpdom_data, 'hbatt', hbatt )
656            CALL iom_get( inum, jpdom_data, 'hbatu', hbatu )
657            CALL iom_get( inum, jpdom_data, 'hbatv', hbatv )
658            CALL iom_get( inum, jpdom_data, 'hbatf', hbatf )
659           
660            CALL iom_get( inum, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef.
661            CALL iom_get( inum, jpdom_unknown, 'gsigw', gsigw )
662            CALL iom_get( inum, jpdom_unknown, 'gsi3w', gsi3w ) 
663            CALL iom_get( inum, jpdom_unknown, 'esigt', esigt )
664            CALL iom_get( inum, jpdom_unknown, 'esigw', esigw )
665
666            CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors
667            CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
668            CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
669            CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
670         ENDIF
671
672 
673         IF( ln_zps ) THEN                                           ! z-coordinate - partial steps
674            !
675            IF( iom_varid( inum, 'e3t_0', ldstop = .FALSE. ) > 0 ) THEN
676               CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) )
677               CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
678               CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
679               CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
680            ELSE                                                        ! 2D bottom scale factors
681               CALL iom_get( inum, jpdom_data, 'e3t_ps', e3tp )
682               CALL iom_get( inum, jpdom_data, 'e3w_ps', e3wp )
683               !                                                        ! deduces the 3D scale factors
684               DO jk = 1, jpk
685                  fse3t_n(:,:,jk) = e3t_1d(jk)                                    ! set to the ref. factors
686                  fse3u_n(:,:,jk) = e3t_1d(jk)
687                  fse3v_n(:,:,jk) = e3t_1d(jk)
688                  fse3w_n(:,:,jk) = e3w_1d(jk)
689               END DO
690               DO jj = 1,jpj                                                  ! adjust the deepest values
691                  DO ji = 1,jpi
692                     ik = mbkt(ji,jj)
693                     fse3t_n(ji,jj,ik) = e3tp(ji,jj) * tmask(ji,jj,1) + e3t_1d(1) * ( 1._wp - tmask(ji,jj,1) )
694                     fse3w_n(ji,jj,ik) = e3wp(ji,jj) * tmask(ji,jj,1) + e3w_1d(1) * ( 1._wp - tmask(ji,jj,1) )
695                  END DO
696               END DO
697               DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
698                  DO jj = 1, jpjm1
699                     DO ji = 1, jpim1
700                        fse3u_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji+1,jj,jk) )
701                        fse3v_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji,jj+1,jk) )
702                     END DO
703                  END DO
704               END DO
705               CALL lbc_lnk( fse3u_n(:,:,:) , 'U', 1._wp )   ;   CALL lbc_lnk( fse3uw_n(:,:,:), 'U', 1._wp )   ! lateral boundary conditions
706               CALL lbc_lnk( fse3v_n(:,:,:) , 'V', 1._wp )   ;   CALL lbc_lnk( fse3vw_n(:,:,:), 'V', 1._wp )
707               !
708               DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
709                  WHERE( fse3u_n(:,:,jk) == 0._wp )   fse3u_n(:,:,jk) = e3t_1d(jk)
710                  WHERE( fse3v_n(:,:,jk) == 0._wp )   fse3v_n(:,:,jk) = e3t_1d(jk)
711               END DO
712            END IF
713
714            IF( iom_varid( inum, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level
715               CALL iom_get( inum, jpdom_data, 'gdept_0', fsdept_n(:,:,:) )
716               CALL iom_get( inum, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) )
717            ELSE                                                           ! 2D bottom depth
718               CALL wrk_alloc( jpi, jpj, zprt, zprw )
719               !
720               CALL iom_get( inum, jpdom_data, 'hdept', zprt )
721               CALL iom_get( inum, jpdom_data, 'hdepw', zprw )
722               !
723               DO jk = 1, jpk                                              ! deduces the 3D depth
724                  fsdept_n(:,:,jk) = gdept_1d(jk)
725                  fsdepw_n(:,:,jk) = gdepw_1d(jk)
726               END DO
727               DO jj = 1, jpj
728                  DO ji = 1, jpi
729                     ik = mbkt(ji,jj)
730                     IF( ik > 0 ) THEN
731                        fsdepw_n(ji,jj,ik+1) = zprw(ji,jj)
732                        fsdept_n(ji,jj,ik  ) = zprt(ji,jj)
733                        fsdept_n(ji,jj,ik+1) = fsdept_n(ji,jj,ik) + fse3t_n(ji,jj,ik)
734                     ENDIF
735                  END DO
736               END DO
737               CALL wrk_dealloc( jpi, jpj, zprt, zprw )
738            ENDIF
739            !
740         ENDIF
741
742         IF( ln_zco ) THEN           ! Vertical coordinates and scales factors
743            DO jk = 1, jpk
744               fse3t_n(:,:,jk) = e3t_1d(jk)                              ! set to the ref. factors
745               fse3u_n(:,:,jk) = e3t_1d(jk)
746               fse3v_n(:,:,jk) = e3t_1d(jk)
747               fse3w_n(:,:,jk) = e3w_1d(jk)
748               fsdept_n(:,:,jk) = gdept_1d(jk)
749               fsdepw_n(:,:,jk) = gdepw_1d(jk)
750            END DO
751         ENDIF
752         !
753      ENDIF
754      !                                     ! ============================
755      !                                     !        close the files
756      !                                     ! ============================
757      CALL iom_close( inum )
758      !
759      !
760   END SUBROUTINE dom_zgr
761
762   SUBROUTINE dom_ctl
763      !!----------------------------------------------------------------------
764      !!                     ***  ROUTINE dom_ctl  ***
765      !!
766      !! ** Purpose :   Domain control.
767      !!
768      !! ** Method  :   compute and print extrema of masked scale factors
769      !!
770      !! History :
771      !!   8.5  !  02-08  (G. Madec)    Original code
772      !!----------------------------------------------------------------------
773      !! * Local declarations
774      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
775      INTEGER, DIMENSION(2) ::   iloc      !
776      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
777      !!----------------------------------------------------------------------
778
779      ! Extrema of the scale factors
780
781      IF(lwp)WRITE(numout,*)
782      IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
783      IF(lwp)WRITE(numout,*) '~~~~~~~'
784
785      IF (lk_mpp) THEN
786         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
787         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
788         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
789         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
790      ELSE
791         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )   
792         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )   
793         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )   
794         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )   
795
796         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
797         iimi1 = iloc(1) + nimpp - 1
798         ijmi1 = iloc(2) + njmpp - 1
799         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
800         iimi2 = iloc(1) + nimpp - 1
801         ijmi2 = iloc(2) + njmpp - 1
802         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
803         iima1 = iloc(1) + nimpp - 1
804         ijma1 = iloc(2) + njmpp - 1
805         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
806         iima2 = iloc(1) + nimpp - 1
807         ijma2 = iloc(2) + njmpp - 1
808      ENDIF
809
810      IF(lwp) THEN
811         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
812         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
813         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
814         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
815      ENDIF
816
817   END SUBROUTINE dom_ctl
818
819   !!======================================================================
820END MODULE domrea
821
Note: See TracBrowser for help on using the repository browser.