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

Last change on this file since 8617 was 8617, checked in by cetlod, 3 years ago

bugfix in offline : remove undefined & unused variable to avoid runtime error in debug mode

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