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

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

source: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OFF_SRC/domrea.F90 @ 5737

Last change on this file since 5737 was 5737, checked in by gm, 9 years ago

#1593: LDF-ADV, step I: Phasing of horizontal scale factors correct 2

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