source: branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 4007

Last change on this file since 4007 was 4007, checked in by davestorkey, 8 years ago
  1. Bug fixes for flagu/flagv calculation in bdyini.F90.
  2. Introduce masking of derivatives in radiation velocity calculation in bdylib.F90
  3. Change relaxation term in Orlanski implementation to explicit timestepping in bdylib.F90.
  4. Remove bdyfmask (redundant).
  • Property svn:keywords set to Id
File size: 80.9 KB
Line 
1MODULE bdyini
2   !!======================================================================
3   !!                       ***  MODULE  bdyini  ***
4   !! Unstructured open boundaries : initialisation
5   !!======================================================================
6   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
7   !!             -   !  2007-01  (D. Storkey) Update to use IOM module
8   !!             -   !  2007-01  (D. Storkey) Tidal forcing
9   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
10   !!            3.3  !  2010-09  (E.O'Dea) updates for Shelf configurations
11   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
12   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
13   !!            3.4  !  2012     (J. Chanut) straight open boundary case update
14   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Updates for the
15   !!                             optimization of BDY communications
16   !!----------------------------------------------------------------------
17#if defined key_bdy
18   !!----------------------------------------------------------------------
19   !!   'key_bdy'                     Unstructured Open Boundary Conditions
20   !!----------------------------------------------------------------------
21   !!   bdy_init       : Initialization of unstructured open boundaries
22   !!----------------------------------------------------------------------
23   USE wrk_nemo        ! Memory Allocation
24   USE timing          ! Timing
25   USE oce             ! ocean dynamics and tracers variables
26   USE dom_oce         ! ocean space and time domain
27   USE bdy_oce         ! unstructured open boundary conditions
28   USE in_out_manager  ! I/O units
29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
30   USE lib_mpp         ! for mpp_sum 
31   USE iom             ! I/O
32   USE sbctide, ONLY: lk_tide ! Tidal forcing or not
33   USE phycst, ONLY: rday
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   bdy_init   ! routine called in nemo_init
39
40   INTEGER, PARAMETER          :: jp_nseg = 100
41   INTEGER, PARAMETER          :: nrimmax = 20 ! maximum rimwidth in structured
42                                               ! open boundary data files
43   ! Straight open boundary segment parameters:
44   INTEGER  :: nbdysege, nbdysegw, nbdysegn, nbdysegs 
45   INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge
46   INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw
47   INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn
48   INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
51   !! $Id$
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55   
56   SUBROUTINE bdy_init
57      !!----------------------------------------------------------------------
58      !!                 ***  ROUTINE bdy_init  ***
59      !!         
60      !! ** Purpose :   Initialization of the dynamics and tracer fields with
61      !!              unstructured open boundaries.
62      !!
63      !! ** Method  :   Read initialization arrays (mask, indices) to identify
64      !!              an unstructured open boundary
65      !!
66      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
67      !!----------------------------------------------------------------------     
68      ! namelist variables
69      !-------------------
70      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile
71      CHARACTER(LEN=1)   ::   ctypebdy
72      INTEGER :: nbdyind, nbdybeg, nbdyend
73
74      ! local variables
75      !-------------------
76      INTEGER  ::   ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
77      INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers
78      INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       -
79      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       -
80      INTEGER  ::   jpbdtau, jpbdtas                       !   -       -
81      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       -
82      INTEGER  ::   i_offset, j_offset                     !   -       -
83      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts
84      REAL(wp), POINTER  ::  flagu, flagv                  !    -   -
85      REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask    ! pointer to 2D mask fields
86      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
87      INTEGER, DIMENSION (2)                  ::   kdimsz
88      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays
89      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
90      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
91      CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid
92      INTEGER :: com_east, com_west, com_south, com_north          ! Flags for boundaries sending
93      INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b  ! Flags for boundaries receiving
94      INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4)                ! Arrays for neighbours coordinates
95      REAL(wp), POINTER, DIMENSION(:,:)       ::   zfmask  ! temporary fmask array excluding coastal boundary condition (shlat)
96
97      !!
98      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,                 &
99         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
100         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             & 
101         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
102#if defined key_lim2
103         &             cn_ice_lim2, nn_ice_lim2_dta,                           &
104#endif
105         &             ln_vol, nn_volctl, nn_rimwidth
106      !!
107      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
108
109      !!----------------------------------------------------------------------
110
111      IF( nn_timing == 1 ) CALL timing_start('bdy_init')
112
113      IF( bdy_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' )
114
115      IF(lwp) WRITE(numout,*)
116      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
117      IF(lwp) WRITE(numout,*) '~~~~~~~~'
118      !
119
120      IF( jperio /= 0 )   CALL ctl_stop( 'Cyclic or symmetric,',   &
121         &                               ' and general open boundary condition are not compatible' )
122
123      cgrid= (/'t','u','v'/)
124     
125      ! -----------------------------------------
126      ! Initialise and read namelist parameters
127      ! -----------------------------------------
128
129      nb_bdy            = 0
130      ln_coords_file(:) = .false.
131      cn_coords_file(:) = ''
132      ln_mask_file      = .false.
133      cn_mask_file(:)   = ''
134      cn_dyn2d(:)       = ''
135      nn_dyn2d_dta(:)   = -1  ! uninitialised flag
136      cn_dyn3d(:)       = ''
137      nn_dyn3d_dta(:)   = -1  ! uninitialised flag
138      cn_tra(:)         = ''
139      nn_tra_dta(:)     = -1  ! uninitialised flag
140      ln_tra_dmp(:)     = .false.
141      ln_dyn3d_dmp(:)   = .false.
142      rn_time_dmp(:)    = 1.
143#if defined key_lim2
144      cn_ice_lim2(:)    = ''
145      nn_ice_lim2_dta(:)= -1  ! uninitialised flag
146#endif
147      ln_vol            = .false.
148      nn_volctl         = -1  ! uninitialised flag
149      nn_rimwidth(:)    = -1  ! uninitialised flag
150
151      REWIND( numnam )                   
152      READ  ( numnam, nambdy )
153
154      ! -----------------------------------------
155      ! Check and write out namelist parameters
156      ! -----------------------------------------
157      !                                   ! control prints
158      IF(lwp) WRITE(numout,*) '         nambdy'
159
160      IF( nb_bdy .eq. 0 ) THEN
161        IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
162      ELSE
163        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy
164      ENDIF
165
166      DO ib_bdy = 1,nb_bdy
167        IF(lwp) WRITE(numout,*) ' ' 
168        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' 
169
170        IF( ln_coords_file(ib_bdy) ) THEN
171           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
172        ELSE
173           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
174        ENDIF
175        IF(lwp) WRITE(numout,*)
176
177        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
178        SELECT CASE( cn_dyn2d(ib_bdy) )                 
179          CASE('none')         
180             IF(lwp) WRITE(numout,*) '      no open boundary condition'       
181             dta_bdy(ib_bdy)%ll_ssh = .false.
182             dta_bdy(ib_bdy)%ll_u2d = .false.
183             dta_bdy(ib_bdy)%ll_v2d = .false.
184          CASE('frs')         
185             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
186             dta_bdy(ib_bdy)%ll_ssh = .false.
187             dta_bdy(ib_bdy)%ll_u2d = .true.
188             dta_bdy(ib_bdy)%ll_v2d = .true.
189          CASE('flather')     
190             IF(lwp) WRITE(numout,*) '      Flather radiation condition'
191             dta_bdy(ib_bdy)%ll_ssh = .true.
192             dta_bdy(ib_bdy)%ll_u2d = .true.
193             dta_bdy(ib_bdy)%ll_v2d = .true.
194          CASE('orlanski')     
195             IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
196             dta_bdy(ib_bdy)%ll_ssh = .false.
197             dta_bdy(ib_bdy)%ll_u2d = .true.
198             dta_bdy(ib_bdy)%ll_v2d = .true.
199          CASE('orlanski_npo') 
200             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
201             dta_bdy(ib_bdy)%ll_ssh = .false.
202             dta_bdy(ib_bdy)%ll_u2d = .true.
203             dta_bdy(ib_bdy)%ll_v2d = .true.
204          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' )
205        END SELECT
206        IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN
207           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !
208              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
209              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
210              CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file'
211              CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files'
212              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
213           END SELECT
214           IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.lk_tide)) THEN
215             CALL ctl_stop( 'You must activate key_tide to add tidal forcing at open boundaries' )
216           ENDIF
217        ENDIF
218        IF(lwp) WRITE(numout,*)
219
220        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
221        SELECT CASE( cn_dyn3d(ib_bdy) )                 
222          CASE('none')
223             IF(lwp) WRITE(numout,*) '      no open boundary condition'       
224             dta_bdy(ib_bdy)%ll_u3d = .false.
225             dta_bdy(ib_bdy)%ll_v3d = .false.
226          CASE('frs')       
227             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
228             dta_bdy(ib_bdy)%ll_u3d = .true.
229             dta_bdy(ib_bdy)%ll_v3d = .true.
230          CASE('specified')
231             IF(lwp) WRITE(numout,*) '      Specified value'
232             dta_bdy(ib_bdy)%ll_u3d = .true.
233             dta_bdy(ib_bdy)%ll_v3d = .true.
234          CASE('zero')
235             IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
236             dta_bdy(ib_bdy)%ll_u3d = .false.
237             dta_bdy(ib_bdy)%ll_v3d = .false.
238          CASE('orlanski')
239             IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
240             dta_bdy(ib_bdy)%ll_u3d = .true.
241             dta_bdy(ib_bdy)%ll_v3d = .true.
242          CASE('orlanski_npo')
243             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
244             dta_bdy(ib_bdy)%ll_u3d = .true.
245             dta_bdy(ib_bdy)%ll_v3d = .true.
246          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn3d' )
247        END SELECT
248        IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
249           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !
250              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
251              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
252              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
253           END SELECT
254        ENDIF
255
256        IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
257           IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN
258              IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
259              ln_dyn3d_dmp(ib_bdy)=.false.
260           ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN
261              CALL ctl_stop( 'Use FRS OR relaxation' )
262           ELSE
263              IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone'
264              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
265              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
266              dta_bdy(ib_bdy)%ll_u3d = .true.
267              dta_bdy(ib_bdy)%ll_v3d = .true.
268           ENDIF
269        ELSE
270           IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities'
271        ENDIF
272        IF(lwp) WRITE(numout,*)
273
274        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
275        SELECT CASE( cn_tra(ib_bdy) )                 
276          CASE('none')
277             IF(lwp) WRITE(numout,*) '      no open boundary condition'       
278             dta_bdy(ib_bdy)%ll_tem = .false.
279             dta_bdy(ib_bdy)%ll_sal = .false.
280          CASE('frs')
281             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
282             dta_bdy(ib_bdy)%ll_tem = .true.
283             dta_bdy(ib_bdy)%ll_sal = .true.
284          CASE('specified')
285             IF(lwp) WRITE(numout,*) '      Specified value'
286             dta_bdy(ib_bdy)%ll_tem = .true.
287             dta_bdy(ib_bdy)%ll_sal = .true.
288          CASE('neumann')
289             IF(lwp) WRITE(numout,*) '      Neumann conditions'
290             dta_bdy(ib_bdy)%ll_tem = .false.
291             dta_bdy(ib_bdy)%ll_sal = .false.
292          CASE('runoff')
293             IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity'
294             dta_bdy(ib_bdy)%ll_tem = .false.
295             dta_bdy(ib_bdy)%ll_sal = .false.
296          CASE('orlanski')
297             IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
298             dta_bdy(ib_bdy)%ll_tem = .true.
299             dta_bdy(ib_bdy)%ll_sal = .true.
300          CASE('orlanski_npo')
301             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
302             dta_bdy(ib_bdy)%ll_tem = .true.
303             dta_bdy(ib_bdy)%ll_sal = .true.
304          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_tra' )
305        END SELECT
306        IF( cn_tra(ib_bdy) /= 'none' ) THEN
307           SELECT CASE( nn_tra_dta(ib_bdy) )                   !
308              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
309              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
310              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
311           END SELECT
312        ENDIF
313
314        IF ( ln_tra_dmp(ib_bdy) ) THEN
315           IF ( cn_tra(ib_bdy) == 'none' ) THEN
316              IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
317              ln_tra_dmp(ib_bdy)=.false.
318           ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN
319              CALL ctl_stop( 'Use FRS OR relaxation' )
320           ELSE
321              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
322              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
323              IF(lwp) WRITE(numout,*) '      Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'
324              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
325              dta_bdy(ib_bdy)%ll_tem = .true.
326              dta_bdy(ib_bdy)%ll_sal = .true.
327           ENDIF
328        ELSE
329           IF(lwp) WRITE(numout,*) '      NO T/S relaxation'
330        ENDIF
331        IF(lwp) WRITE(numout,*)
332
333#if defined key_lim2
334        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
335        SELECT CASE( cn_ice_lim2(ib_bdy) )                 
336          CASE('none')
337             IF(lwp) WRITE(numout,*) '      no open boundary condition'       
338             dta_bdy(ib_bdy)%ll_frld  = .false.
339             dta_bdy(ib_bdy)%ll_hicif = .false.
340             dta_bdy(ib_bdy)%ll_hsnif = .false.
341          CASE('frs')
342             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
343             dta_bdy(ib_bdy)%ll_frld  = .true.
344             dta_bdy(ib_bdy)%ll_hicif = .true.
345             dta_bdy(ib_bdy)%ll_hsnif = .true.
346          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice_lim2' )
347        END SELECT
348        IF( cn_ice_lim2(ib_bdy) /= 'none' ) THEN
349           SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   !
350              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
351              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
352              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' )
353           END SELECT
354        ENDIF
355        IF(lwp) WRITE(numout,*)
356#endif
357
358        IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
359        IF(lwp) WRITE(numout,*)
360
361      ENDDO
362
363     IF (nb_bdy .gt. 0) THEN
364        IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
365          IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
366          IF(lwp) WRITE(numout,*)
367          SELECT CASE ( nn_volctl )
368            CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
369            CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
370            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
371          END SELECT
372          IF(lwp) WRITE(numout,*)
373        ELSE
374          IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
375          IF(lwp) WRITE(numout,*)
376        ENDIF
377     ENDIF
378
379      ! -------------------------------------------------
380      ! Initialise indices arrays for open boundaries
381      ! -------------------------------------------------
382
383      ! Work out global dimensions of boundary data
384      ! ---------------------------------------------
385      REWIND( numnam )                   
386               
387      nblendta(:,:) = 0
388      nbdysege = 0
389      nbdysegw = 0
390      nbdysegn = 0
391      nbdysegs = 0
392      icount   = 0 ! count user defined segments
393      ! Dimensions below are used to allocate arrays to read external data
394      jpbdtas = 1 ! Maximum size of boundary data (structured case)
395      jpbdtau = 1 ! Maximum size of boundary data (unstructured case)
396
397      DO ib_bdy = 1, nb_bdy
398
399         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
400 
401            icount = icount + 1
402            ! No REWIND here because may need to read more than one nambdy_index namelist.
403            READ  ( numnam, nambdy_index )
404
405            SELECT CASE ( TRIM(ctypebdy) )
406              CASE( 'N' )
407                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
408                    nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain.
409                    nbdybeg  = 2
410                    nbdyend  = jpiglo - 1
411                 ENDIF
412                 nbdysegn = nbdysegn + 1
413                 npckgn(nbdysegn) = ib_bdy ! Save bdy package number
414                 jpjnob(nbdysegn) = nbdyind
415                 jpindt(nbdysegn) = nbdybeg
416                 jpinft(nbdysegn) = nbdyend
417                 !
418              CASE( 'S' )
419                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
420                    nbdyind  = 2           ! set boundary to whole side of model domain.
421                    nbdybeg  = 2
422                    nbdyend  = jpiglo - 1
423                 ENDIF
424                 nbdysegs = nbdysegs + 1
425                 npckgs(nbdysegs) = ib_bdy ! Save bdy package number
426                 jpjsob(nbdysegs) = nbdyind
427                 jpisdt(nbdysegs) = nbdybeg
428                 jpisft(nbdysegs) = nbdyend
429                 !
430              CASE( 'E' )
431                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
432                    nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain.
433                    nbdybeg  = 2
434                    nbdyend  = jpjglo - 1
435                 ENDIF
436                 nbdysege = nbdysege + 1 
437                 npckge(nbdysege) = ib_bdy ! Save bdy package number
438                 jpieob(nbdysege) = nbdyind
439                 jpjedt(nbdysege) = nbdybeg
440                 jpjeft(nbdysege) = nbdyend
441                 !
442              CASE( 'W' )
443                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
444                    nbdyind  = 2           ! set boundary to whole side of model domain.
445                    nbdybeg  = 2
446                    nbdyend  = jpjglo - 1
447                 ENDIF
448                 nbdysegw = nbdysegw + 1
449                 npckgw(nbdysegw) = ib_bdy ! Save bdy package number
450                 jpiwob(nbdysegw) = nbdyind
451                 jpjwdt(nbdysegw) = nbdybeg
452                 jpjwft(nbdysegw) = nbdyend
453                 !
454              CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
455            END SELECT
456
457            ! For simplicity we assume that in case of straight bdy, arrays have the same length
458            ! (even if it is true that last tangential velocity points
459            ! are useless). This simplifies a little bit boundary data format (and agrees with format
460            ! used so far in obc package)
461
462            nblendta(1:jpbgrd,ib_bdy) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy)
463            jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1))
464            IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) &
465            & CALL ctl_stop( 'rimwidth must be lower than nrimmax' )
466
467         ELSE            ! Read size of arrays in boundary coordinates file.
468            CALL iom_open( cn_coords_file(ib_bdy), inum )
469            DO igrd = 1, jpbgrd
470               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 
471               nblendta(igrd,ib_bdy) = kdimsz(1)
472               jpbdtau = MAX(jpbdtau, kdimsz(1))
473            ENDDO
474            CALL iom_close( inum )
475
476         ENDIF
477
478      ENDDO ! ib_bdy
479
480      IF (nb_bdy>0) THEN
481         jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
482
483         ! Allocate arrays
484         !---------------
485         ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
486            &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
487
488         ALLOCATE( dta_global(jpbdtau, 1, jpk) )
489         IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) )
490         !
491      ENDIF
492
493      ! Now look for crossings in user (namelist) defined open boundary segments:
494      !--------------------------------------------------------------------------
495      IF ( icount>0 ) CALL bdy_ctl_seg
496
497      ! Calculate global boundary index arrays or read in from file
498      !------------------------------------------------------------               
499      ! 1. Read global index arrays from boundary coordinates file.
500      DO ib_bdy = 1, nb_bdy
501
502         IF( ln_coords_file(ib_bdy) ) THEN
503
504            CALL iom_open( cn_coords_file(ib_bdy), inum )
505            DO igrd = 1, jpbgrd
506               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
507               DO ii = 1,nblendta(igrd,ib_bdy)
508                  nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
509               END DO
510               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
511               DO ii = 1,nblendta(igrd,ib_bdy)
512                  nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
513               END DO
514               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
515               DO ii = 1,nblendta(igrd,ib_bdy)
516                  nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
517               END DO
518
519               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
520               IF(lwp) WRITE(numout,*)
521               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
522               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
523               IF (ibr_max < nn_rimwidth(ib_bdy))   &
524                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
525            END DO
526            CALL iom_close( inum )
527
528         ENDIF
529
530      ENDDO     
531   
532      ! 2. Now fill indices corresponding to straight open boundary arrays:
533      ! East
534      !-----
535      DO iseg = 1, nbdysege
536         ib_bdy = npckge(iseg)
537         !
538         ! ------------ T points -------------
539         igrd=1
540         icount=0
541         DO ir = 1, nn_rimwidth(ib_bdy)
542            DO ij = jpjedt(iseg), jpjeft(iseg)
543               icount = icount + 1
544               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
545               nbjdta(icount, igrd, ib_bdy) = ij
546               nbrdta(icount, igrd, ib_bdy) = ir
547            ENDDO
548         ENDDO
549         !
550         ! ------------ U points -------------
551         igrd=2
552         icount=0
553         DO ir = 1, nn_rimwidth(ib_bdy)
554            DO ij = jpjedt(iseg), jpjeft(iseg)
555               icount = icount + 1
556               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir
557               nbjdta(icount, igrd, ib_bdy) = ij
558               nbrdta(icount, igrd, ib_bdy) = ir
559            ENDDO
560         ENDDO
561         !
562         ! ------------ V points -------------
563         igrd=3
564         icount=0
565         DO ir = 1, nn_rimwidth(ib_bdy)
566!            DO ij = jpjedt(iseg), jpjeft(iseg) - 1
567            DO ij = jpjedt(iseg), jpjeft(iseg)
568               icount = icount + 1
569               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
570               nbjdta(icount, igrd, ib_bdy) = ij
571               nbrdta(icount, igrd, ib_bdy) = ir
572            ENDDO
573            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
574            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
575         ENDDO
576      ENDDO
577      !
578      ! West
579      !-----
580      DO iseg = 1, nbdysegw
581         ib_bdy = npckgw(iseg)
582         !
583         ! ------------ T points -------------
584         igrd=1
585         icount=0
586         DO ir = 1, nn_rimwidth(ib_bdy)
587            DO ij = jpjwdt(iseg), jpjwft(iseg)
588               icount = icount + 1
589               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
590               nbjdta(icount, igrd, ib_bdy) = ij
591               nbrdta(icount, igrd, ib_bdy) = ir
592            ENDDO
593         ENDDO
594         !
595         ! ------------ U points -------------
596         igrd=2
597         icount=0
598         DO ir = 1, nn_rimwidth(ib_bdy)
599            DO ij = jpjwdt(iseg), jpjwft(iseg)
600               icount = icount + 1
601               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
602               nbjdta(icount, igrd, ib_bdy) = ij
603               nbrdta(icount, igrd, ib_bdy) = ir
604            ENDDO
605         ENDDO
606         !
607         ! ------------ V points -------------
608         igrd=3
609         icount=0
610         DO ir = 1, nn_rimwidth(ib_bdy)
611!            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
612            DO ij = jpjwdt(iseg), jpjwft(iseg)
613               icount = icount + 1
614               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
615               nbjdta(icount, igrd, ib_bdy) = ij
616               nbrdta(icount, igrd, ib_bdy) = ir
617            ENDDO
618            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
619            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
620         ENDDO
621      ENDDO
622      !
623      ! North
624      !-----
625      DO iseg = 1, nbdysegn
626         ib_bdy = npckgn(iseg)
627         !
628         ! ------------ T points -------------
629         igrd=1
630         icount=0
631         DO ir = 1, nn_rimwidth(ib_bdy)
632            DO ii = jpindt(iseg), jpinft(iseg)
633               icount = icount + 1
634               nbidta(icount, igrd, ib_bdy) = ii
635               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
636               nbrdta(icount, igrd, ib_bdy) = ir
637            ENDDO
638         ENDDO
639         !
640         ! ------------ U points -------------
641         igrd=2
642         icount=0
643         DO ir = 1, nn_rimwidth(ib_bdy)
644!            DO ii = jpindt(iseg), jpinft(iseg) - 1
645            DO ii = jpindt(iseg), jpinft(iseg)
646               icount = icount + 1
647               nbidta(icount, igrd, ib_bdy) = ii
648               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
649               nbrdta(icount, igrd, ib_bdy) = ir
650            ENDDO
651            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
652            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
653         ENDDO
654         !
655         ! ------------ V points -------------
656         igrd=3
657         icount=0
658         DO ir = 1, nn_rimwidth(ib_bdy)
659            DO ii = jpindt(iseg), jpinft(iseg)
660               icount = icount + 1
661               nbidta(icount, igrd, ib_bdy) = ii
662               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir
663               nbrdta(icount, igrd, ib_bdy) = ir
664            ENDDO
665         ENDDO
666      ENDDO
667      !
668      ! South
669      !-----
670      DO iseg = 1, nbdysegs
671         ib_bdy = npckgs(iseg)
672         !
673         ! ------------ T points -------------
674         igrd=1
675         icount=0
676         DO ir = 1, nn_rimwidth(ib_bdy)
677            DO ii = jpisdt(iseg), jpisft(iseg)
678               icount = icount + 1
679               nbidta(icount, igrd, ib_bdy) = ii
680               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
681               nbrdta(icount, igrd, ib_bdy) = ir
682            ENDDO
683         ENDDO
684         !
685         ! ------------ U points -------------
686         igrd=2
687         icount=0
688         DO ir = 1, nn_rimwidth(ib_bdy)
689!            DO ii = jpisdt(iseg), jpisft(iseg) - 1
690            DO ii = jpisdt(iseg), jpisft(iseg)
691               icount = icount + 1
692               nbidta(icount, igrd, ib_bdy) = ii
693               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
694               nbrdta(icount, igrd, ib_bdy) = ir
695            ENDDO
696            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
697            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
698         ENDDO
699         !
700         ! ------------ V points -------------
701         igrd=3
702         icount=0
703         DO ir = 1, nn_rimwidth(ib_bdy)
704            DO ii = jpisdt(iseg), jpisft(iseg)
705               icount = icount + 1
706               nbidta(icount, igrd, ib_bdy) = ii
707               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
708               nbrdta(icount, igrd, ib_bdy) = ir
709            ENDDO
710         ENDDO
711      ENDDO
712
713      !  Deal with duplicated points
714      !-----------------------------
715      ! We assign negative indices to duplicated points (to remove them from bdy points to be updated)
716      ! if their distance to the bdy is greater than the other
717      ! If their distance are the same, just keep only one to avoid updating a point twice
718      DO igrd = 1, jpbgrd
719         DO ib_bdy1 = 1, nb_bdy
720            DO ib_bdy2 = 1, nb_bdy
721               IF (ib_bdy1/=ib_bdy2) THEN
722                  DO ib1 = 1, nblendta(igrd,ib_bdy1)
723                     DO ib2 = 1, nblendta(igrd,ib_bdy2)
724                        IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
725                        &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
726!                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &
727!                                                       &              nbidta(ib1, igrd, ib_bdy1),      &
728!                                                       &              nbjdta(ib2, igrd, ib_bdy2)
729                           ! keep only points with the lowest distance to boundary:
730                           IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN
731                             nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
732                             nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
733                           ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN
734                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
735                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
736                           ! Arbitrary choice if distances are the same:
737                           ELSE
738                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
739                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
740                           ENDIF
741                        END IF
742                     END DO
743                  END DO
744               ENDIF
745            END DO
746         END DO
747      END DO
748
749      ! Work out dimensions of boundary data on each processor
750      ! ------------------------------------------------------
751
752      ! Rather assume that boundary data indices are given on global domain
753      ! TO BE DISCUSSED ?
754!      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
755!      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
756!      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
757!      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1     
758      iw = mig(1) - jpizoom + 2         ! if monotasking and no zoom, iw=2
759      ie = mig(1) + nlci - jpizoom - 1  ! if monotasking and no zoom, ie=jpim1
760      is = mjg(1) - jpjzoom + 2         ! if monotasking and no zoom, is=2
761      in = mjg(1) + nlcj - jpjzoom - 1  ! if monotasking and no zoom, in=jpjm1
762
763      ALLOCATE( nbondi_bdy(nb_bdy))
764      ALLOCATE( nbondj_bdy(nb_bdy))
765      nbondi_bdy(:)=2
766      nbondj_bdy(:)=2
767      ALLOCATE( nbondi_bdy_b(nb_bdy))
768      ALLOCATE( nbondj_bdy_b(nb_bdy))
769      nbondi_bdy_b(:)=2
770      nbondj_bdy_b(:)=2
771
772      ! Work out dimensions of boundary data on each neighbour process
773      IF(nbondi .eq. 0) THEN
774         iw_b(1) = jpizoom + nimppt(nowe+1)
775         ie_b(1) = jpizoom + nimppt(nowe+1)+nlcit(nowe+1)-3
776         is_b(1) = jpjzoom + njmppt(nowe+1)
777         in_b(1) = jpjzoom + njmppt(nowe+1)+nlcjt(nowe+1)-3
778
779         iw_b(2) = jpizoom + nimppt(noea+1)
780         ie_b(2) = jpizoom + nimppt(noea+1)+nlcit(noea+1)-3
781         is_b(2) = jpjzoom + njmppt(noea+1)
782         in_b(2) = jpjzoom + njmppt(noea+1)+nlcjt(noea+1)-3
783      ELSEIF(nbondi .eq. 1) THEN
784         iw_b(1) = jpizoom + nimppt(nowe+1)
785         ie_b(1) = jpizoom + nimppt(nowe+1)+nlcit(nowe+1)-3
786         is_b(1) = jpjzoom + njmppt(nowe+1)
787         in_b(1) = jpjzoom + njmppt(nowe+1)+nlcjt(nowe+1)-3
788      ELSEIF(nbondi .eq. -1) THEN
789         iw_b(2) = jpizoom + nimppt(noea+1)
790         ie_b(2) = jpizoom + nimppt(noea+1)+nlcit(noea+1)-3
791         is_b(2) = jpjzoom + njmppt(noea+1)
792         in_b(2) = jpjzoom + njmppt(noea+1)+nlcjt(noea+1)-3
793      ENDIF
794
795      IF(nbondj .eq. 0) THEN
796         iw_b(3) = jpizoom + nimppt(noso+1)
797         ie_b(3) = jpizoom + nimppt(noso+1)+nlcit(noso+1)-3
798         is_b(3) = jpjzoom + njmppt(noso+1)
799         in_b(3) = jpjzoom + njmppt(noso+1)+nlcjt(noso+1)-3
800
801         iw_b(4) = jpizoom + nimppt(nono+1)
802         ie_b(4) = jpizoom + nimppt(nono+1)+nlcit(nono+1)-3
803         is_b(4) = jpjzoom + njmppt(nono+1)
804         in_b(4) = jpjzoom + njmppt(nono+1)+nlcjt(nono+1)-3
805      ELSEIF(nbondj .eq. 1) THEN
806         iw_b(3) = jpizoom + nimppt(noso+1)
807         ie_b(3) = jpizoom + nimppt(noso+1)+nlcit(noso+1)-3
808         is_b(3) = jpjzoom + njmppt(noso+1)
809         in_b(3) = jpjzoom + njmppt(noso+1)+nlcjt(noso+1)-3
810      ELSEIF(nbondj .eq. -1) THEN
811         iw_b(4) = jpizoom + nimppt(nono+1)
812         ie_b(4) = jpizoom + nimppt(nono+1)+nlcit(nono+1)-3
813         is_b(4) = jpjzoom + njmppt(nono+1)
814         in_b(4) = jpjzoom + njmppt(nono+1)+nlcjt(nono+1)-3
815      ENDIF
816
817      DO ib_bdy = 1, nb_bdy
818         DO igrd = 1, jpbgrd
819            icount  = 0
820            icountr = 0
821            idx_bdy(ib_bdy)%nblen(igrd)    = 0
822            idx_bdy(ib_bdy)%nblenrim(igrd) = 0
823            DO ib = 1, nblendta(igrd,ib_bdy)
824               ! check that data is in correct order in file
825               ibm1 = MAX(1,ib-1)
826               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc...
827                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
828                     CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', &
829                                   'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory')
830                  ENDIF   
831               ENDIF
832               ! check if point is in local domain
833               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
834                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN
835                  !
836                  icount = icount  + 1
837                  !
838                  IF( nbrdta(ib,igrd,ib_bdy) == 1 )   icountr = icountr+1
839               ENDIF
840            ENDDO
841            idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
842            idx_bdy(ib_bdy)%nblen   (igrd) = icount  !: length of boundary data on each proc       
843         ENDDO  ! igrd
844
845         ! Allocate index arrays for this boundary set
846         !--------------------------------------------
847         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
848         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
849         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
850         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
851         ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) )
852         ALLOCATE( idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) )
853         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
854         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
855         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1,jpbgrd) )
856         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1,jpbgrd) )
857
858         ! Dispatch mapping indices and discrete distances on each processor
859         ! -----------------------------------------------------------------
860
861         com_east = 0
862         com_west = 0
863         com_south = 0
864         com_north = 0
865
866         com_east_b = 0
867         com_west_b = 0
868         com_south_b = 0
869         com_north_b = 0
870         DO igrd = 1, jpbgrd
871            icount  = 0
872            ! Loop on rimwidth to ensure outermost points come first in the local arrays.
873            DO ir=1, nn_rimwidth(ib_bdy)
874               DO ib = 1, nblendta(igrd,ib_bdy)
875                  ! check if point is in local domain and equals ir
876                  IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
877                     & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   &
878                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
879                     !
880                     icount = icount  + 1
881
882                     ! Rather assume that boundary data indices are given on global domain
883                     ! TO BE DISCUSSED ?
884!                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1
885!                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
886                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+jpizoom
887                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+jpjzoom
888                     ! check if point has to be sent
889                     ii = idx_bdy(ib_bdy)%nbi(icount,igrd)
890                     ij = idx_bdy(ib_bdy)%nbj(icount,igrd)
891                     if((com_east .ne. 1) .and. (ii .eq. (nlci-1)) .and. (nbondi .le. 0)) then
892                        com_east = 1
893                     elseif((com_west .ne. 1) .and. (ii .eq. 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then
894                        com_west = 1
895                     endif
896                     if((com_south .ne. 1) .and. (ij .eq. 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then
897                        com_south = 1
898                     elseif((com_north .ne. 1) .and. (ij .eq. (nlcj-1)) .and. (nbondj .le. 0)) then
899                        com_north = 1
900                     endif
901                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
902                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
903                  ENDIF
904                  ! check if point has to be received from a neighbour
905                  IF(nbondi .eq. 0) THEN
906                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   &
907                       & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   &
908                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
909                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
910                       if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then
911                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
912                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
913                            com_south = 1
914                          elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
915                            com_north = 1
916                          endif
917                          com_west_b = 1
918                       endif
919                     ENDIF
920                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND.   &
921                       & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND.   &
922                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
923                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
924                       if((com_east_b .ne. 1) .and. (ii .eq. 2)) then
925                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
926                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
927                            com_south = 1
928                          elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
929                            com_north = 1
930                          endif
931                          com_east_b = 1
932                       endif
933                     ENDIF
934                  ELSEIF(nbondi .eq. 1) THEN
935                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   &
936                       & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   &
937                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
938                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
939                       if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then
940                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
941                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
942                            com_south = 1
943                          elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
944                            com_north = 1
945                          endif
946                          com_west_b = 1
947                       endif
948                     ENDIF
949                  ELSEIF(nbondi .eq. -1) THEN
950                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND.   &
951                       & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND.   &
952                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
953                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
954                       if((com_east_b .ne. 1) .and. (ii .eq. 2)) then
955                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
956                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
957                            com_south = 1
958                          elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
959                            com_north = 1
960                          endif
961                          com_east_b = 1
962                       endif
963                     ENDIF
964                  ENDIF
965                  IF(nbondj .eq. 0) THEN
966                     IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1  &
967                       & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
968                       & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
969                       com_north_b = 1 
970                     ENDIF
971                     IF(com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1  &
972                       &.OR. nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
973                       & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
974                       com_south_b = 1 
975                     ENDIF
976                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND.   &
977                       & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND.   &
978                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
979                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
980                       if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then
981                          com_south_b = 1
982                       endif
983                     ENDIF
984                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND.   &
985                       & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND.   &
986                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
987                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
988                       if((com_north_b .ne. 1) .and. (ij .eq. 2)) then
989                          com_north_b = 1
990                       endif
991                     ENDIF
992                  ELSEIF(nbondj .eq. 1) THEN
993                     IF( com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. &
994                       & nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
995                       & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
996                       com_south_b = 1 
997                     ENDIF
998                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND.   &
999                       & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND.   &
1000                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
1001                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
1002                       if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then
1003                          com_south_b = 1
1004                       endif
1005                     ENDIF
1006                  ELSEIF(nbondj .eq. -1) THEN
1007                     IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1  &
1008                       & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
1009                       & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
1010                       com_north_b = 1 
1011                     ENDIF
1012                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND.   &
1013                       & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND.   &
1014                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
1015                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
1016                       if((com_north_b .ne. 1) .and. (ij .eq. 2)) then
1017                          com_north_b = 1
1018                       endif
1019                     ENDIF
1020                  ENDIF
1021               ENDDO
1022            ENDDO
1023         ENDDO 
1024
1025         ! definition of the i- and j- direction local boundaries arrays
1026         ! used for sending the boudaries
1027         IF((com_east .eq. 1) .and. (com_west .eq. 1)) THEN
1028            nbondi_bdy(ib_bdy) = 0
1029         ELSEIF ((com_east .eq. 1) .and. (com_west .eq. 0)) THEN
1030            nbondi_bdy(ib_bdy) = -1
1031         ELSEIF ((com_east .eq. 0) .and. (com_west .eq. 1)) THEN
1032            nbondi_bdy(ib_bdy) = 1
1033         ENDIF
1034
1035         IF((com_north .eq. 1) .and. (com_south .eq. 1)) THEN
1036            nbondj_bdy(ib_bdy) = 0
1037         ELSEIF ((com_north .eq. 1) .and. (com_south .eq. 0)) THEN
1038            nbondj_bdy(ib_bdy) = -1
1039         ELSEIF ((com_north .eq. 0) .and. (com_south .eq. 1)) THEN
1040            nbondj_bdy(ib_bdy) = 1
1041         ENDIF
1042
1043         ! definition of the i- and j- direction local boundaries arrays
1044         ! used for receiving the boudaries
1045         IF((com_east_b .eq. 1) .and. (com_west_b .eq. 1)) THEN
1046            nbondi_bdy_b(ib_bdy) = 0
1047         ELSEIF ((com_east_b .eq. 1) .and. (com_west_b .eq. 0)) THEN
1048            nbondi_bdy_b(ib_bdy) = -1
1049         ELSEIF ((com_east_b .eq. 0) .and. (com_west_b .eq. 1)) THEN
1050            nbondi_bdy_b(ib_bdy) = 1
1051         ENDIF
1052
1053         IF((com_north_b .eq. 1) .and. (com_south_b .eq. 1)) THEN
1054            nbondj_bdy_b(ib_bdy) = 0
1055         ELSEIF ((com_north_b .eq. 1) .and. (com_south_b .eq. 0)) THEN
1056            nbondj_bdy_b(ib_bdy) = -1
1057         ELSEIF ((com_north_b .eq. 0) .and. (com_south_b .eq. 1)) THEN
1058            nbondj_bdy_b(ib_bdy) = 1
1059         ENDIF
1060
1061         ! Compute rim weights for FRS scheme
1062         ! ----------------------------------
1063         DO igrd = 1, jpbgrd
1064            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
1065               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
1066               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
1067!               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.  ! quadratic
1068!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy))       ! linear
1069            END DO
1070         END DO 
1071
1072         ! Compute damping coefficients
1073         ! ----------------------------
1074         DO igrd = 1, jpbgrd
1075            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
1076               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
1077               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 
1078               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic
1079               idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) & 
1080               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic
1081            END DO
1082         END DO
1083
1084      ENDDO
1085
1086      ! ------------------------------------------------------
1087      ! Initialise masks and find normal/tangential directions
1088      ! ------------------------------------------------------
1089
1090      ! Read global 2D mask at T-points: bdytmask
1091      ! -----------------------------------------
1092      ! bdytmask = 1  on the computational domain AND on open boundaries
1093      !          = 0  elsewhere   
1094 
1095      IF( ln_mask_file ) THEN
1096         CALL iom_open( cn_mask_file, inum )
1097         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
1098         CALL iom_close( inum )
1099
1100         ! Derive mask on U and V grid from mask on T grid
1101         bdyumask(:,:) = 0.e0
1102         bdyvmask(:,:) = 0.e0
1103         DO ij=1, jpjm1
1104            DO ii=1, jpim1
1105               bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
1106               bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1) 
1107            END DO
1108         END DO
1109         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
1110
1111
1112         ! Mask corrections
1113         ! ----------------
1114         DO ik = 1, jpkm1
1115            DO ij = 1, jpj
1116               DO ii = 1, jpi
1117                  tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
1118                  umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
1119                  vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
1120                  bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij)
1121               END DO     
1122            END DO
1123         END DO
1124
1125         DO ik = 1, jpkm1
1126            DO ij = 2, jpjm1
1127               DO ii = 2, jpim1
1128                  fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   &
1129                     &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
1130               END DO     
1131            END DO
1132         END DO
1133
1134         tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)
1135
1136      ENDIF ! ln_mask_file=.TRUE.
1137     
1138      bdytmask(:,:) = tmask(:,:,1)
1139      IF( .not. ln_mask_file ) THEN
1140         ! If .not. ln_mask_file then we need to derive mask on U and V grid
1141         ! from mask on T grid here.
1142         bdyumask(:,:) = 0.e0
1143         bdyvmask(:,:) = 0.e0
1144         DO ij=1, jpjm1
1145            DO ii=1, jpim1
1146               bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
1147               bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1) 
1148            END DO
1149         END DO
1150         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
1151      ENDIF
1152
1153      ! bdy masks and bmask are now set to zero on boundary points:
1154      igrd = 1       ! In the free surface case, bmask is at T-points
1155      DO ib_bdy = 1, nb_bdy
1156        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
1157          bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1158        ENDDO
1159      ENDDO
1160      !
1161      igrd = 1
1162      DO ib_bdy = 1, nb_bdy
1163        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
1164          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1165        ENDDO
1166      ENDDO
1167      !
1168      igrd = 2
1169      DO ib_bdy = 1, nb_bdy
1170        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1171          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1172        ENDDO
1173      ENDDO
1174      !
1175      igrd = 3
1176      DO ib_bdy = 1, nb_bdy
1177        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1178          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1179        ENDDO
1180      ENDDO
1181
1182      ! For the flagu/flagv calculation below we require a version of fmask without
1183      ! the land boundary condition (shlat) included:
1184      CALL wrk_alloc(jpi,jpj,zfmask) 
1185      DO ij = 2, jpjm1
1186         DO ii = 2, jpim1
1187            zfmask(ii,ij) = tmask(ii,ij  ,1) * tmask(ii+1,ij  ,1)   &
1188           &              * tmask(ii,ij+1,1) * tmask(ii+1,ij+1,1)
1189         END DO     
1190      END DO
1191
1192      ! Lateral boundary conditions
1193      CALL lbc_lnk( zfmask       , 'F', 1. )
1194      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
1195      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
1196
1197      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
1198
1199         idx_bdy(ib_bdy)%flagu(:,:) = 0.e0
1200         idx_bdy(ib_bdy)%flagv(:,:) = 0.e0
1201         icount = 0 
1202
1203         ! Calculate relationship of U direction to the local orientation of the boundary
1204         ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
1205         ! flagu =  0 : u is tangential
1206         ! flagu =  1 : u is normal to the boundary and is direction is inward
1207 
1208         DO igrd = 1,jpbgrd 
1209            SELECT CASE( igrd )
1210               CASE( 1 )
1211                  pmask => umask(:,:,1)
1212                  i_offset = 0
1213               CASE( 2 ) 
1214                  pmask => bdytmask
1215                  i_offset = 1
1216               CASE( 3 ) 
1217                  pmask => zfmask(:,:)
1218                  i_offset = 0
1219            END SELECT
1220            icount = 0
1221            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
1222               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1223               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1224               zefl = pmask(nbi+i_offset-1,nbj)
1225               zwfl = pmask(nbi+i_offset,nbj)
1226               ! This error check only works if you are using the bdyXmask arrays
1227               IF( i_offset == 1 .and. zefl + zwfl == 2 ) THEN
1228                  icount = icount + 1
1229                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj)
1230               ELSE
1231                  idx_bdy(ib_bdy)%flagu(ib,igrd) = -zefl + zwfl
1232               ENDIF
1233            END DO
1234            IF( icount /= 0 ) THEN
1235               IF(lwp) WRITE(numout,*)
1236               IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   &
1237                  ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1238               IF(lwp) WRITE(numout,*) ' ========== '
1239               IF(lwp) WRITE(numout,*)
1240               nstop = nstop + 1
1241            ENDIF
1242         END DO
1243
1244         ! Calculate relationship of V direction to the local orientation of the boundary
1245         ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
1246         ! flagv =  0 : v is tangential
1247         ! flagv =  1 : v is normal to the boundary and is direction is inward
1248
1249         DO igrd = 1,jpbgrd 
1250            SELECT CASE( igrd )
1251               CASE( 1 )
1252                  pmask => vmask(:,:,1)
1253                  j_offset = 0
1254               CASE( 2 )
1255                  pmask => zfmask(:,:)
1256                  j_offset = 0
1257               CASE( 3 )
1258                  pmask => bdytmask
1259                  j_offset = 1
1260            END SELECT
1261            icount = 0
1262            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
1263               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1264               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1265               znfl = pmask(nbi,nbj+j_offset-1  )
1266               zsfl = pmask(nbi,nbj+j_offset)
1267               ! This error check only works if you are using the bdyXmask arrays
1268               IF( j_offset == 1 .and. znfl + zsfl == 2 ) THEN
1269                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj)
1270                  icount = icount + 1
1271               ELSE
1272                  idx_bdy(ib_bdy)%flagv(ib,igrd) = -znfl + zsfl
1273               END IF
1274            END DO
1275            IF( icount /= 0 ) THEN
1276               IF(lwp) WRITE(numout,*)
1277               IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   &
1278                  ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1279               IF(lwp) WRITE(numout,*) ' ========== '
1280               IF(lwp) WRITE(numout,*)
1281               nstop = nstop + 1
1282            ENDIF
1283         END DO
1284
1285      END DO
1286
1287      ! Compute total lateral surface for volume correction:
1288      ! ----------------------------------------------------
1289      ! JC: this must be done at each time step with key_vvl
1290      bdysurftot = 0.e0 
1291      IF( ln_vol ) THEN 
1292         igrd = 2      ! Lateral surface at U-points
1293         DO ib_bdy = 1, nb_bdy
1294            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1295               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1296               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1297               flagu => idx_bdy(ib_bdy)%flagu(ib,igrd)
1298               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           &
1299                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
1300                  &                    * tmask_i(nbi  , nbj)                           &
1301                  &                    * tmask_i(nbi+1, nbj)                   
1302            ENDDO
1303         ENDDO
1304
1305         igrd=3 ! Add lateral surface at V-points
1306         DO ib_bdy = 1, nb_bdy
1307            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1308               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1309               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1310               flagv => idx_bdy(ib_bdy)%flagv(ib,igrd)
1311               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           &
1312                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
1313                  &                    * tmask_i(nbi, nbj  )                           &
1314                  &                    * tmask_i(nbi, nbj+1)
1315            ENDDO
1316         ENDDO
1317         !
1318         IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain
1319      END IF   
1320      !
1321      ! Tidy up
1322      !--------
1323      IF (nb_bdy>0) THEN
1324         DEALLOCATE(nbidta, nbjdta, nbrdta)
1325      ENDIF
1326
1327      CALL wrk_dealloc(jpi,jpj,zfmask) 
1328
1329      IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
1330
1331   END SUBROUTINE bdy_init
1332
1333   SUBROUTINE bdy_ctl_seg
1334      !!----------------------------------------------------------------------
1335      !!                 ***  ROUTINE bdy_ctl_seg  ***
1336      !!
1337      !! ** Purpose :   Check straight open boundary segments location
1338      !!
1339      !! ** Method  :   - Look for open boundary corners
1340      !!                - Check that segments start or end on land
1341      !!----------------------------------------------------------------------
1342      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest 
1343      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns 
1344      REAL(wp), DIMENSION(2) ::   ztestmask
1345      !!----------------------------------------------------------------------
1346      !
1347      IF (lwp) WRITE(numout,*) ' '
1348      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1349      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1350      !
1351      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege
1352      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw
1353      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn
1354      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs
1355      ! 1. Check bounds
1356      !----------------
1357      DO ib = 1, nbdysegn
1358         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
1359         IF ((jpjnob(ib).ge.jpjglo-1).or.& 
1360            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1361         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1362         IF (jpindt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1363         IF (jpinft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
1364      END DO
1365      !
1366      DO ib = 1, nbdysegs
1367         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
1368         IF ((jpjsob(ib).ge.jpjglo-1).or.& 
1369            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1370         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1371         IF (jpisdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1372         IF (jpisft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
1373      END DO
1374      !
1375      DO ib = 1, nbdysege
1376         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib)
1377         IF ((jpieob(ib).ge.jpiglo-1).or.& 
1378            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1379         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1380         IF (jpjedt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1381         IF (jpjeft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
1382      END DO
1383      !
1384      DO ib = 1, nbdysegw
1385         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib)
1386         IF ((jpiwob(ib).ge.jpiglo-1).or.& 
1387            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1388         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1389         IF (jpjwdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1390         IF (jpjwft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
1391      ENDDO
1392      !
1393      !     
1394      ! 2. Look for segment crossings
1395      !------------------------------
1396      IF (lwp) WRITE(numout,*) '**Look for segments corners  :'
1397      !
1398      itest = 0 ! corner number
1399      !
1400      ! flag to detect if start or end of open boundary belongs to a corner
1401      ! if not (=0), it must be on land.
1402      ! if a corner is detected, save bdy package number for further tests
1403      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1404      ! South/West crossings
1405      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1406         DO ib1 = 1, nbdysegw       
1407            DO ib2 = 1, nbdysegs
1408               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1409                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1410                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1411                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1412                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN 
1413                     ! We have a possible South-West corner                     
1414!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1415!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1416                     icornw(ib1,1) = npckgs(ib2)
1417                     icorns(ib2,1) = npckgw(ib1)
1418                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1419                     IF(lwp) WRITE(numout,*)
1420                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1421                     &                                     jpisft(ib2), jpjwft(ib1)
1422                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1423                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), & 
1424                     &                                                    ' and South segment: ',npckgs(ib2)
1425                     IF(lwp) WRITE(numout,*)
1426                     nstop = nstop + 1
1427                  ELSE
1428                     IF(lwp) WRITE(numout,*)
1429                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices'
1430                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1) , &
1431                     &                                                    ' and South segment: ',npckgs(ib2)
1432                     IF(lwp) WRITE(numout,*)
1433                     nstop = nstop+1
1434                  END IF
1435               END IF
1436            END DO
1437         END DO
1438      END IF
1439      !
1440      ! South/East crossings
1441      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1442         DO ib1 = 1, nbdysege
1443            DO ib2 = 1, nbdysegs
1444               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1445                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1446                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. &
1447                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN
1448                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1449                     ! We have a possible South-East corner
1450!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1451!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1452                     icorne(ib1,1) = npckgs(ib2)
1453                     icorns(ib2,2) = npckge(ib1)
1454                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1455                     IF(lwp) WRITE(numout,*)
1456                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1457                     &                                     jpisdt(ib2), jpjeft(ib1)
1458                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1459                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), &
1460                     &                                                    ' and South segment: ',npckgs(ib2)
1461                     IF(lwp) WRITE(numout,*)
1462                     nstop = nstop + 1
1463                  ELSE
1464                     IF(lwp) WRITE(numout,*)
1465                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices'
1466                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
1467                     &                                                    ' and South segment: ',npckgs(ib2)
1468                     IF(lwp) WRITE(numout,*)
1469                     nstop = nstop + 1
1470                  END IF
1471               END IF
1472            END DO
1473         END DO
1474      END IF
1475      !
1476      ! North/West crossings
1477      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1478         DO ib1 = 1, nbdysegw       
1479            DO ib2 = 1, nbdysegn
1480               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. &
1481                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. &
1482                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1483                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1484                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1485                     ! We have a possible North-West corner
1486!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1487!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1488                     icornw(ib1,2) = npckgn(ib2)
1489                     icornn(ib2,1) = npckgw(ib1)
1490                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1491                     IF(lwp) WRITE(numout,*)
1492                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1493                     &                                     jpinft(ib2), jpjwdt(ib1)
1494                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1495                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), &
1496                     &                                                    ' and North segment: ',npckgn(ib2)
1497                     IF(lwp) WRITE(numout,*)
1498                     nstop = nstop + 1
1499                  ELSE
1500                     IF(lwp) WRITE(numout,*)
1501                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices'
1502                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1), &
1503                     &                                                    ' and North segment: ',npckgn(ib2)
1504                     IF(lwp) WRITE(numout,*)
1505                     nstop = nstop + 1
1506                  END IF
1507               END IF
1508            END DO
1509         END DO
1510      END IF
1511      !
1512      ! North/East crossings
1513      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1514         DO ib1 = 1, nbdysege       
1515            DO ib2 = 1, nbdysegn
1516               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1517                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1518                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1519                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1520                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1521                     ! We have a possible North-East corner
1522!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1523!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1524                     icorne(ib1,2) = npckgn(ib2)
1525                     icornn(ib2,2) = npckge(ib1)
1526                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1527                     IF(lwp) WRITE(numout,*)
1528                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1529                     &                                     jpindt(ib2), jpjedt(ib1)
1530                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1531                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), &
1532                     &                                                    ' and North segment: ',npckgn(ib2)
1533                     IF(lwp) WRITE(numout,*)
1534                     nstop = nstop + 1
1535                  ELSE
1536                     IF(lwp) WRITE(numout,*)
1537                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices'
1538                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
1539                     &                                                    ' and North segment: ',npckgn(ib2)
1540                     IF(lwp) WRITE(numout,*)
1541                     nstop = nstop + 1
1542                  END IF
1543               END IF
1544            END DO
1545         END DO
1546      END IF
1547      !
1548      ! 3. Check if segment extremities are on land
1549      !--------------------------------------------
1550      !
1551      ! West segments
1552      DO ib = 1, nbdysegw
1553         ! get mask at boundary extremities:
1554         ztestmask(1:2)=0.
1555         DO ji = 1, jpi
1556            DO jj = 1, jpj             
1557              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1558               &  ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1559              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1560               &  ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1561            END DO
1562         END DO
1563         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1564
1565         IF (ztestmask(1)==1) THEN
1566            IF (icornw(ib,1)==0) THEN
1567               IF(lwp) WRITE(numout,*)
1568               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
1569               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1570               IF(lwp) WRITE(numout,*)
1571               nstop = nstop + 1
1572            ELSE
1573               ! This is a corner
1574               WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1575               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1576               itest=itest+1
1577            ENDIF
1578         ENDIF
1579         IF (ztestmask(2)==1) THEN
1580            IF (icornw(ib,2)==0) THEN
1581               IF(lwp) WRITE(numout,*)
1582               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
1583               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1584               IF(lwp) WRITE(numout,*)
1585               nstop = nstop + 1
1586            ELSE
1587               ! This is a corner
1588               WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1589               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1590               itest=itest+1
1591            ENDIF
1592         ENDIF
1593      END DO
1594      !
1595      ! East segments
1596      DO ib = 1, nbdysege
1597         ! get mask at boundary extremities:
1598         ztestmask(1:2)=0.
1599         DO ji = 1, jpi
1600            DO jj = 1, jpj             
1601              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1602               &  ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
1603              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1604               &  ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1605            END DO
1606         END DO
1607         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1608
1609         IF (ztestmask(1)==1) THEN
1610            IF (icorne(ib,1)==0) THEN
1611               IF(lwp) WRITE(numout,*)
1612               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
1613               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1614               IF(lwp) WRITE(numout,*)
1615               nstop = nstop + 1 
1616            ELSE
1617               ! This is a corner
1618               WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1619               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1620               itest=itest+1
1621            ENDIF
1622         ENDIF
1623         IF (ztestmask(2)==1) THEN
1624            IF (icorne(ib,2)==0) THEN
1625               IF(lwp) WRITE(numout,*)
1626               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
1627               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1628               IF(lwp) WRITE(numout,*)
1629               nstop = nstop + 1
1630            ELSE
1631               ! This is a corner
1632               WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1633               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1634               itest=itest+1
1635            ENDIF
1636         ENDIF
1637      END DO
1638      !
1639      ! South segments
1640      DO ib = 1, nbdysegs
1641         ! get mask at boundary extremities:
1642         ztestmask(1:2)=0.
1643         DO ji = 1, jpi
1644            DO jj = 1, jpj             
1645              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1646               &  ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1647              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1648               &  ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1649            END DO
1650         END DO
1651         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1652
1653         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1654            IF(lwp) WRITE(numout,*)
1655            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
1656            IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1657            IF(lwp) WRITE(numout,*)
1658            nstop = nstop + 1
1659         ENDIF
1660         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1661            IF(lwp) WRITE(numout,*)
1662            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
1663            IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1664            IF(lwp) WRITE(numout,*)
1665            nstop = nstop + 1
1666         ENDIF
1667      END DO
1668      !
1669      ! North segments
1670      DO ib = 1, nbdysegn
1671         ! get mask at boundary extremities:
1672         ztestmask(1:2)=0.
1673         DO ji = 1, jpi
1674            DO jj = 1, jpj             
1675              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1676               &  ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
1677              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1678               &  ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1679            END DO
1680         END DO
1681         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1682
1683         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1684            IF(lwp) WRITE(numout,*)
1685            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
1686            IF(lwp) WRITE(numout,*) ' ==========  does not start on land'                                                 
1687            IF(lwp) WRITE(numout,*)
1688            nstop = nstop + 1
1689         ENDIF
1690         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1691            IF(lwp) WRITE(numout,*)
1692            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
1693            IF(lwp) WRITE(numout,*) ' ==========  does not end on land'                                                 
1694            IF(lwp) WRITE(numout,*)
1695            nstop = nstop + 1
1696         ENDIF
1697      END DO
1698      !
1699      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1700      !
1701      ! Other tests TBD:
1702      ! segments completly on land
1703      ! optimized open boundary array length according to landmask
1704      ! Nudging layers that overlap with interior domain
1705      !
1706   END SUBROUTINE bdy_ctl_seg
1707
1708   SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1709      !!----------------------------------------------------------------------
1710      !!                 ***  ROUTINE bdy_ctl_corn  ***
1711      !!
1712      !! ** Purpose :   Check numerical schemes consistency between
1713      !!                segments having a common corner
1714      !!
1715      !! ** Method  :   
1716      !!----------------------------------------------------------------------
1717      INTEGER, INTENT(in)  ::   ib1, ib2
1718      INTEGER :: itest
1719      !!----------------------------------------------------------------------
1720      itest = 0
1721
1722      IF (cn_dyn2d(ib1)/=cn_dyn2d(ib2)) itest = itest + 1
1723      IF (cn_dyn3d(ib1)/=cn_dyn3d(ib2)) itest = itest + 1
1724      IF (cn_tra(ib1)/=cn_tra(ib2)) itest = itest + 1
1725      !
1726      IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1
1727      IF (nn_dyn3d_dta(ib1)/=nn_dyn3d_dta(ib2)) itest = itest + 1
1728      IF (nn_tra_dta(ib1)/=nn_tra_dta(ib2)) itest = itest + 1
1729      !
1730      IF (nn_rimwidth(ib1)/=nn_rimwidth(ib2)) itest = itest + 1   
1731      !
1732      IF ( itest>0 ) THEN
1733         IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2
1734         IF(lwp) WRITE(numout,*) ' ==========  have different open bdy schemes'                                                 
1735         IF(lwp) WRITE(numout,*)
1736         nstop = nstop + 1
1737      ENDIF
1738      !
1739   END SUBROUTINE bdy_ctl_corn
1740
1741#else
1742   !!---------------------------------------------------------------------------------
1743   !!   Dummy module                                   NO open boundaries
1744   !!---------------------------------------------------------------------------------
1745CONTAINS
1746   SUBROUTINE bdy_init      ! Dummy routine
1747   END SUBROUTINE bdy_init
1748#endif
1749
1750   !!=================================================================================
1751END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.