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

source: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 6069

Last change on this file since 6069 was 6069, checked in by timgraham, 8 years ago

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