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

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

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 4694

Last change on this file since 4694 was 4694, checked in by jamesharle, 10 years ago

Update of fldread to handle depth information in BDY files and addition of an interpolation routine. Updated BDY code to handle T/S BDY interpolation on the fly. Conservative remapping of U/V still to be coded. Not compiled or test yet.

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