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/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 6973

Last change on this file since 6973 was 6596, checked in by gm, 8 years ago

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