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/UKMO/dev_r6393_CO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/UKMO/dev_r6393_CO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 7019

Last change on this file since 7019 was 7019, checked in by deazer, 7 years ago

Cleared svn keywords

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