source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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