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 @ 6095

Last change on this file since 6095 was 6095, checked in by cetlod, 9 years ago

Add time varying part of vertical coordinate system in offline
  • Property svn:keywords set to Id
File size: 42.7 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, ln_linssh
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         WRITE(numout,*) '             Linear free surface            ln_linssh = ', ln_linssh
343      ENDIF
344
345      ioptio = 0                       ! Check Vertical coordinate options
346      IF( ln_zco ) ioptio = ioptio + 1
347      IF( ln_zps ) ioptio = ioptio + 1
348      IF( ln_sco ) ioptio = ioptio + 1
349      IF( ln_isfcav ) ioptio = 33
350      IF ( ioptio /= 1  )   CALL ctl_stop( ' none or several vertical coordinate options used' )
351      IF ( ioptio == 33 )   CALL ctl_stop( ' isf cavity with off line module not yet done    ' )
352
353   END SUBROUTINE dom_zgr
354
355
356   SUBROUTINE dom_ctl
357      !!----------------------------------------------------------------------
358      !!                     ***  ROUTINE dom_ctl  ***
359      !!
360      !! ** Purpose :   Domain control.
361      !!
362      !! ** Method  :   compute and print extrema of masked scale factors
363      !!
364      !!----------------------------------------------------------------------
365      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
366      INTEGER, DIMENSION(2) ::   iloc      !
367      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
368      !!----------------------------------------------------------------------
369
370      ! Extrema of the scale factors
371
372      IF(lwp)WRITE(numout,*)
373      IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
374      IF(lwp)WRITE(numout,*) '~~~~~~~'
375
376      IF (lk_mpp) THEN
377         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
378         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
379         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
380         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
381      ELSE
382         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )   
383         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )   
384         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )   
385         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )   
386
387         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
388         iimi1 = iloc(1) + nimpp - 1
389         ijmi1 = iloc(2) + njmpp - 1
390         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
391         iimi2 = iloc(1) + nimpp - 1
392         ijmi2 = iloc(2) + njmpp - 1
393         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
394         iima1 = iloc(1) + nimpp - 1
395         ijma1 = iloc(2) + njmpp - 1
396         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
397         iima2 = iloc(1) + nimpp - 1
398         ijma2 = iloc(2) + njmpp - 1
399      ENDIF
400      !
401      IF(lwp) THEN
402         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
403         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
404         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
405         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
406      ENDIF
407      !
408   END SUBROUTINE dom_ctl
409
410
411   SUBROUTINE dom_grd
412      !!----------------------------------------------------------------------
413      !!                  ***  ROUTINE dom_grd  ***
414      !!                   
415      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
416      !!      ocean domain informations (mesh and mask arrays). This (these)
417      !!      file(s) is (are) used for visualisation (SAXO software) and
418      !!      diagnostic computation.
419      !!
420      !! ** Method  :   Read in a file all the arrays generated in routines
421      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
422      !!      the vertical coord. used (z-coord, partial steps, s-coord)
423      !!                    nmsh = 1  :   'mesh_mask.nc' file
424      !!                         = 2  :   'mesh.nc' and mask.nc' files
425      !!                         = 3  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
426      !!                                  'mask.nc' files
427      !!      For huge size domain, use option 2 or 3 depending on your
428      !!      vertical coordinate.
429      !!
430      !! ** input file :
431      !!      meshmask.nc  : domain size, horizontal grid-point position,
432      !!                     masks, depth and vertical scale factors
433      !!----------------------------------------------------------------------
434      USE iom
435      !!
436      INTEGER  ::   ji, jj, jk   ! dummy loop indices
437      INTEGER  ::   ik, inum0 , inum1 , inum2 , inum3 , inum4   ! local integers
438      REAL(wp) ::   zrefdep         ! local real
439      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk, zprt, zprw
440      !!----------------------------------------------------------------------
441
442      IF(lwp) WRITE(numout,*)
443      IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)'
444      IF(lwp) WRITE(numout,*) '~~~~~~~'
445
446      CALL wrk_alloc( jpi, jpj, zmbk, zprt, zprw )
447
448      zmbk(:,:) = 0._wp
449
450      SELECT CASE (nmsh)
451         !                                     ! ============================
452         CASE ( 1 )                            !  create 'mesh_mask.nc' file
453            !                                  ! ============================
454
455            IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
456            CALL iom_open( 'mesh_mask', inum0 )
457
458            inum2 = inum0                                            ! put all the informations
459            inum3 = inum0                                            ! in unit inum0
460            inum4 = inum0
461
462            !                                  ! ============================
463         CASE ( 2 )                            !  create 'mesh.nc' and
464            !                                  !         'mask.nc' files
465            !                                  ! ============================
466
467            IF(lwp) WRITE(numout,*) '          two files in "mesh.nc" and "mask.nc" '
468            CALL iom_open( 'mesh', inum1 )
469            CALL iom_open( 'mask', inum2 )
470
471            inum3 = inum1                                            ! put mesh informations
472            inum4 = inum1                                            ! in unit inum1
473
474            !                                  ! ============================
475         CASE ( 3 )                            !  create 'mesh_hgr.nc'
476            !                                  !         'mesh_zgr.nc' and
477            !                                  !         'mask.nc'     files
478            !                                  ! ============================
479
480            IF(lwp) WRITE(numout,*) '          three files in "mesh_hgr.nc" , "mesh_zgr.nc" and "mask.nc" '
481            CALL iom_open( 'mesh_hgr', inum3 ) ! create 'mesh_hgr.nc'
482            CALL iom_open( 'mesh_zgr', inum4 ) ! create 'mesh_zgr.nc'
483            CALL iom_open( 'mask'    , inum2 ) ! create 'mask.nc'
484
485            !                                  ! ===========================
486         CASE DEFAULT                          ! return error
487            !                                  ! mesh has to be provided
488            !                                  ! ===========================
489            CALL ctl_stop( ' OFFLINE mode requires the input mesh mask(s). ',   &
490            &                                 ' Invalid nn_msh value in the namelist (0 is not allowed)' )
491
492         END SELECT
493
494         !                                                         ! masks (inum2)
495         CALL iom_get( inum2, jpdom_data, 'tmask', tmask )
496         CALL iom_get( inum2, jpdom_data, 'umask', umask )
497         CALL iom_get( inum2, jpdom_data, 'vmask', vmask )
498         CALL iom_get( inum2, jpdom_data, 'fmask', fmask )
499
500         CALL lbc_lnk( tmask, 'T', 1._wp )    ! Lateral boundary conditions
501         CALL lbc_lnk( umask, 'U', 1._wp )     
502         CALL lbc_lnk( vmask, 'V', 1._wp )
503         CALL lbc_lnk( fmask, 'F', 1._wp )
504
505#if defined key_c1d
506         ! set umask and vmask equal tmask in 1D configuration
507         IF(lwp) WRITE(numout,*)
508         IF(lwp) WRITE(numout,*) '**********  1D configuration : set umask and vmask equal tmask ********'
509         IF(lwp) WRITE(numout,*) '**********                                                     ********'
510
511         umask(:,:,:) = tmask(:,:,:)
512         vmask(:,:,:) = tmask(:,:,:)
513#endif
514
515#if defined key_degrad
516         CALL iom_get( inum2, jpdom_data, 'facvolt', facvol )
517#endif
518         !                                                         ! horizontal mesh (inum3)
519         CALL iom_get( inum3, jpdom_data, 'glamt', glamt )
520         CALL iom_get( inum3, jpdom_data, 'glamu', glamu )
521         CALL iom_get( inum3, jpdom_data, 'glamv', glamv )
522         CALL iom_get( inum3, jpdom_data, 'glamf', glamf )
523
524         CALL iom_get( inum3, jpdom_data, 'gphit', gphit )
525         CALL iom_get( inum3, jpdom_data, 'gphiu', gphiu )
526         CALL iom_get( inum3, jpdom_data, 'gphiv', gphiv )
527         CALL iom_get( inum3, jpdom_data, 'gphif', gphif )
528
529         CALL iom_get( inum3, jpdom_data, 'e1t', e1t )
530         CALL iom_get( inum3, jpdom_data, 'e1u', e1u )
531         CALL iom_get( inum3, jpdom_data, 'e1v', e1v )
532         
533         CALL iom_get( inum3, jpdom_data, 'e2t', e2t )
534         CALL iom_get( inum3, jpdom_data, 'e2u', e2u )
535         CALL iom_get( inum3, jpdom_data, 'e2v', e2v )
536
537         CALL iom_get( inum3, jpdom_data, 'ff', ff )
538
539         CALL iom_get( inum4, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points
540         mbathy (:,:) = INT( zmbk(:,:) )
541         misfdep(:,:) = 1                                               ! ice shelf case not yet done
542         
543         CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points
544         !
545         IF( ln_sco ) THEN                                         ! s-coordinate
546            CALL iom_get( inum4, jpdom_data, 'hbatt', hbatt )
547            CALL iom_get( inum4, jpdom_data, 'hbatu', hbatu )
548            CALL iom_get( inum4, jpdom_data, 'hbatv', hbatv )
549            CALL iom_get( inum4, jpdom_data, 'hbatf', hbatf )
550           
551            CALL iom_get( inum4, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef.
552            CALL iom_get( inum4, jpdom_unknown, 'gsigw', gsigw )
553            CALL iom_get( inum4, jpdom_unknown, 'gsi3w', gsi3w ) 
554            CALL iom_get( inum4, jpdom_unknown, 'esigt', esigt )
555            CALL iom_get( inum4, jpdom_unknown, 'esigw', esigw )
556
557            CALL iom_get( inum4, jpdom_data, 'e3t_0', e3t_0(:,:,:) ) ! scale factors
558            CALL iom_get( inum4, jpdom_data, 'e3u_0', e3u_0(:,:,:) )
559            CALL iom_get( inum4, jpdom_data, 'e3v_0', e3v_0(:,:,:) )
560            CALL iom_get( inum4, jpdom_data, 'e3w_0', e3w_0(:,:,:) )
561
562            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
563            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
564         ENDIF
565
566 
567         IF( ln_zps ) THEN                                           ! z-coordinate - partial steps
568            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d )  ! reference depth
569            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
570            CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors
571            CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
572            !
573            IF( nmsh <= 6 ) THEN                                        ! 3D vertical scale factors
574               CALL iom_get( inum4, jpdom_data, 'e3t_0', e3t_0(:,:,:) )
575               CALL iom_get( inum4, jpdom_data, 'e3u_0', e3u_0(:,:,:) )
576               CALL iom_get( inum4, jpdom_data, 'e3v_0', e3v_0(:,:,:) )
577               CALL iom_get( inum4, jpdom_data, 'e3w_0', e3w_0(:,:,:) )
578            ELSE                                                        ! 2D bottom scale factors
579               CALL iom_get( inum4, jpdom_data, 'e3t_ps', e3tp )
580               CALL iom_get( inum4, jpdom_data, 'e3w_ps', e3wp )
581               !                                                        ! deduces the 3D scale factors
582               DO jk = 1, jpk
583                  e3t_0(:,:,jk) = e3t_1d(jk)                                    ! set to the ref. factors
584                  e3u_0(:,:,jk) = e3t_1d(jk)
585                  e3v_0(:,:,jk) = e3t_1d(jk)
586                  e3w_0(:,:,jk) = e3w_1d(jk)
587               END DO
588               DO jj = 1,jpj                                                  ! adjust the deepest values
589                  DO ji = 1,jpi
590                     ik = mbkt(ji,jj)
591                     e3t_0(ji,jj,ik) = e3tp(ji,jj) * tmask(ji,jj,1) + e3t_1d(1) * ( 1._wp - tmask(ji,jj,1) )
592                     e3w_0(ji,jj,ik) = e3wp(ji,jj) * tmask(ji,jj,1) + e3w_1d(1) * ( 1._wp - tmask(ji,jj,1) )
593                  END DO
594               END DO
595               DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
596                  DO jj = 1, jpjm1
597                     DO ji = 1, jpim1
598                        e3u_0(ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
599                        e3v_0(ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
600                     END DO
601                  END DO
602               END DO
603               CALL lbc_lnk( e3u_0(:,:,:) , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0(:,:,:), 'U', 1._wp )   ! lateral boundary conditions
604               CALL lbc_lnk( e3v_0(:,:,:) , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0(:,:,:), 'V', 1._wp )
605               !
606               DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
607                  WHERE( e3u_0(:,:,jk) == 0._wp )   e3u_0(:,:,jk) = e3t_1d(jk)
608                  WHERE( e3v_0(:,:,jk) == 0._wp )   e3v_0(:,:,jk) = e3t_1d(jk)
609               END DO
610            END IF
611
612            IF( iom_varid( inum4, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level
613               CALL iom_get( inum4, jpdom_data, 'gdept_0', gdept_0(:,:,:) )
614               CALL iom_get( inum4, jpdom_data, 'gdepw_0', gdepw_0(:,:,:) )
615            ELSE                                                           ! 2D bottom depth
616               CALL iom_get( inum4, jpdom_data, 'hdept', zprt )
617               CALL iom_get( inum4, jpdom_data, 'hdepw', zprw )
618               !
619               DO jk = 1, jpk                                              ! deduces the 3D depth
620                  gdept_0(:,:,jk) = gdept_1d(jk)
621                  gdepw_0(:,:,jk) = gdepw_1d(jk)
622               END DO
623               DO jj = 1, jpj
624                  DO ji = 1, jpi
625                     ik = mbkt(ji,jj)
626                     IF( ik > 0 ) THEN
627                        gdepw_0(ji,jj,ik+1) = zprw(ji,jj)
628                        gdept_0(ji,jj,ik  ) = zprt(ji,jj)
629                        gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
630                     ENDIF
631                  END DO
632               END DO
633            ENDIF
634            !
635         ENDIF
636
637         IF( ln_zco ) THEN           ! Vertical coordinates and scales factors
638            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
639            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
640            CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )
641            CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
642            DO jk = 1, jpk
643               e3t_0(:,:,jk) = e3t_1d(jk)                              ! set to the ref. factors
644               e3u_0(:,:,jk) = e3t_1d(jk)
645               e3v_0(:,:,jk) = e3t_1d(jk)
646               e3w_0(:,:,jk) = e3w_1d(jk)
647               gdept_0(:,:,jk) = gdept_1d(jk)
648               gdepw_0(:,:,jk) = gdepw_1d(jk)
649            END DO
650         ENDIF
651
652      !
653      !              !==  time varying part of coordinate system  ==!
654      !
655      !       before        !          now          !       after         !
656      ;  gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
657      ;  gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
658      ;                     ;   gde3w_n = gde3w_0   !        ---          !
659      !
660      ;    e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
661      ;    e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
662      ;    e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
663      ;                     ;     e3f_n =   e3f_0   !        ---          !
664      ;    e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
665      ;   e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
666      ;   e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
667      !
668
669!!gm BUG in s-coordinate this does not work!
670      ! deepest/shallowest W level Above/Below ~10m
671      zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness)
672      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
673      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
674!!gm end bug
675
676      ! Control printing : Grid informations (if not restart)
677      ! ----------------
678
679      IF(lwp .AND. .NOT.ln_rstart ) THEN
680         WRITE(numout,*)
681         WRITE(numout,*) '          longitude and e1 scale factors'
682         WRITE(numout,*) '          ------------------------------'
683         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
684            glamv(ji,1), glamf(ji,1),   &
685            e1t(ji,1), e1u(ji,1),   &
686            e1v(ji,1), ji = 1, jpi,10)
6879300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
688            f19.10, 1x, f19.10, 1x, f19.10 )
689
690         WRITE(numout,*)
691         WRITE(numout,*) '          latitude and e2 scale factors'
692         WRITE(numout,*) '          -----------------------------'
693         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
694            &                     gphiv(1,jj), gphif(1,jj),   &
695            &                     e2t  (1,jj), e2u  (1,jj),   &
696            &                     e2v  (1,jj), jj = 1, jpj, 10 )
697      ENDIF
698
699      IF(lwp) THEN
700         WRITE(numout,*)
701         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
702         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
703         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
704      ENDIF
705
706      DO jk = 1, jpk
707         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' )
708         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' )
709      END DO
710      !                                     ! ============================
711      !                                     !        close the files
712      !                                     ! ============================
713      SELECT CASE ( nmsh )
714         CASE ( 1 )               
715            CALL iom_close( inum0 )
716         CASE ( 2 )
717            CALL iom_close( inum1 )
718            CALL iom_close( inum2 )
719         CASE ( 3 )
720            CALL iom_close( inum2 )
721            CALL iom_close( inum3 )
722            CALL iom_close( inum4 )
723      END SELECT
724      !
725      CALL wrk_dealloc( jpi, jpj, zmbk, zprt, zprw )
726      !
727   END SUBROUTINE dom_grd
728
729
730   SUBROUTINE zgr_bot_level
731      !!----------------------------------------------------------------------
732      !!                    ***  ROUTINE zgr_bot_level  ***
733      !!
734      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
735      !!
736      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
737      !!
738      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
739      !!                                     ocean level at t-, u- & v-points
740      !!                                     (min value = 1 over land)
741      !!----------------------------------------------------------------------
742      INTEGER ::   ji, jj   ! dummy loop indices
743      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
744      !!----------------------------------------------------------------------
745
746      !
747      IF(lwp) WRITE(numout,*)
748      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
749      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
750      !
751      CALL wrk_alloc( jpi, jpj, zmbk )
752      !
753      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
754      mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf)
755      !                                     ! bottom k-index of W-level = mbkt+1
756      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
757         DO ji = 1, jpim1
758            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
759            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
760         END DO
761      END DO
762      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
763      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
764      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
765      !
766      CALL wrk_dealloc( jpi, jpj, zmbk )
767      !
768   END SUBROUTINE zgr_bot_level
769
770
771   SUBROUTINE dom_msk
772      !!---------------------------------------------------------------------
773      !!                 ***  ROUTINE dom_msk  ***
774      !!
775      !! ** Purpose :   Off-line case: defines the interior domain T-mask.
776      !!
777      !! ** Method  :   The interior ocean/land mask is computed from tmask
778      !!              setting to zero the duplicated row and lines due to
779      !!              MPP exchange halos, est-west cyclic and north fold
780      !!              boundary conditions.
781      !!
782      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point
783      !!               tpol     : ???
784      !!----------------------------------------------------------------------
785      INTEGER  ::   ji, jj, jk           ! dummy loop indices
786      INTEGER  ::   iif, iil, ijf, ijl   ! local integers
787      INTEGER, POINTER, DIMENSION(:,:) ::  imsk 
788      !!---------------------------------------------------------------------
789     
790      CALL wrk_alloc( jpi, jpj, imsk )
791      !
792      ! Interior domain mask (used for global sum)
793      ! --------------------
794      ssmask(:,:)  = tmask(:,:,1)
795      tmask_i(:,:) = tmask(:,:,1)
796      iif = jpreci                        ! thickness of exchange halos in i-axis
797      iil = nlci - jpreci + 1
798      ijf = jprecj                        ! thickness of exchange halos in j-axis
799      ijl = nlcj - jprecj + 1
800      !
801      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns
802      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns)
803      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows
804      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows)
805      !
806      !                                   ! north fold mask
807      tpol(1:jpiglo) = 1._wp
808      !                               
809      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot
810      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot
811      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row
812         IF( mjg(ijl-1) == jpjglo-1 ) THEN
813            DO ji = iif+1, iil-1
814               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji))
815            END DO
816         ENDIF
817      ENDIF 
818      !
819      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at
820      ! least 1 wet u point
821      DO jj = 1, jpjm1
822         DO ji = 1, fs_jpim1   ! vector loop
823            ssumask(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
824            ssvmask(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
825         END DO
826         DO ji = 1, jpim1      ! NO vector opt.
827            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
828               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
829         END DO
830      END DO
831      CALL lbc_lnk( ssumask, 'U', 1._wp )      ! Lateral boundary conditions
832      CALL lbc_lnk( ssvmask, 'V', 1._wp )
833      CALL lbc_lnk( ssfmask, 'F', 1._wp )
834
835      ! 3. Ocean/land mask at wu-, wv- and w points
836      !----------------------------------------------
837      wmask (:,:,1) = tmask(:,:,1)        ! surface value
838      wumask(:,:,1) = umask(:,:,1) 
839      wvmask(:,:,1) = vmask(:,:,1)
840      DO jk = 2, jpk                      ! deeper value
841         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
842         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
843         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
844      END DO
845      !
846      CALL wrk_dealloc( jpi, jpj, imsk )
847      !
848   END SUBROUTINE dom_msk
849
850   !!======================================================================
851END MODULE domrea
852
Note: See TracBrowser for help on using the repository browser.