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

source: branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 6808

Last change on this file since 6808 was 6808, checked in by jamesharle, 8 years ago

merge with trunk@6232 for consistency with SSB code

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