New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domrea.F90 in branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OFF_SRC – NEMO

source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 7293

Last change on this file since 7293 was 7158, checked in by clem, 8 years ago

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