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

source: branches/UKMO/r5936_hadgem3_cplfld/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 7138

Last change on this file since 7138 was 7138, checked in by jcastill, 7 years ago

Remove svn keywords

File size: 42.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 "domzgr_substitute.h90"
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OFF 3.7 , NEMO Consortium (2015)
40   !! $Id$
41   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE dom_rea
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE dom_rea  ***
48      !!                   
49      !! ** Purpose :   Domain initialization. Call the routines that are
50      !!      required to create the arrays which define the space and time
51      !!      domain of the ocean model.
52      !!
53      !! ** Method  :
54      !!      - dom_stp: defined the model time step
55      !!      - dom_rea: read the meshmask file if nmsh=1
56      !!----------------------------------------------------------------------
57      INTEGER ::   jk          ! dummy loop index
58      INTEGER ::   iconf = 0   ! local integers
59      !!----------------------------------------------------------------------
60      !
61      IF(lwp) THEN
62         WRITE(numout,*)
63         WRITE(numout,*) 'dom_init : domain initialization'
64         WRITE(numout,*) '~~~~~~~~'
65      ENDIF
66      !
67      CALL dom_nam      ! read namelist ( namrun, namdom )
68      CALL dom_zgr      ! Vertical mesh and bathymetry option
69      CALL dom_grd      ! Create a domain file
70      !
71      !                                      ! associated horizontal metrics
72      !
73      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:)
74      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:)
75      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:)
76      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:)
77      !
78      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:)
79      e1e2u (:,:) = e1u(:,:) * e2u(:,:)   ;   r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)
80      e1e2v (:,:) = e1v(:,:) * e2v(:,:)   ;   r1_e1e2v(:,:) = 1._wp / e1e2v(:,:)
81      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:)
82      !   
83      e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
84      e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
85      !
86      hu(:,:) = 0._wp                        ! Ocean depth at U- and V-points
87      hv(:,:) = 0._wp
88      DO jk = 1, jpk
89         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
90         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
91      END DO
92      !                                        ! Inverse of the local depth
93      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
94      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
95      !
96      CALL dom_stp      ! Time step
97      CALL dom_msk      ! Masks
98      CALL dom_ctl      ! Domain control
99      !
100   END SUBROUTINE dom_rea
101
102
103   SUBROUTINE dom_nam
104      !!----------------------------------------------------------------------
105      !!                     ***  ROUTINE dom_nam  ***
106      !!                   
107      !! ** Purpose :   read domaine namelists and print the variables.
108      !!
109      !! ** input   : - namrun namelist
110      !!              - namdom namelist
111      !!----------------------------------------------------------------------
112      USE ioipsl
113      INTEGER  ::   ios   ! Local integer output status for namelist read
114      !
115      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
116         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
117         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
118         &             nn_write, ln_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,   &
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', fse3t_n(:,:,:) ) ! scale factors
557            CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
558            CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
559            CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_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', fse3t_n(:,:,:) )
574               CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
575               CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
576               CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_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                  fse3t_n(:,:,jk) = e3t_1d(jk)                                    ! set to the ref. factors
583                  fse3u_n(:,:,jk) = e3t_1d(jk)
584                  fse3v_n(:,:,jk) = e3t_1d(jk)
585                  fse3w_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                     fse3t_n(ji,jj,ik) = e3tp(ji,jj) * tmask(ji,jj,1) + e3t_1d(1) * ( 1._wp - tmask(ji,jj,1) )
591                     fse3w_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                        fse3u_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji+1,jj,jk) )
598                        fse3v_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji,jj+1,jk) )
599                     END DO
600                  END DO
601               END DO
602               CALL lbc_lnk( fse3u_n(:,:,:) , 'U', 1._wp )   ;   CALL lbc_lnk( fse3uw_n(:,:,:), 'U', 1._wp )   ! lateral boundary conditions
603               CALL lbc_lnk( fse3v_n(:,:,:) , 'V', 1._wp )   ;   CALL lbc_lnk( fse3vw_n(:,:,:), 'V', 1._wp )
604               !
605               DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
606                  WHERE( fse3u_n(:,:,jk) == 0._wp )   fse3u_n(:,:,jk) = e3t_1d(jk)
607                  WHERE( fse3v_n(:,:,jk) == 0._wp )   fse3v_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', fsdept_n(:,:,:) )
613               CALL iom_get( inum4, jpdom_data, 'gdepw_0', fsdepw_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                  fsdept_n(:,:,jk) = gdept_1d(jk)
620                  fsdepw_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                        fsdepw_n(ji,jj,ik+1) = zprw(ji,jj)
627                        fsdept_n(ji,jj,ik  ) = zprt(ji,jj)
628                        fsdept_n(ji,jj,ik+1) = fsdept_n(ji,jj,ik) + fse3t_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               fse3t_n(:,:,jk) = e3t_1d(jk)                              ! set to the ref. factors
643               fse3u_n(:,:,jk) = e3t_1d(jk)
644               fse3v_n(:,:,jk) = e3t_1d(jk)
645               fse3w_n(:,:,jk) = e3w_1d(jk)
646               fsdept_n(:,:,jk) = gdept_1d(jk)
647               fsdepw_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
682      IF( nprint == 1 .AND. lwp ) THEN
683         WRITE(numout,*) '          e1u e2u '
684         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
685         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
686         WRITE(numout,*) '          e1v e2v  '
687         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
688         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
689      ENDIF
690
691      IF(lwp) THEN
692         WRITE(numout,*)
693         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
694         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
695         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
696      ENDIF
697
698      DO jk = 1, jpk
699         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' )
700         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' )
701      END DO
702      !                                     ! ============================
703      !                                     !        close the files
704      !                                     ! ============================
705      SELECT CASE ( nmsh )
706         CASE ( 1 )               
707            CALL iom_close( inum0 )
708         CASE ( 2 )
709            CALL iom_close( inum1 )
710            CALL iom_close( inum2 )
711         CASE ( 3 )
712            CALL iom_close( inum2 )
713            CALL iom_close( inum3 )
714            CALL iom_close( inum4 )
715      END SELECT
716      !
717      CALL wrk_dealloc( jpi, jpj, zmbk, zprt, zprw )
718      !
719   END SUBROUTINE dom_grd
720
721
722   SUBROUTINE zgr_bot_level
723      !!----------------------------------------------------------------------
724      !!                    ***  ROUTINE zgr_bot_level  ***
725      !!
726      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
727      !!
728      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
729      !!
730      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
731      !!                                     ocean level at t-, u- & v-points
732      !!                                     (min value = 1 over land)
733      !!----------------------------------------------------------------------
734      INTEGER ::   ji, jj   ! dummy loop indices
735      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
736      !!----------------------------------------------------------------------
737
738      !
739      IF(lwp) WRITE(numout,*)
740      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
741      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
742      !
743      CALL wrk_alloc( jpi, jpj, zmbk )
744      !
745      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
746      mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf)
747      !                                     ! bottom k-index of W-level = mbkt+1
748      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
749         DO ji = 1, jpim1
750            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
751            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
752         END DO
753      END DO
754      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
755      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
756      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
757      !
758      CALL wrk_dealloc( jpi, jpj, zmbk )
759      !
760   END SUBROUTINE zgr_bot_level
761
762
763   SUBROUTINE dom_msk
764      !!---------------------------------------------------------------------
765      !!                 ***  ROUTINE dom_msk  ***
766      !!
767      !! ** Purpose :   Off-line case: defines the interior domain T-mask.
768      !!
769      !! ** Method  :   The interior ocean/land mask is computed from tmask
770      !!              setting to zero the duplicated row and lines due to
771      !!              MPP exchange halos, est-west cyclic and north fold
772      !!              boundary conditions.
773      !!
774      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point
775      !!               tpol     : ???
776      !!----------------------------------------------------------------------
777      INTEGER  ::   ji, jj, jk           ! dummy loop indices
778      INTEGER  ::   iif, iil, ijf, ijl   ! local integers
779      INTEGER, POINTER, DIMENSION(:,:) ::  imsk 
780      !!---------------------------------------------------------------------
781     
782      CALL wrk_alloc( jpi, jpj, imsk )
783      !
784      ! Interior domain mask (used for global sum)
785      ! --------------------
786      ssmask(:,:)  = tmask(:,:,1)
787      tmask_i(:,:) = tmask(:,:,1)
788      iif = jpreci                        ! thickness of exchange halos in i-axis
789      iil = nlci - jpreci + 1
790      ijf = jprecj                        ! thickness of exchange halos in j-axis
791      ijl = nlcj - jprecj + 1
792      !
793      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns
794      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns)
795      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows
796      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows)
797      !
798      !                                   ! north fold mask
799      tpol(1:jpiglo) = 1._wp
800      !                               
801      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot
802      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot
803      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row
804         IF( mjg(ijl-1) == jpjglo-1 ) THEN
805            DO ji = iif+1, iil-1
806               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji))
807            END DO
808         ENDIF
809      ENDIF 
810      !
811      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at
812      ! least 1 wet u point
813      DO jj = 1, jpjm1
814         DO ji = 1, fs_jpim1   ! vector loop
815            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
816            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
817         END DO
818         DO ji = 1, jpim1      ! NO vector opt.
819            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
820               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
821         END DO
822      END DO
823      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
824      CALL lbc_lnk( vmask_i, 'V', 1._wp )
825      CALL lbc_lnk( fmask_i, 'F', 1._wp )
826
827      ! 3. Ocean/land mask at wu-, wv- and w points
828      !----------------------------------------------
829      wmask (:,:,1) = tmask(:,:,1)        ! surface value
830      wumask(:,:,1) = umask(:,:,1) 
831      wvmask(:,:,1) = vmask(:,:,1)
832      DO jk = 2, jpk                      ! deeper value
833         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
834         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
835         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
836      END DO
837      !
838      IF( nprint == 1 .AND. lwp ) THEN    ! Control print
839         imsk(:,:) = INT( tmask_i(:,:) )
840         WRITE(numout,*) ' tmask_i : '
841         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
842         WRITE (numout,*)
843         WRITE (numout,*) ' dommsk: tmask for each level'
844         WRITE (numout,*) ' ----------------------------'
845         DO jk = 1, jpk
846            imsk(:,:) = INT( tmask(:,:,jk) )
847            WRITE(numout,*)
848            WRITE(numout,*) ' level = ',jk
849            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
850         END DO
851      ENDIF
852      !
853      CALL wrk_dealloc( jpi, jpj, imsk )
854      !
855   END SUBROUTINE dom_msk
856
857   !!======================================================================
858END MODULE domrea
859
Note: See TracBrowser for help on using the repository browser.