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

source: trunk/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 5504

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

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