source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 18 months ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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