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.
domrea.F90 in branches/2015/dev_CMCC_merge_2015/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/2015/dev_CMCC_merge_2015/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 6051

Last change on this file since 6051 was 6051, checked in by lovato, 8 years ago

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