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 NEMO/trunk/src/OCE/BDY – NEMO

source: NEMO/trunk/src/OCE/BDY/bdyini.F90 @ 15363

Last change on this file since 15363 was 15363, checked in by smasson, 12 months ago

trunk: continuing with ticket #2731

  • Property svn:keywords set to Id
File size: 107.8 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 sbc_oce , ONLY: nn_ice
22   USE bdy_oce        ! unstructured open boundary conditions
23   USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine)
24   USE bdytides       ! open boundary cond. setting   (bdytide_init routine)
25   USE tide_mod, ONLY: ln_tide ! tidal forcing
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
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   bdy_init    ! routine called in nemo_init
37   PUBLIC   find_neib   ! routine called in bdy_nmn
38
39   INTEGER, PARAMETER ::   jp_nseg = 100   !
40   ! Straight open boundary segment parameters:
41   INTEGER  ::   nbdysege, nbdysegw, nbdysegn, nbdysegs 
42   INTEGER, DIMENSION(jp_nseg) ::   jpieob, jpjedt, jpjeft, npckge   !
43   INTEGER, DIMENSION(jp_nseg) ::   jpiwob, jpjwdt, jpjwft, npckgw   !
44   INTEGER, DIMENSION(jp_nseg) ::   jpjnob, jpindt, jpinft, npckgn   !
45   INTEGER, DIMENSION(jp_nseg) ::   jpjsob, jpisdt, jpisft, npckgs   !
46   
47   !! * Substitutions
48#  include "do_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE bdy_init
57      !!----------------------------------------------------------------------
58      !!                 ***  ROUTINE bdy_init  ***
59      !!
60      !! ** Purpose :   Initialization of the dynamics and tracer fields with
61      !!              unstructured open boundaries.
62      !!
63      !! ** Method  :   Read initialization arrays (mask, indices) to identify
64      !!              an unstructured open boundary
65      !!
66      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
67      !!----------------------------------------------------------------------
68      NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,         &
69         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
70         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &
71         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
72         &             cn_ice, nn_ice_dta,                                     &
73         &             ln_vol, nn_volctl, nn_rimwidth
74         !
75      INTEGER  ::   ios                 ! Local integer output status for namelist read
76      !!----------------------------------------------------------------------
77
78      ! ------------------------
79      ! Read namelist parameters
80      ! ------------------------
81      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901)
82901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist' )
83      ! make sur that all elements of the namelist variables have a default definition from namelist_ref
84      ln_coords_file (2:jp_bdy) = ln_coords_file (1)
85      cn_coords_file (2:jp_bdy) = cn_coords_file (1)
86      cn_dyn2d       (2:jp_bdy) = cn_dyn2d       (1)
87      nn_dyn2d_dta   (2:jp_bdy) = nn_dyn2d_dta   (1)
88      cn_dyn3d       (2:jp_bdy) = cn_dyn3d       (1)
89      nn_dyn3d_dta   (2:jp_bdy) = nn_dyn3d_dta   (1)
90      cn_tra         (2:jp_bdy) = cn_tra         (1)
91      nn_tra_dta     (2:jp_bdy) = nn_tra_dta     (1)   
92      ln_tra_dmp     (2:jp_bdy) = ln_tra_dmp     (1)
93      ln_dyn3d_dmp   (2:jp_bdy) = ln_dyn3d_dmp   (1)
94      rn_time_dmp    (2:jp_bdy) = rn_time_dmp    (1)
95      rn_time_dmp_out(2:jp_bdy) = rn_time_dmp_out(1)
96      cn_ice         (2:jp_bdy) = cn_ice         (1)
97      nn_ice_dta     (2:jp_bdy) = nn_ice_dta     (1)
98      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
99902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist' )
100      IF(lwm) WRITE ( numond, nambdy )
101
102      IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE.   ! forced for Agrif children
103
104      IF( nb_bdy == 0 ) ln_bdy = .FALSE.
105     
106      ! -----------------------------------------
107      ! unstructured open boundaries use control
108      ! -----------------------------------------
109      IF ( ln_bdy ) THEN
110         IF(lwp) WRITE(numout,*)
111         IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
112         IF(lwp) WRITE(numout,*) '~~~~~~~~'
113         !
114         ! Open boundaries definition (arrays and masks)
115         CALL bdy_def
116         IF( ln_meshmask )   CALL bdy_meshwri()
117         !
118         ! Open boundaries initialisation of external data arrays
119         CALL bdy_dta_init
120         !
121         ! Open boundaries initialisation of tidal harmonic forcing
122         IF( ln_tide ) CALL bdytide_init
123         !
124      ELSE
125         IF(lwp) WRITE(numout,*)
126         IF(lwp) WRITE(numout,*) 'bdy_init : open boundaries not used (ln_bdy = F)'
127         IF(lwp) WRITE(numout,*) '~~~~~~~~'
128         !
129      ENDIF
130      !
131   END SUBROUTINE bdy_init
132
133
134   SUBROUTINE bdy_def
135      !!----------------------------------------------------------------------
136      !!                 ***  ROUTINE bdy_init  ***
137      !!         
138      !! ** Purpose :   Definition of unstructured open boundaries.
139      !!
140      !! ** Method  :   Read initialization arrays (mask, indices) to identify
141      !!              an unstructured open boundary
142      !!
143      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
144      !!----------------------------------------------------------------------     
145      INTEGER  ::   ji, jj                                 ! dummy loop indices
146      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, ir, iseg     ! dummy loop indices
147      INTEGER  ::   icount, icountr, icountr0, ibr_max     ! local integers
148      INTEGER  ::   ilen1                           !   -       -
149      INTEGER  ::   iiRst, iiRnd, iiSst, iiSnd, iiRcorn, iiSdiag, iiSsono
150      INTEGER  ::   ijRst, ijRnd, ijSst, ijSnd, ijRcorn, ijSdiag, ijSsono
151      INTEGER  ::   iiout, ijout, iioutdir, ijoutdir, icnt
152      INTEGER  ::   iRnei, iRdiag, iRsono
153      INTEGER  ::   iSnei, iSdiag, iSsono                  !   -       -
154      INTEGER  ::   iwe, ies, iso, ino, inum, id_dummy     !   -       -
155      INTEGER  ::   jpbdta                                 !   -       -
156      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       -
157      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3           !   -       -
158      INTEGER  ::   iibe, ijbe, iibi, ijbi                 !   -       -
159      INTEGER  ::   flagu, flagv                           ! short cuts
160      INTEGER  ::   nbdyind, nbdybeg, nbdyend
161      INTEGER              , DIMENSION(4)             ::   kdimsz
162      INTEGER              , DIMENSION(jpbgrd,jp_bdy) ::   nblendta          ! Length of index arrays
163      INTEGER,  ALLOCATABLE, DIMENSION(:,:,:)         ::   nbidta, nbjdta    ! Index arrays: i and j indices of bdy dta
164      INTEGER,  ALLOCATABLE, DIMENSION(:,:,:)         ::   nbrdta            ! Discrete distance from rim points
165      CHARACTER(LEN=1)     , DIMENSION(jpbgrd)        ::   cgrid
166      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zz_read                 ! work space for 2D global boundary data
167      REAL(wp), POINTER    , DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields
168      REAL(wp)             , DIMENSION(jpi,jpj) ::   zfmask   ! temporary fmask array excluding coastal boundary condition (shlat)
169      REAL(wp)             , DIMENSION(jpi,jpj) ::   ztmask, zumask, zvmask  ! temporary u/v mask array
170      REAL(wp)             , DIMENSION(jpi,jpj) ::   zzbdy
171      !!----------------------------------------------------------------------
172      !
173      cgrid = (/'t','u','v'/)
174
175      ! -----------------------------------------
176      ! Check and write out namelist parameters
177      ! -----------------------------------------
178     
179      IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy
180
181      DO ib_bdy = 1,nb_bdy
182
183         IF(lwp) THEN
184            WRITE(numout,*) ' ' 
185            WRITE(numout,*) '------ Open boundary data set ',ib_bdy,' ------' 
186            IF( ln_coords_file(ib_bdy) ) THEN
187               WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
188            ELSE
189               WRITE(numout,*) 'Boundary defined in namelist.'
190            ENDIF
191            WRITE(numout,*)
192         ENDIF
193
194         ! barotropic bdy
195         !----------------
196         IF(lwp) THEN
197            WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
198            SELECT CASE( cn_dyn2d(ib_bdy) )                 
199            CASE( 'none' )           ;   WRITE(numout,*) '      no open boundary condition'       
200            CASE( 'frs' )            ;   WRITE(numout,*) '      Flow Relaxation Scheme'
201            CASE( 'flather' )        ;   WRITE(numout,*) '      Flather radiation condition'
202            CASE( 'orlanski' )       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
203            CASE( 'orlanski_npo' )   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
204            CASE DEFAULT             ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' )
205            END SELECT
206         ENDIF
207
208         dta_bdy(ib_bdy)%lneed_ssh   = cn_dyn2d(ib_bdy) == 'flather'
209         dta_bdy(ib_bdy)%lneed_dyn2d = cn_dyn2d(ib_bdy) /= 'none'
210
211         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
212            SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !
213            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'       
214            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
215            CASE( 2 )      ;   WRITE(numout,*) '      tidal harmonic forcing taken from file'
216            CASE( 3 )      ;   WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files'
217            CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
218            END SELECT
219         ENDIF
220         IF ( dta_bdy(ib_bdy)%lneed_dyn2d .AND. nn_dyn2d_dta(ib_bdy) .GE. 2  .AND. .NOT.ln_tide ) THEN
221            CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' )
222         ENDIF
223         IF(lwp) WRITE(numout,*)
224
225         ! baroclinic bdy
226         !----------------
227         IF(lwp) THEN
228            WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
229            SELECT CASE( cn_dyn3d(ib_bdy) )                 
230            CASE('none')           ;   WRITE(numout,*) '      no open boundary condition'       
231            CASE('frs')            ;   WRITE(numout,*) '      Flow Relaxation Scheme'
232            CASE('specified')      ;   WRITE(numout,*) '      Specified value'
233            CASE('neumann')        ;   WRITE(numout,*) '      Neumann conditions'
234            CASE('zerograd')       ;   WRITE(numout,*) '      Zero gradient for baroclinic velocities'
235            CASE('zero')           ;   WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
236            CASE('orlanski')       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
237            CASE('orlanski_npo')   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
238            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn3d' )
239            END SELECT
240         ENDIF
241
242         dta_bdy(ib_bdy)%lneed_dyn3d = cn_dyn3d(ib_bdy) == 'frs'      .OR. cn_dyn3d(ib_bdy) == 'specified'   &
243            &                     .OR. cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo'
244
245         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn3d ) THEN
246            SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !
247            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'       
248            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
249            CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
250            END SELECT
251         END IF
252
253         IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
254            IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN
255               IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
256               ln_dyn3d_dmp(ib_bdy) = .false.
257            ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN
258               CALL ctl_stop( 'Use FRS OR relaxation' )
259            ELSE
260               IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone'
261               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
262               IF(rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
263               dta_bdy(ib_bdy)%lneed_dyn3d = .TRUE.
264            ENDIF
265         ELSE
266            IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities'
267         ENDIF
268         IF(lwp) WRITE(numout,*)
269
270         !    tra bdy
271         !----------------
272         IF(lwp) THEN
273            WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
274            SELECT CASE( cn_tra(ib_bdy) )                 
275            CASE('none')           ;   WRITE(numout,*) '      no open boundary condition'       
276            CASE('frs')            ;   WRITE(numout,*) '      Flow Relaxation Scheme'
277            CASE('specified')      ;   WRITE(numout,*) '      Specified value'
278            CASE('neumann')        ;   WRITE(numout,*) '      Neumann conditions'
279            CASE('runoff')         ;   WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity'
280            CASE('orlanski')       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
281            CASE('orlanski_npo')   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
282            CASE DEFAULT           ;   CALL ctl_stop( 'unrecognised value for cn_tra' )
283            END SELECT
284         ENDIF
285
286         dta_bdy(ib_bdy)%lneed_tra = cn_tra(ib_bdy) == 'frs'       .OR. cn_tra(ib_bdy) == 'specified'   &
287            &                   .OR. cn_tra(ib_bdy) == 'orlanski'  .OR. cn_tra(ib_bdy) == 'orlanski_npo' 
288
289         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_tra ) THEN
290            SELECT CASE( nn_tra_dta(ib_bdy) )                   !
291            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'       
292            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
293            CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
294            END SELECT
295         ENDIF
296
297         IF ( ln_tra_dmp(ib_bdy) ) THEN
298            IF ( cn_tra(ib_bdy) == 'none' ) THEN
299               IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
300               ln_tra_dmp(ib_bdy) = .false.
301            ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN
302               CALL ctl_stop( 'Use FRS OR relaxation' )
303            ELSE
304               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
305               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
306               IF(lwp) WRITE(numout,*) '      Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'
307               IF(lwp.AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
308               dta_bdy(ib_bdy)%lneed_tra = .TRUE.
309            ENDIF
310         ELSE
311            IF(lwp) WRITE(numout,*) '      NO T/S relaxation'
312         ENDIF
313         IF(lwp) WRITE(numout,*)
314
315#if defined key_si3
316         IF(lwp) THEN
317            WRITE(numout,*) 'Boundary conditions for sea ice:  '
318            SELECT CASE( cn_ice(ib_bdy) )                 
319            CASE('none')   ;   WRITE(numout,*) '      no open boundary condition'       
320            CASE('frs')    ;   WRITE(numout,*) '      Flow Relaxation Scheme'
321            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice' )
322            END SELECT
323         ENDIF
324
325         dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none'
326
327         IF( dta_bdy(ib_bdy)%lneed_ice .AND. nn_ice /= 2 ) THEN
328            WRITE(ctmp1,*) 'bdy number ', ib_bdy,', needs ice model but nn_ice = ', nn_ice
329            CALL ctl_stop( ctmp1 )
330         ENDIF
331
332         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN
333            SELECT CASE( nn_ice_dta(ib_bdy) )                   !
334            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'       
335            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
336            CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_dta must be 0 or 1' )
337            END SELECT
338         ENDIF
339#else
340         dta_bdy(ib_bdy)%lneed_ice = .FALSE.
341#endif
342         !
343         IF(lwp) WRITE(numout,*)
344         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
345         IF(lwp) WRITE(numout,*)
346         !
347      END DO   ! nb_bdy
348
349      IF( lwp ) THEN
350         IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
351            WRITE(numout,*) 'Volume correction applied at open boundaries'
352            WRITE(numout,*)
353            SELECT CASE ( nn_volctl )
354            CASE( 1 )      ;   WRITE(numout,*) '      The total volume will be constant'
355            CASE( 0 )      ;   WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
356            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
357            END SELECT
358            WRITE(numout,*)
359            !
360            ! sanity check if used with tides       
361            IF( ln_tide ) THEN
362               WRITE(numout,*) ' The total volume correction is not working with tides. '
363               WRITE(numout,*) ' Set ln_vol to .FALSE. '
364               WRITE(numout,*) ' or '
365               WRITE(numout,*) ' equilibriate your bdy input files '
366               CALL ctl_stop( 'The total volume correction is not working with tides.' )
367            END IF
368         ELSE
369            WRITE(numout,*) 'No volume correction applied at open boundaries'
370            WRITE(numout,*)
371         ENDIF
372      ENDIF
373
374      ! -------------------------------------------------
375      ! Initialise indices arrays for open boundaries
376      ! -------------------------------------------------
377
378      nblendta(:,:) = 0
379      nbdysege = 0
380      nbdysegw = 0
381      nbdysegn = 0
382      nbdysegs = 0
383
384      ! Define all boundaries
385      ! ---------------------
386      DO ib_bdy = 1, nb_bdy
387         !
388         IF( .NOT. ln_coords_file(ib_bdy) ) THEN     ! build bdy coordinates with segments defined in namelist
389
390            CALL bdy_read_seg( ib_bdy, nblendta(:,ib_bdy) )
391
392         ELSE                                        ! Read size of arrays in boundary coordinates file.
393           
394            CALL iom_open( cn_coords_file(ib_bdy), inum )
395            DO igrd = 1, jpbgrd
396               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 
397               nblendta(igrd,ib_bdy) = MAXVAL(kdimsz)
398            END DO
399            CALL iom_close( inum )
400         ENDIF
401         !
402      END DO ! ib_bdy
403
404      ! Now look for crossings in user (namelist) defined open boundary segments:
405      IF( nbdysege > 0 .OR. nbdysegw > 0 .OR. nbdysegn > 0 .OR. nbdysegs > 0)   CALL bdy_ctl_seg
406     
407      ! Allocate arrays
408      !---------------
409      jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
410      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), nbrdta(jpbdta, jpbgrd, nb_bdy) )
411      nbrdta(:,:,:) = 0   ! initialize nbrdta as it may not be completely defined for each bdy
412     
413      ! Calculate global boundary index arrays or read in from file
414      !------------------------------------------------------------               
415      ! 1. Read global index arrays from boundary coordinates file.
416      DO ib_bdy = 1, nb_bdy
417         !
418         IF( ln_coords_file(ib_bdy) ) THEN
419            !
420            ALLOCATE( zz_read( MAXVAL(nblendta), 1 ) )         
421            CALL iom_open( cn_coords_file(ib_bdy), inum )
422            !
423            DO igrd = 1, jpbgrd
424               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
425               DO ii = 1,nblendta(igrd,ib_bdy)
426                  nbidta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) + nn_hls
427               END DO
428               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
429               DO ii = 1,nblendta(igrd,ib_bdy)
430                  nbjdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) + nn_hls
431               END DO
432               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
433               DO ii = 1,nblendta(igrd,ib_bdy)
434                  nbrdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
435               END DO
436               !
437               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
438               IF(lwp) WRITE(numout,*)
439               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
440               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
441               IF (ibr_max < nn_rimwidth(ib_bdy))   &
442                  CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
443            END DO
444            !
445            CALL iom_close( inum )
446            DEALLOCATE( zz_read )
447            !
448         ENDIF
449         !
450      END DO
451
452      ! 2. Now fill indices corresponding to straight open boundary arrays:
453      CALL bdy_coords_seg( nbidta, nbjdta, nbrdta )
454
455      !  Deal with duplicated points
456      !-----------------------------
457      ! We assign negative indices to duplicated points (to remove them from bdy points to be updated)
458      ! if their distance to the bdy is greater than the other
459      ! If their distance are the same, just keep only one to avoid updating a point twice
460      DO igrd = 1, jpbgrd
461         DO ib_bdy1 = 1, nb_bdy
462            DO ib_bdy2 = 1, nb_bdy
463               IF (ib_bdy1/=ib_bdy2) THEN
464                  DO ib1 = 1, nblendta(igrd,ib_bdy1)
465                     DO ib2 = 1, nblendta(igrd,ib_bdy2)
466                        IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
467                           &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
468                           !                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &
469                           !                                                       &              nbidta(ib1, igrd, ib_bdy1),      &
470                           !                                                       &              nbjdta(ib2, igrd, ib_bdy2)
471                           ! keep only points with the lowest distance to boundary:
472                           IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN
473                              nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
474                              nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
475                           ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN
476                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
477                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
478                              ! Arbitrary choice if distances are the same:
479                           ELSE
480                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
481                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
482                           ENDIF
483                        END IF
484                     END DO
485                  END DO
486               ENDIF
487            END DO
488         END DO
489      END DO
490      !
491      ! Find lenght of boundaries and rim on local mpi domain
492      !------------------------------------------------------
493      !
494      iwe = mig(1)
495      ies = mig(jpi)
496      iso = mjg(1) 
497      ino = mjg(jpj) 
498      !
499      DO ib_bdy = 1, nb_bdy
500         DO igrd = 1, jpbgrd
501            icount   = 0   ! initialization of local bdy length
502            icountr  = 0   ! initialization of local rim 0 and rim 1 bdy length
503            icountr0 = 0   ! initialization of local rim 0 bdy length
504            idx_bdy(ib_bdy)%nblen(igrd)     = 0
505            idx_bdy(ib_bdy)%nblenrim(igrd)  = 0
506            idx_bdy(ib_bdy)%nblenrim0(igrd) = 0
507            DO ib = 1, nblendta(igrd,ib_bdy)
508               ! check that data is in correct order in file
509               IF( ib > 1 ) THEN
510                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ib-1,igrd,ib_bdy) ) THEN
511                     CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', &
512                        &        ' in order of distance from edge nbr A utility for re-ordering ', &
513                        &        ' boundary coordinates and data files exists in the TOOLS/OBC directory')
514                  ENDIF
515               ENDIF
516               ! check if point is in local domain
517               IF(  nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND.   &
518                  & nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino      ) THEN
519                  !
520                  icount = icount + 1
521                  IF( nbrdta(ib,igrd,ib_bdy) == 1 .OR. nbrdta(ib,igrd,ib_bdy) == 0 )   icountr = icountr + 1
522                  IF( nbrdta(ib,igrd,ib_bdy) == 0 )   icountr0 = icountr0 + 1
523               ENDIF
524            END DO
525            idx_bdy(ib_bdy)%nblen    (igrd) = icount   !: length of boundary data on each proc
526            idx_bdy(ib_bdy)%nblenrim (igrd) = icountr  !: length of rim 0 and rim 1 boundary data on each proc   
527            idx_bdy(ib_bdy)%nblenrim0(igrd) = icountr0 !: length of rim 0 boundary data on each proc     
528         END DO   ! igrd
529
530         ! Allocate index arrays for this boundary set
531         !--------------------------------------------
532         ilen1 = MAXVAL( idx_bdy(ib_bdy)%nblen(:) )
533         ALLOCATE( idx_bdy(ib_bdy)%nbi   (ilen1,jpbgrd) ,   &
534            &      idx_bdy(ib_bdy)%nbj   (ilen1,jpbgrd) ,   &
535            &      idx_bdy(ib_bdy)%nbr   (ilen1,jpbgrd) ,   &
536            &      idx_bdy(ib_bdy)%nbd   (ilen1,jpbgrd) ,   &
537            &      idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) ,   &
538            &      idx_bdy(ib_bdy)%ntreat(ilen1,jpbgrd) ,   &
539            &      idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) ,   &
540            &      idx_bdy(ib_bdy)%nbw   (ilen1,jpbgrd) ,   &
541            &      idx_bdy(ib_bdy)%flagu (ilen1,jpbgrd) ,   &
542            &      idx_bdy(ib_bdy)%flagv (ilen1,jpbgrd) )
543
544         ! Dispatch mapping indices and discrete distances on each processor
545         ! -----------------------------------------------------------------
546         DO igrd = 1, jpbgrd
547            icount  = 0
548            ! Outer loop on rimwidth to ensure outermost points come first in the local arrays.
549            DO ir = 0, nn_rimwidth(ib_bdy)
550               DO ib = 1, nblendta(igrd,ib_bdy)
551                  ! check if point is in local domain and equals ir
552                  IF(  nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND.   &
553                     & nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino .AND.   &
554                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
555                     !
556                     icount = icount  + 1
557                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy) - mig(1) + 1   ! global to local indexes
558                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy) - mjg(1) + 1   ! global to local indexes
559                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
560                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
561                  ENDIF
562               END DO
563            END DO
564         END DO   ! igrd
565
566      END DO   ! ib_bdy
567
568      ! Initialize array indicating communications in bdy
569      ! -------------------------------------------------
570      ALLOCATE( lsend_bdyolr(nb_bdy,jpbgrd,8,0:1), lrecv_bdyolr(nb_bdy,jpbgrd,8,0:1) )
571      lsend_bdyolr(:,:,:,:) = .false.
572      lrecv_bdyolr(:,:,:,:) = .false. 
573
574      DO ib_bdy = 1, nb_bdy
575         DO igrd = 1, jpbgrd
576            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! only the rim triggers communications, see bdy routines
577               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
578               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
579               IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN   ;   ir = 0
580               ELSE                                                 ;   ir = 1
581               END IF
582               !
583               ! check if point has to be sent     to   a neighbour
584               IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. ij >= Njs0 .AND. ij <= Nje0         ) THEN   ! we inner side
585                  IF( mpiSnei(nn_hls,jpwe) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpwe,ir) = .TRUE.
586               ENDIF
587               IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. ij >= Njs0 .AND. ij <= Nje0         ) THEN   ! ea inner side
588                  IF( mpiSnei(nn_hls,jpea) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpea,ir) = .TRUE.
589               ENDIF
590               IF( ii >= Nis0 .AND. ii <= Nie0         .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! so inner side
591                  IF( mpiSnei(nn_hls,jpso) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
592               ENDIF
593               IF( ii  < Nis0                          .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! so side we-halo
594                  IF( mpiSnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 )   lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
595               ENDIF
596               IF( ii  > Nie0                          .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! so side ea-halo
597                  IF( mpiSnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 )   lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
598               ENDIF
599               IF( ii >= Nis0 .AND. ii <= Nie0         .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! no inner side
600                  IF( mpiSnei(nn_hls,jpno) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
601               ENDIF
602               IF( ii  < Nis0                          .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! no side we-halo
603                  IF( mpiSnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 )   lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
604               ENDIF
605               IF( ii  > Nie0                          .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! no side ea-halo
606                  IF( mpiSnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 )   lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
607               ENDIF
608               IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! sw inner corner
609                  IF( mpiSnei(nn_hls,jpsw) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpsw,ir) = .TRUE.
610               ENDIF
611               IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! se inner corner
612                  IF( mpiSnei(nn_hls,jpse) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpse,ir) = .TRUE.
613               ENDIF
614               IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! nw inner corner
615                  IF( mpiSnei(nn_hls,jpnw) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpnw,ir) = .TRUE.
616               ENDIF
617               IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! ne inner corner
618                  IF( mpiSnei(nn_hls,jpne) > -1                    )   lsend_bdyolr(ib_bdy,igrd,jpne,ir) = .TRUE.
619               ENDIF
620               !
621               ! check if point has to be received from a neighbour
622               IF( ii  < Nis0                  .AND. ij >= Njs0 .AND. ij <= Nje0 ) THEN   ! we side
623                  IF( mpiRnei(nn_hls,jpwe) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpwe,ir) = .TRUE.
624               ENDIF
625               IF( ii  > Nie0                  .AND. ij >= Njs0 .AND. ij <= Nje0 ) THEN   ! ea side
626                  IF( mpiRnei(nn_hls,jpea) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpea,ir) = .TRUE.
627               ENDIF
628               IF( ii >= Nis0 .AND. ii <= Nie0 .AND. ij  < Njs0                  ) THEN   ! so side
629                  IF( mpiRnei(nn_hls,jpso) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
630               ENDIF
631               IF( ii >= Nis0 .AND. ii <= Nie0 .AND. ij  > Nje0                  ) THEN   ! no side
632                  IF( mpiRnei(nn_hls,jpno) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
633               ENDIF
634               IF( ii  < Nis0                  .AND. ij  < Njs0                  ) THEN   ! sw corner
635                  IF( mpiRnei(nn_hls,jpsw) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpsw,ir) = .TRUE.
636                  IF( mpiRnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 )   lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
637               ENDIF
638               IF( ii  > Nie0                  .AND. ij  < Njs0                  ) THEN   ! se corner
639                  IF( mpiRnei(nn_hls,jpse) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpse,ir) = .TRUE.
640                  IF( mpiRnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 )   lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
641               ENDIF
642               IF( ii  < Nis0                  .AND. ij  > Nje0                  ) THEN   ! nw corner
643                  IF( mpiRnei(nn_hls,jpnw) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpnw,ir) = .TRUE.
644                  IF( mpiRnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 )   lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
645               ENDIF
646               IF( ii  > Nie0                  .AND. ij  > Nje0                  ) THEN   ! ne corner
647                  IF( mpiRnei(nn_hls,jpne) > -1                    )   lrecv_bdyolr(ib_bdy,igrd,jpne,ir) = .TRUE.
648                  IF( mpiRnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 )   lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
649               ENDIF
650               !
651            END DO
652         END DO   !   igrd
653         
654         ! Comment out for debug
655!!$         DO ir = 0,1
656!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing,   &
657!!$               &                              lsend = lsend_bdyolr(ib_bdy,1,:,ir), lrecv = lrecv_bdyolr(ib_bdy,1,:,ir) )
658!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug olr T', ir ; CALL FLUSH(numout)
659!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing,   &
660!!$               &                              lsend = lsend_bdyolr(ib_bdy,2,:,ir), lrecv = lrecv_bdyolr(ib_bdy,2,:,ir) )
661!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug olr U', ir ; CALL FLUSH(numout)
662!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing,   &
663!!$               &                              lsend = lsend_bdyolr(ib_bdy,3,:,ir), lrecv = lrecv_bdyolr(ib_bdy,3,:,ir) )   
664!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug olr V', ir ; CALL FLUSH(numout)
665!!$         END DO
666         
667         ! Compute rim weights for FRS scheme
668         ! ----------------------------------
669         DO igrd = 1, jpbgrd
670            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
671               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same weights
672               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 )      ! tanh formulation
673               !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic
674               !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy))       ! linear
675            END DO
676         END DO
677
678         ! Compute damping coefficients
679         ! ----------------------------
680         DO igrd = 1, jpbgrd
681            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
682               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same damping coefficients
683               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 
684                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
685               idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) & 
686                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
687            END DO
688         END DO
689
690      END DO ! ib_bdy
691
692      ! ------------------------------------------------------
693      ! Initialise masks and find normal/tangential directions
694      ! ------------------------------------------------------
695
696      ! ------------------------------------------
697      ! handle rim0, do as if rim 1 was free ocean
698      ! ------------------------------------------
699
700      ztmask(:,:) = tmask(:,:,1)   ;   zumask(:,:) = umask(:,:,1)   ;   zvmask(:,:) = vmask(:,:,1)
701      ! For the flagu/flagv calculation below we require a version of fmask without
702      ! the land boundary condition (shlat) included:
703      DO_2D( 0, 0, 0, 0 )
704         zfmask(ji,jj) =  ztmask(ji,jj  ) * ztmask(ji+1,jj  )   &
705            &           * ztmask(ji,jj+1) * ztmask(ji+1,jj+1)
706      END_2D
707      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp )
708
709      ! Read global 2D mask at T-points: bdytmask
710      ! -----------------------------------------
711      ! bdytmask = 1  on the computational domain but not on open boundaries
712      !          = 0  elsewhere   
713
714      bdytmask(:,:) = ssmask(:,:)
715
716      ! Derive mask on U and V grid from mask on T grid
717      DO_2D( 0, 0, 0, 0 )
718            bdyumask(ji,jj) = bdytmask(ji,jj) * bdytmask(ji+1,jj  )
719            bdyvmask(ji,jj) = bdytmask(ji,jj) * bdytmask(ji  ,jj+1) 
720      END_2D
721      CALL lbc_lnk( 'bdyini', bdyumask, 'U', 1.0_wp , bdyvmask, 'V', 1.0_wp )   ! Lateral boundary cond.
722
723      ! bdy masks are now set to zero on rim 0 points:
724      DO ib_bdy = 1, nb_bdy
725         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0
726            bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
727         END DO
728         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2)   ! extent of rim 0
729            bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
730         END DO
731         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3)   ! extent of rim 0
732            bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
733         END DO
734      END DO
735
736      CALL bdy_rim_treat( zumask, zvmask, zfmask, .true. )   ! compute flagu, flagv, ntreat on rim 0
737
738      ! ------------------------------------
739      ! handle rim1, do as if rim 0 was land
740      ! ------------------------------------
741     
742      ! z[tuv]mask are now set to zero on rim 0 points:
743      DO ib_bdy = 1, nb_bdy
744         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0
745            ztmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
746         END DO
747         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2)   ! extent of rim 0
748            zumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
749         END DO
750         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3)   ! extent of rim 0
751            zvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
752         END DO
753      END DO
754
755      ! Recompute zfmask
756      DO_2D( 0, 0, 0, 0 )
757         zfmask(ji,jj) =  ztmask(ji,jj  ) * ztmask(ji+1,jj  )   &
758            &           * ztmask(ji,jj+1) * ztmask(ji+1,jj+1)
759      END_2D
760      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp )
761
762      ! bdy masks are now set to zero on rim1 points:
763      DO ib_bdy = 1, nb_bdy
764         DO ib = idx_bdy(ib_bdy)%nblenrim0(1) + 1,  idx_bdy(ib_bdy)%nblenrim(1)   ! extent of rim 1
765            bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
766         END DO
767         DO ib = idx_bdy(ib_bdy)%nblenrim0(2) + 1,  idx_bdy(ib_bdy)%nblenrim(2)   ! extent of rim 1
768            bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
769         END DO
770         DO ib = idx_bdy(ib_bdy)%nblenrim0(3) + 1,  idx_bdy(ib_bdy)%nblenrim(3)   ! extent of rim 1
771            bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
772         END DO
773      END DO
774
775      CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. )   ! compute flagu, flagv, ntreat on rim 1
776      !
777      ! Check which boundaries might need communication
778      ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,8,0:1), lrecv_bdyint(nb_bdy,jpbgrd,8,0:1) )
779      lsend_bdyint(:,:,:,:) = .false.
780      lrecv_bdyint(:,:,:,:) = .false. 
781      ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,8,0:1), lrecv_bdyext(nb_bdy,jpbgrd,8,0:1) )
782      lsend_bdyext(:,:,:,:) = .false.
783      lrecv_bdyext(:,:,:,:) = .false.
784      !
785      DO ib_bdy = 1, nb_bdy
786         DO igrd = 1, jpbgrd
787            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
788               IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE
789               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
790               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
791               ir = idx_bdy(ib_bdy)%nbr(ib,igrd)
792               flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd))
793               flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd))
794               iibe = ii - flagu   ! neighbouring point towards the exterior of the computational domain
795               ijbe = ij - flagv
796               iibi = ii + flagu   ! neighbouring point towards the interior of the computational domain
797               ijbi = ij + flagv
798               CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 )   ! free ocean neighbours
799               !
800               !  take care of the 4 sides
801               !
802               DO icnt = 1, 4
803                  SELECT CASE( icnt )
804                     !                                           ... _____
805                  CASE( 1 )   ! x: rim on rcvwe/sndea-side         o|  :
806                     !          o: potential neighbour(s)          o|x :
807                     !             outside of the MPI domain     ..o|__:__
808                     iRnei    = jpwe             ;   iSnei    = jpea
809                     iiRst    = 1                ;   ijRst    = 2               ! Rcv we-side starting point, excluding sw-corner
810                     iiRnd    = 1                ;   ijRnd    = jpj-1           ! Rcv we-side   ending point, excluding nw-corner
811                     iiSst    = jpi-2*nn_hls+1   ;   ijSst    = 2               ! Snd ea-side starting point, excluding se-corner
812                     iiSnd    = jpi-2*nn_hls+1   ;   ijSnd    = jpj-1           ! Snd ea-side   ending point, excluding ne-corner
813                     iioutdir = -1               ;   ijoutdir = -999            ! outside MPI domain: westward
814                     !                                           ______....
815                  CASE( 2 )   ! x: rim on rcvea/sndwe-side          :  |o
816                     !          o: potential neighbour(s)           : x|o
817                     !             outside of the MPI domain     ___:__|o..
818                     iRnei    = jpea             ;   iSnei    = jpwe
819                     iiRst    = jpi              ;   ijRst    = 2                ! Rcv ea-side starting point, excluding se-corner
820                     iiRnd    = jpi              ;   ijRnd    = jpj-1            ! Rcv ea-side   ending point, excluding ne-corner
821                     iiSst    = 2*nn_hls         ;   ijSst    = 2                ! Snd we-side starting point, excluding sw-corner
822                     iiSnd    = 2*nn_hls         ;   ijSnd    = jpj-1            ! Snd we-side   ending point, excluding nw-corner
823                     iioutdir = 1                ;   ijoutdir = -999             ! outside MPI domain: eastward
824                     !
825                  CASE( 3 )   ! x: rim on rcvso/sndno-side       |       |
826                     !          o: potential neighbour(s)        |¨¨¨¨¨¨¨|
827                     !             outside of the MPI domain     |___x___|
828                     !                                           : o o o :
829                     !                                           :       :
830                     iRnei    = jpso             ;   iSnei    = jpno
831                     iiRst    = 2                ;   ijRst    = 1                ! Rcv so-side starting point, excluding sw-corner
832                     iiRnd    = jpi-1            ;   ijRnd    = 1                ! Rcv so-side   ending point, excluding se-corner
833                     iiSst    = 2                ;   ijSst    = jpj-2*nn_hls+1   ! Snd no-side starting point, excluding nw-corner
834                     iiSnd    = jpi-1            ;   ijSnd    = jpj-2*nn_hls+1   ! Snd no-side   ending point, excluding ne-corner
835                     iioutdir = -999             ;   ijoutdir = -1               ! outside MPI domain: southward
836                     !                                           :       :
837                  CASE( 4 )   ! x: rim on rcvno/sndso-side       :_o_o_o_:
838                     !          o: potential neighbour(s)        |   x   |
839                     !             outside of the MPI domain     |       |
840                     !                                           |¨¨¨¨¨¨¨|
841                     iRnei    = jpno             ;   iSnei    = jpso
842                     iiRst    = 2                ;   ijRst    = jpj              ! Rcv no-side starting point, excluding nw-corner
843                     iiRnd    = jpi-1            ;   ijRnd    = jpj              ! Rcv no-side   ending point, excluding ne-corner
844                     iiSst    = 2                ;   ijSst    =     2*nn_hls     ! Snd so-side starting point, excluding sw-corner
845                     iiSnd    = jpi-1            ;   ijSnd    =     2*nn_hls     ! Snd so-side   ending point, excluding se-corner
846                     iioutdir = -999             ;   ijoutdir = 1                ! outside MPI domain: northward
847                  END SELECT
848                  !
849                  IF( ii >= iiRst .AND. ii <= iiRnd .AND. ij >= ijRst .AND. ij <= ijRnd ) THEN   ! rim point in recv side
850                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the MPI domain?
851                     ! take care of neighbourg(s) in the interior of the computational domain
852                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of the MPI domain
853                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN     ! -> I cannot compute it -> recv it
854                        IF( mpiRnei(nn_hls,iRnei) > -1 )   lrecv_bdyint(ib_bdy,igrd,iRnei,ir) = .TRUE.
855                     ENDIF
856                     ! take care of neighbourg in the exterior of the computational domain
857                     IF(  iibe==iiout .OR. ijbe==ijout ) THEN   ! Neib outside of the MPI domain -> I cannot compute it -> recv it
858                        IF( mpiRnei(nn_hls,iRnei) > -1 )   lrecv_bdyext(ib_bdy,igrd,iRnei,ir) = .TRUE.
859                     ENDIF
860                  ENDIF
861                 
862                  IF( ii >= iiSst .AND. ii <= iiSnd .AND. ij >= ijSst .AND. ij <= ijSnd ) THEN   ! rim point in send side
863                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the nei MPI domain?
864                     ! take care of neighbourg(s) in the interior of the computational domain
865                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of nei MPI domain
866                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN     ! -> nei cannot compute it
867                        IF( mpiSnei(nn_hls,iSnei) > -1 )   lsend_bdyint(ib_bdy,igrd,iSnei,ir) = .TRUE.   ! -> send to nei
868                     ENDIF
869                     ! take care of neighbourg in the exterior of the computational domain
870                     IF( iibe == iiout .OR. ijbe == ijout ) THEN   ! Neib outside of the nei MPI domain -> nei cannot compute it
871                        IF( mpiSnei(nn_hls,iSnei) > -1 )   lsend_bdyext(ib_bdy,igrd,iSnei,ir) = .TRUE.   ! -> send to nei
872                     ENDIF
873                  END IF
874
875               END DO   ! 4 sides
876               !
877               ! specific treatment for the corners
878               !
879               DO icnt = 1, 4
880                  SELECT CASE( icnt )
881                     !                                       ...|....
882                  CASE( 1 )   ! x: rim on sw-corner            o|   :
883                     !          o: potential neighbour(s)      o|x__:__
884                     !             outside of the MPI domain   o o o:
885                     !                                              :
886                     iRdiag   = jpsw             ;   iRsono   = jpso             ! Recv: for sw or so
887                     iSdiag   = jpne             ;   iSsono   = jpno             ! Send: to ne or no
888                     iiRcorn  = 1                ;   ijRcorn  = 1                ! receiving sw-corner
889                     iiSdiag  = jpi-2*nn_hls+1   ;   ijSdiag  = jpj-2*nn_hls+1   ! send to sw-corner of ne neighbourg
890                     iiSsono  = 1                ;   ijSsono  = jpj-2*nn_hls+1   ! send to sw-corner of no neighbourg
891                     iioutdir = -1               ;   ijoutdir = -1               ! outside MPI domain: westward or southward
892                     !                                          ....|...
893                  CASE( 2 )   ! x: rim on se-corner             :   |o
894                     !          o: potential neighbour(s)     __:__x|o
895                     !             outside of the MPI domain    :o o o
896                     !                                          :   
897                     iRdiag   = jpse             ;   iRsono   = jpso             ! Recv: for se or so
898                     iSdiag   = jpnw             ;   iSsono   = jpno             ! Send: to nw or no
899                     iiRcorn  = jpi              ;   ijRcorn  = 1                ! receiving se-corner
900                     iiSdiag  =     2*nn_hls     ;   ijSdiag  = jpj-2*nn_hls+1   ! send to se-corner of nw neighbourg
901                     iiSsono  = jpi              ;   ijSsono  = jpj-2*nn_hls+1   ! send to se-corner of no neighbourg
902                     iioutdir = 1                ;   ijoutdir = -1               ! outside MPI domain: eastward or southward
903                     !                                              :       
904                     !                                         o o_o:___
905                  CASE( 3 )   ! x: rim on nw-corner            o|x  :
906                     !          o: potential neighbour(s)    ..o|...:
907                     !             outside of the MPI domain    |
908                     iRdiag   = jpnw             ;   iRsono   = jpno             ! Recv: for nw or no
909                     iSdiag   = jpse             ;   iSsono   = jpso             ! Send: to se or so
910                     iiRcorn  = 1                ;   ijRcorn  = jpj              ! receiving nw-corner
911                     iiSdiag  = jpi-2*nn_hls+1   ;   ijSdiag  =     2*nn_hls     ! send to nw-corner of se neighbourg
912                     iiSsono  = 1                ;   ijSsono  =     2*nn_hls     ! send to nw-corner of so neighbourg
913                     iioutdir = -1               ;   ijoutdir =  1               ! outside MPI domain: westward or northward
914                     !                                          :       
915                     !                                       ___:o_o o
916                  CASE( 4 )   ! x: rim on ne-corner             :  x|o
917                     !          o: potential neighbour(s)       :...|o...
918                     !             outside of the MPI domain        |
919                     iRdiag   = jpne             ;   iRsono   = jpno             ! Recv: for ne or no
920                     iSdiag   = jpsw             ;   iSsono   = jpso             ! Send: to sw or so
921                     iiRcorn  = jpi              ;   ijRcorn  = jpj              ! receiving ne-corner
922                     iiSdiag  =     2*nn_hls     ;   ijSdiag  =     2*nn_hls     ! send to ne-corner of sw neighbourg
923                     iiSsono  = jpi              ;   ijSsono  =     2*nn_hls     ! send to ne-corner of so neighbourg
924                     iioutdir = 1                ;   ijoutdir = 1                ! outside MPI domain: eastward or southward
925                  END SELECT
926                  !
927                  ! Check if we need to receive data for this rim point
928                  IF( ii == iiRcorn .AND. ij == ijRcorn ) THEN   ! the rim point is located on the corner for the MPI domain
929                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the MPI domain?
930                     ! take care of neighbourg(s) in the interior of the computational domain
931                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of the MPI domain
932                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN     ! -> I cannot compute it -> recv it
933                        IF( mpiRnei(nn_hls,iRdiag) > -1                    )   lrecv_bdyint(ib_bdy,igrd,iRdiag,ir) = .TRUE.   ! Receive directly from diagonal neighbourg
934                        IF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 )   lrecv_bdyint(ib_bdy,igrd,iRsono,ir) = .TRUE.   ! Receive through the South/North neighbourg
935                     ENDIF
936                     ! take care of neighbourg in the exterior of the computational domain
937                     IF(  iibe==iiout .OR. ijbe==ijout ) THEN   ! Neib outside of the MPI domain -> I cannot compute it -> recv it
938                        IF( mpiRnei(nn_hls,iRdiag) > -1                    )   lrecv_bdyext(ib_bdy,igrd,iRdiag,ir) = .TRUE.   ! Receive directly from diagonal neighbourg
939                        IF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 )   lrecv_bdyext(ib_bdy,igrd,iRsono,ir) = .TRUE.   ! Receive through the South/North neighbourg
940                     ENDIF
941                  ENDIF
942                  !
943                  ! Check if this rim point corresponds to the corner of one neighbourg. if yes, do we need to send data?
944                  ! Direct send to diag: Is this rim point the corner point of a diag neighbour with which we communicate?
945                  IF( ii == iiSdiag .AND. ij == ijSdiag .AND. mpiSnei(nn_hls,iSdiag) > -1 ) THEN
946                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the nei MPI domain?
947                     ! take care of neighbourg(s) in the interior of the computational domain
948                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of diag nei MPI
949                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout )      &   ! domain -> nei cannot compute it
950                        &    lsend_bdyint(ib_bdy,igrd,iSdiag,ir) = .TRUE.                        ! send rim point data to diag nei
951                     ! take care of neighbourg in the exterior of the computational domain
952                     IF(  iibe==iiout .OR. ijbe==ijout )   &                                 
953                        &    lsend_bdyext(ib_bdy,igrd,iSdiag,ir) = .TRUE.
954                  ENDIF
955                  ! Indirect send to diag (through so/no): rim point is the corner point of a so/no nei with which we communicate
956                  IF( ii == iiSsono .AND. ij == ijSsono .AND. mpiSnei(nn_hls,iSsono) > -1 .AND. nn_comm == 1 ) THEN
957                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the nei MPI domain?
958                     ! take care of neighbourg(s) in the interior of the computational domain
959                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of so/no nei MPI
960                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout )      &   ! domain -> nei cannot compute it
961                        &    lsend_bdyint(ib_bdy,igrd,iSsono,ir) = .TRUE.                        ! send rim point data to so/no nei
962                     ! take care of neighbourg in the exterior of the computational domain
963                     IF(  iibe==iiout .OR. ijbe==ijout )   &
964                        &    lsend_bdyext(ib_bdy,igrd,iSsono,ir) = .TRUE.
965                  ENDIF
966                  !
967               END DO   ! 4 corners
968            END DO   ! ib
969         END DO   ! igrd
970
971         ! Comment out for debug
972!!$         DO ir = 0,1
973!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing,   &
974!!$               &                              lsend = lsend_bdyint(ib_bdy,1,:,ir), lrecv = lrecv_bdyint(ib_bdy,1,:,ir) )
975!!$            IF(lwp) WRITE(numout,*) ' bdy debug int T', ir ; CALL FLUSH(numout)
976!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing,   &
977!!$               &                              lsend = lsend_bdyint(ib_bdy,2,:,ir), lrecv = lrecv_bdyint(ib_bdy,2,:,ir) )
978!!$            IF(lwp) WRITE(numout,*) ' bdy debug int U', ir ; CALL FLUSH(numout)
979!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing,   &
980!!$               &                              lsend = lsend_bdyint(ib_bdy,3,:,ir), lrecv = lrecv_bdyint(ib_bdy,3,:,ir) )   
981!!$            IF(lwp) WRITE(numout,*) ' bdy debug int V', ir ; CALL FLUSH(numout)
982!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing,   &
983!!$               &                              lsend = lsend_bdyext(ib_bdy,1,:,ir), lrecv = lrecv_bdyext(ib_bdy,1,:,ir) )
984!!$            IF(lwp) WRITE(numout,*) ' bdy debug ext T', ir ; CALL FLUSH(numout)
985!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing,   &
986!!$               &                              lsend = lsend_bdyext(ib_bdy,2,:,ir), lrecv = lrecv_bdyext(ib_bdy,2,:,ir) )
987!!$            IF(lwp) WRITE(numout,*) ' bdy debug ext U', ir ; CALL FLUSH(numout)
988!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing,   &
989!!$               &                              lsend = lsend_bdyext(ib_bdy,3,:,ir), lrecv = lrecv_bdyext(ib_bdy,3,:,ir) )   
990!!$            IF(lwp) WRITE(numout,*) ' bdy debug ext V', ir ; CALL FLUSH(numout)
991!!$         END DO
992         
993      END DO   ! ib_bdy
994
995      DO ib_bdy = 1,nb_bdy
996         IF(  cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &
997            & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &
998            & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN
999            DO igrd = 1, jpbgrd
1000               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1001                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1002                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1003                  IF(  mig0(ii) > 2 .AND. mig0(ii) < Ni0glo-2 .AND. mjg0(ij) > 2 .AND. mjg0(ij) < Nj0glo-2  ) THEN
1004                     WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain'
1005                     CALL ctl_stop( ctmp1 )
1006                  END IF
1007               END DO
1008            END DO
1009         END IF
1010      END DO
1011      !
1012      DEALLOCATE( nbidta, nbjdta, nbrdta )
1013      !
1014   END SUBROUTINE bdy_def
1015
1016
1017   SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 )
1018      !!----------------------------------------------------------------------
1019      !!                 ***  ROUTINE bdy_rim_treat  ***
1020      !!
1021      !! ** Purpose :   Initialize structures ( flagu, flagv, ntreat ) indicating how rim points
1022      !!                  are to be handled in the boundary condition treatment
1023      !!
1024      !! ** Method  :   - to handle rim 0 zmasks must indicate ocean points      (set at one on rim 0 and rim 1 and interior)
1025      !!                            and bdymasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
1026      !!                    (as if rim 1 was free ocean)
1027      !!                - to handle rim 1 zmasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
1028      !!                            and bdymasks must indicate free ocean points (set at one on interior)
1029      !!                    (as if rim 0 was land)
1030      !!                - we can then check in which direction the interior of the computational domain is with the difference
1031      !!                         mask array values on both sides to compute flagu and flagv
1032      !!                - and look at the ocean neighbours to compute ntreat
1033      !!----------------------------------------------------------------------
1034      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pumask, pvmask   ! temporary u/v mask array
1035      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pfmask           ! temporary fmask excluding coastal boundary condition (shlat)
1036      LOGICAL                             , INTENT (in   ) :: lrim0            ! .true. -> rim 0   .false. -> rim 1
1037      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, icount       ! dummy loop indices
1038      INTEGER  ::   i_offset, j_offset, inn                ! local integer
1039      INTEGER  ::   ibeg, iend                             ! local integer
1040      LOGICAL  ::   llnon, llson, llean, llwen             ! local logicals indicating the presence of a ocean neighbour
1041      REAL(wp), POINTER, DIMENSION(:,:)       ::   zmask   ! pointer to 2D mask fields
1042      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
1043      CHARACTER(LEN=1), DIMENSION(jpbgrd)     ::   cgrid
1044      REAL(wp)        , DIMENSION(jpi,jpj)    ::   ztmp
1045      !!----------------------------------------------------------------------
1046
1047      cgrid = (/'t','u','v'/)
1048
1049      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
1050
1051         DO igrd = 1, jpbgrd
1052           
1053            IF( lrim0 ) THEN   ! extent of rim 0
1054               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
1055            ELSE               ! extent of rim 1
1056               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
1057            END IF
1058
1059            ! Calculate relationship of U direction to the local orientation of the boundary
1060            ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
1061            ! flagu =  0 : u is tangential
1062            ! flagu =  1 : u is normal to the boundary and is direction is inward
1063            SELECT CASE( igrd )
1064               CASE( 1 )   ;   zmask => pumask     ;   i_offset = 0   ! U(i-1)   T(i)   U(i  )
1065               CASE( 2 )   ;   zmask => bdytmask   ;   i_offset = 1   ! T(i  )   U(i)   T(i+1)
1066               CASE( 3 )   ;   zmask => pfmask     ;   i_offset = 0   ! F(i-1)   V(i)   F(i  )
1067            END SELECT
1068            icount = 0
1069            ztmp(:,:) = -999._wp
1070            DO ib = ibeg, iend 
1071               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1072               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1073               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
1074               zwfl = zmask(ii+i_offset-1,ij)
1075               zefl = zmask(ii+i_offset  ,ij)
1076               ! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
1077               IF( i_offset == 1 .and. zefl + zwfl == 2._wp ) THEN
1078                  icount = icount + 1
1079                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
1080               ELSE
1081                  ztmp(ii,ij) = -zwfl + zefl
1082               ENDIF
1083            END DO
1084            IF( icount /= 0 ) THEN
1085               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
1086                  ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1087               CALL ctl_stop( ctmp1 )
1088            ENDIF
1089            SELECT CASE( igrd )
1090               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1091               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1092               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1093            END SELECT
1094            DO ib = ibeg, iend
1095               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1096               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1097               idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij)
1098            END DO
1099           
1100            ! Calculate relationship of V direction to the local orientation of the boundary
1101            ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
1102            ! flagv =  0 : v is tangential
1103            ! flagv =  1 : v is normal to the boundary and is direction is inward
1104            SELECT CASE( igrd )
1105               CASE( 1 )   ;   zmask => pvmask     ;   j_offset = 0
1106               CASE( 2 )   ;   zmask => pfmask     ;   j_offset = 0
1107               CASE( 3 )   ;   zmask => bdytmask   ;   j_offset = 1
1108            END SELECT
1109            icount = 0
1110            ztmp(:,:) = -999._wp
1111            DO ib = ibeg, iend
1112               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1113               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1114               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
1115               zsfl = zmask(ii,ij+j_offset-1)
1116               znfl = zmask(ii,ij+j_offset  )
1117               ! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
1118               IF( j_offset == 1 .and. znfl + zsfl == 2._wp ) THEN
1119                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
1120                  icount = icount + 1
1121               ELSE
1122                  ztmp(ii,ij) = -zsfl + znfl
1123               END IF
1124            END DO
1125            IF( icount /= 0 ) THEN
1126               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
1127                  ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1128               CALL ctl_stop( ctmp1 )
1129            ENDIF
1130            SELECT CASE( igrd )
1131               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1132               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1133               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1134            END SELECT
1135            DO ib = ibeg, iend
1136               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1137               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1138               idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij)
1139            END DO
1140     
1141            ! Calculate ntreat
1142            SELECT CASE( igrd )
1143               CASE( 1 )   ;   zmask => bdytmask 
1144               CASE( 2 )   ;   zmask => bdyumask 
1145               CASE( 3 )   ;   zmask => bdyvmask 
1146            END SELECT
1147            ztmp(:,:) = -999._wp
1148            DO ib = ibeg, iend
1149               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1150               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1151               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
1152               llnon = zmask(ii  ,ij+1) == 1._wp 
1153               llson = zmask(ii  ,ij-1) == 1._wp 
1154               llean = zmask(ii+1,ij  ) == 1._wp 
1155               llwen = zmask(ii-1,ij  ) == 1._wp 
1156               inn  = COUNT( (/ llnon, llson, llean, llwen /) )
1157               IF( inn == 0 ) THEN   ! no neighbours -> interior of a corner  or  cluster of rim points
1158                  !               !              !     _____     !     _____    !    __     __
1159                  !  1 |   o      !  2  o   |    !  3 | x        !  4     x |   !      |   |   -> error
1160                  !    |_x_ _     !    _ _x_|    !    |   o      !      o   |   !      |x_x|
1161                  IF(     zmask(ii+1,ij+1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 1._wp
1162                  ELSEIF( zmask(ii-1,ij+1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 2._wp
1163                  ELSEIF( zmask(ii+1,ij-1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 3._wp
1164                  ELSEIF( zmask(ii-1,ij-1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 4._wp
1165                  ELSE                                       ;   ztmp(ii,ij) = -1._wp
1166                     WRITE(ctmp1,*) 'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
1167                       ' on boundary set ', ib_bdy, ' has no free ocean neighbour'
1168                     IF( lrim0 ) THEN
1169                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 0 points.'
1170                     ELSE
1171                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.'
1172                     END IF
1173                     CALL ctl_warn( ctmp1, ctmp2 )
1174                  END IF
1175               END IF
1176               IF( inn == 1 ) THEN   ! middle of linear bdy  or incomplete corner  ! ___ o
1177                  !    |         !         |   !      o     !    ______            !    |x___
1178                  ! 5  | x o     ! 6   o x |   ! 7  __x__   ! 8    x
1179                  !    |         !         |   !            !      o
1180                  IF( llean )   ztmp(ii,ij) = 5._wp
1181                  IF( llwen )   ztmp(ii,ij) = 6._wp
1182                  IF( llnon )   ztmp(ii,ij) = 7._wp
1183                  IF( llson )   ztmp(ii,ij) = 8._wp
1184               END IF
1185               IF( inn == 2 ) THEN   ! exterior of a corner
1186                  !        o      !        o      !    _____|       !       |_____ 
1187                  !  9 ____x o    ! 10   o x___   ! 11     x o      ! 12   o x     
1188                  !         |     !       |       !        o        !        o
1189                  IF( llnon .AND. llean )   ztmp(ii,ij) =  9._wp
1190                  IF( llnon .AND. llwen )   ztmp(ii,ij) = 10._wp
1191                  IF( llson .AND. llean )   ztmp(ii,ij) = 11._wp
1192                  IF( llson .AND. llwen )   ztmp(ii,ij) = 12._wp
1193               END IF
1194               IF( inn == 3 ) THEN   ! 3 neighbours     __   __
1195                  !    |_  o      !        o  _|  !       |_|     !       o         
1196                  ! 13  _| x o    ! 14   o x |_   ! 15   o x o    ! 16  o x o       
1197                  !    |   o      !        o   |  !        o      !    __|¨|__   
1198                  IF( llnon .AND. llean .AND. llson )   ztmp(ii,ij) = 13._wp
1199                  IF( llnon .AND. llwen .AND. llson )   ztmp(ii,ij) = 14._wp
1200                  IF( llwen .AND. llson .AND. llean )   ztmp(ii,ij) = 15._wp
1201                  IF( llwen .AND. llnon .AND. llean )   ztmp(ii,ij) = 16._wp
1202               END IF
1203               IF( inn == 4 ) THEN
1204                  WRITE(ctmp1,*)  'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
1205                       ' on boundary set ', ib_bdy, ' have 4 neighbours'
1206                  CALL ctl_stop( ctmp1 )
1207               END IF
1208            END DO
1209            SELECT CASE( igrd )
1210               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1211               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1212               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1213            END SELECT
1214            DO ib = ibeg, iend
1215               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1216               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1217               idx_bdy(ib_bdy)%ntreat(ib,igrd) = NINT(ztmp(ii,ij))
1218            END DO
1219            !
1220         END DO   ! jpbgrd
1221         !
1222      END DO   ! ib_bdy
1223
1224    END SUBROUTINE bdy_rim_treat
1225
1226   
1227    SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )
1228      !!----------------------------------------------------------------------
1229      !!                 ***  ROUTINE find_neib  ***
1230      !!
1231      !! ** Purpose :   get ii1, ij1, ii2, ij2, ii3, ij3, the indices of
1232      !!               the free ocean neighbours of (ii,ij) for bdy treatment
1233      !!
1234      !! ** Method  :  use itreat input to select a case
1235      !!               N.B. ntreat is defined for all bdy points in routine bdy_rim_treat
1236      !!
1237      !!----------------------------------------------------------------------
1238      INTEGER, INTENT(in   )      ::   ii, ij, itreat
1239      INTEGER, INTENT(  out)      ::   ii1, ij1, ii2, ij2, ii3, ij3
1240      !!----------------------------------------------------------------------
1241      SELECT CASE( itreat )   ! points that will be used by bdy routines, -1 will be discarded
1242         !               !               !     _____     !     _____     
1243         !  1 |   o      !  2  o   |     !  3 | x        !  4     x |   
1244         !    |_x_ _     !    _ _x_|     !    |   o      !      o   |
1245      CASE( 1 )    ;   ii1 = ii+1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1246      CASE( 2 )    ;   ii1 = ii-1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1247      CASE( 3 )    ;   ii1 = ii+1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1248      CASE( 4 )    ;   ii1 = ii-1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1249         !    |          !         |     !      o        !    ______                   ! or incomplete corner
1250         ! 5  | x o      ! 6   o x |     ! 7  __x__      ! 8    x                      !  7  ____ o
1251         !    |          !         |     !               !      o                      !         |x___
1252      CASE( 5 )    ;   ii1 = ii+1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1253      CASE( 6 )    ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1254      CASE( 7 )    ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1255      CASE( 8 )    ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1256         !        o      !        o      !    _____|     !       |_____ 
1257         !  9 ____x o    ! 10   o x___   ! 11     x o    ! 12   o x     
1258         !         |     !       |       !        o      !        o     
1259      CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
1260      CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1261      CASE( 11 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1262      CASE( 12 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1263         !    |_  o      !        o  _|  !     ¨¨|_|¨¨   !       o         
1264         ! 13  _| x o    !  14  o x |_   !  15  o x o    ! 16  o x o       
1265         !    |   o      !        o   |  !        o      !    __|¨|__
1266      CASE( 13 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1   
1267      CASE( 14 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1 
1268      CASE( 15 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij   
1269      CASE( 16 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij
1270      END SELECT
1271   END SUBROUTINE find_neib
1272   
1273
1274   SUBROUTINE bdy_read_seg( kb_bdy, knblendta ) 
1275      !!----------------------------------------------------------------------
1276      !!                 ***  ROUTINE bdy_read_seg  ***
1277      !!
1278      !! ** Purpose :  build bdy coordinates with segments defined in namelist
1279      !!
1280      !! ** Method  :  read namelist nambdy_index blocks
1281      !!
1282      !!----------------------------------------------------------------------
1283      INTEGER                   , INTENT (in   ) ::   kb_bdy           ! bdy number
1284      INTEGER, DIMENSION(jpbgrd), INTENT (  out) ::   knblendta        ! length of index arrays
1285      !!
1286      INTEGER          ::   ios                 ! Local integer output status for namelist read
1287      INTEGER          ::   nbdyind, nbdybeg, nbdyend
1288      INTEGER          ::   nbdy_count, nbdy_rdstart, nbdy_loc
1289      CHARACTER(LEN=1) ::   ctypebdy   !     -        -
1290      CHARACTER(LEN=50)::   cerrmsg    !     -        -
1291      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
1292      !!----------------------------------------------------------------------
1293      ! Need to support possibility of reading more than one nambdy_index from
1294      ! the namelist_cfg internal file.
1295      ! Do this by finding the kb_bdy'th occurence of nambdy_index in the
1296      ! character buffer as the starting point.
1297      nbdy_rdstart = 1
1298      DO nbdy_count = 1, kb_bdy
1299       nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_index' )
1300       IF( nbdy_loc .GT. 0 ) THEN
1301          nbdy_rdstart = nbdy_rdstart + nbdy_loc
1302       ELSE
1303          WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',kb_bdy,' of nambdy_index not found'
1304          ios = -1
1305          CALL ctl_nam ( ios , cerrmsg )
1306       ENDIF
1307      END DO
1308      nbdy_rdstart = MAX( 1, nbdy_rdstart - 2 )
1309      READ  ( numnam_cfg( nbdy_rdstart: ), nambdy_index, IOSTAT = ios, ERR = 904)
1310904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' )
1311      IF(lwm) WRITE ( numond, nambdy_index )
1312     
1313      SELECT CASE ( TRIM(ctypebdy) )
1314      CASE( 'N' )
1315         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1316            nbdyind  = Nj0glo - 2  ! set boundary to whole side of model domain.
1317            nbdybeg  = 2
1318            nbdyend  = Ni0glo - 1
1319         ENDIF
1320         nbdysegn = nbdysegn + 1
1321         npckgn(nbdysegn) = kb_bdy ! Save bdy package number
1322         jpjnob(nbdysegn) = nbdyind 
1323         jpindt(nbdysegn) = nbdybeg
1324         jpinft(nbdysegn) = nbdyend
1325         !
1326      CASE( 'S' )
1327         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1328            nbdyind  = 2           ! set boundary to whole side of model domain.
1329            nbdybeg  = 2
1330            nbdyend  = Ni0glo - 1
1331         ENDIF
1332         nbdysegs = nbdysegs + 1
1333         npckgs(nbdysegs) = kb_bdy ! Save bdy package number
1334         jpjsob(nbdysegs) = nbdyind
1335         jpisdt(nbdysegs) = nbdybeg
1336         jpisft(nbdysegs) = nbdyend
1337         !
1338      CASE( 'E' )
1339         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1340            nbdyind  = Ni0glo - 2  ! set boundary to whole side of model domain.
1341            nbdybeg  = 2
1342            nbdyend  = Nj0glo - 1
1343         ENDIF
1344         nbdysege = nbdysege + 1 
1345         npckge(nbdysege) = kb_bdy ! Save bdy package number
1346         jpieob(nbdysege) = nbdyind
1347         jpjedt(nbdysege) = nbdybeg
1348         jpjeft(nbdysege) = nbdyend
1349         !
1350      CASE( 'W' )
1351         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1352            nbdyind  = 2           ! set boundary to whole side of model domain.
1353            nbdybeg  = 2
1354            nbdyend  = Nj0glo - 1
1355         ENDIF
1356         nbdysegw = nbdysegw + 1
1357         npckgw(nbdysegw) = kb_bdy ! Save bdy package number
1358         jpiwob(nbdysegw) = nbdyind
1359         jpjwdt(nbdysegw) = nbdybeg
1360         jpjwft(nbdysegw) = nbdyend
1361         !
1362      CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
1363      END SELECT
1364     
1365      ! For simplicity we assume that in case of straight bdy, arrays have the same length
1366      ! (even if it is true that last tangential velocity points
1367      ! are useless). This simplifies a little bit boundary data format (and agrees with format
1368      ! used so far in obc package)
1369     
1370      knblendta(1:jpbgrd) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy)
1371     
1372   END SUBROUTINE bdy_read_seg
1373
1374   
1375   SUBROUTINE bdy_ctl_seg
1376      !!----------------------------------------------------------------------
1377      !!                 ***  ROUTINE bdy_ctl_seg  ***
1378      !!
1379      !! ** Purpose :   Check straight open boundary segments location
1380      !!
1381      !! ** Method  :   - Look for open boundary corners
1382      !!                - Check that segments start or end on land
1383      !!----------------------------------------------------------------------
1384      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest 
1385      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns 
1386      REAL(wp), DIMENSION(2) ::   ztestmask
1387      !!----------------------------------------------------------------------
1388      !
1389      IF (lwp) WRITE(numout,*) ' '
1390      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1391      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1392      !
1393      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege
1394      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw
1395      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn
1396      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs
1397      !
1398      ! 1. Check bounds
1399      !----------------
1400      DO ib = 1, nbdysegn
1401         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
1402         IF ((jpjnob(ib).ge.Nj0glo-1).or.& 
1403            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1404         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1405         IF (jpindt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1406         IF (jpinft(ib).gt.Ni0glo)     CALL ctl_stop( 'End index out of domain' )
1407      END DO
1408      !
1409      DO ib = 1, nbdysegs
1410         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
1411         IF ((jpjsob(ib).ge.Nj0glo-1).or.& 
1412            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1413         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1414         IF (jpisdt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1415         IF (jpisft(ib).gt.Ni0glo)     CALL ctl_stop( 'End index out of domain' )
1416      END DO
1417      !
1418      DO ib = 1, nbdysege
1419         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib)
1420         IF ((jpieob(ib).ge.Ni0glo-1).or.& 
1421            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1422         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1423         IF (jpjedt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1424         IF (jpjeft(ib).gt.Nj0glo)     CALL ctl_stop( 'End index out of domain' )
1425      END DO
1426      !
1427      DO ib = 1, nbdysegw
1428         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib)
1429         IF ((jpiwob(ib).ge.Ni0glo-1).or.& 
1430            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1431         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1432         IF (jpjwdt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1433         IF (jpjwft(ib).gt.Nj0glo)     CALL ctl_stop( 'End index out of domain' )
1434      ENDDO
1435      !     
1436      ! 2. Look for segment crossings
1437      !------------------------------
1438      IF (lwp) WRITE(numout,*) '**Look for segments corners  :'
1439      !
1440      itest = 0 ! corner number
1441      !
1442      ! flag to detect if start or end of open boundary belongs to a corner
1443      ! if not (=0), it must be on land.
1444      ! if a corner is detected, save bdy package number for further tests
1445      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1446      ! South/West crossings
1447      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1448         DO ib1 = 1, nbdysegw       
1449            DO ib2 = 1, nbdysegs
1450               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1451                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1452                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1453                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1454                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN 
1455                     ! We have a possible South-West corner                     
1456!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1457!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1458                     icornw(ib1,1) = npckgs(ib2)
1459                     icorns(ib2,1) = npckgw(ib1)
1460                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1461                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1462                        &                                     jpisft(ib2), jpjwft(ib1)
1463                     WRITE(ctmp2,*) ' Not allowed yet'
1464                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 
1465                        &                            ' and South segment: ',npckgs(ib2)
1466                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1467                  ELSE
1468                     WRITE(ctmp1,*) ' Check South and West Open boundary indices'
1469                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , &
1470                        &                            ' and South segment: ',npckgs(ib2)
1471                     CALL ctl_stop( ctmp1, ctmp2 )
1472                  END IF
1473               END IF
1474            END DO
1475         END DO
1476      END IF
1477      !
1478      ! South/East crossings
1479      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1480         DO ib1 = 1, nbdysege
1481            DO ib2 = 1, nbdysegs
1482               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1483                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1484                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. &
1485                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN
1486                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1487                     ! We have a possible South-East corner
1488!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1489!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1490                     icorne(ib1,1) = npckgs(ib2)
1491                     icorns(ib2,2) = npckge(ib1)
1492                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1493                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1494                        &                                     jpisdt(ib2), jpjeft(ib1)
1495                     WRITE(ctmp2,*) ' Not allowed yet'
1496                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1497                        &                            ' and South segment: ',npckgs(ib2)
1498                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1499                  ELSE
1500                     WRITE(ctmp1,*) ' Check South and East Open boundary indices'
1501                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1502                     &                               ' and South segment: ',npckgs(ib2)
1503                     CALL ctl_stop( ctmp1, ctmp2 )
1504                  END IF
1505               END IF
1506            END DO
1507         END DO
1508      END IF
1509      !
1510      ! North/West crossings
1511      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1512         DO ib1 = 1, nbdysegw       
1513            DO ib2 = 1, nbdysegn
1514               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. &
1515                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. &
1516                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1517                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1518                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1519                     ! We have a possible North-West corner
1520!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1521!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1522                     icornw(ib1,2) = npckgn(ib2)
1523                     icornn(ib2,1) = npckgw(ib1)
1524                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1525                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1526                     &                                     jpinft(ib2), jpjwdt(ib1)
1527                     WRITE(ctmp2,*) ' Not allowed yet'
1528                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1529                     &                               ' and North segment: ',npckgn(ib2)
1530                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1531                  ELSE
1532                     WRITE(ctmp1,*) ' Check North and West Open boundary indices'
1533                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1534                     &                               ' and North segment: ',npckgn(ib2)
1535                     CALL ctl_stop( ctmp1, ctmp2 )
1536                  END IF
1537               END IF
1538            END DO
1539         END DO
1540      END IF
1541      !
1542      ! North/East crossings
1543      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1544         DO ib1 = 1, nbdysege       
1545            DO ib2 = 1, nbdysegn
1546               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1547                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1548                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1549                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1550                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1551                     ! We have a possible North-East corner
1552!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1553!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1554                     icorne(ib1,2) = npckgn(ib2)
1555                     icornn(ib2,2) = npckge(ib1)
1556                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1557                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1558                     &                                     jpindt(ib2), jpjedt(ib1)
1559                     WRITE(ctmp2,*) ' Not allowed yet'
1560                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1561                     &                               ' and North segment: ',npckgn(ib2)
1562                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1563                  ELSE
1564                     WRITE(ctmp1,*) ' Check North and East Open boundary indices'
1565                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1566                     &                               ' and North segment: ',npckgn(ib2)
1567                     CALL ctl_stop( ctmp1, ctmp2 )
1568                  END IF
1569               END IF
1570            END DO
1571         END DO
1572      END IF
1573      !
1574      ! 3. Check if segment extremities are on land
1575      !--------------------------------------------
1576      !
1577      ! West segments
1578      DO ib = 1, nbdysegw
1579         ! get mask at boundary extremities:
1580         ztestmask(1:2)=0.
1581         DO ji = 1, jpi
1582            DO jj = 1, jpj             
1583              IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1584              IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1585            END DO
1586         END DO
1587         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1588
1589         IF (ztestmask(1)==1) THEN
1590            IF (icornw(ib,1)==0) THEN
1591               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1592               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1593            ELSE
1594               ! This is a corner
1595               IF(lwp) WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1596               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1597               itest=itest+1
1598            ENDIF
1599         ENDIF
1600         IF (ztestmask(2)==1) THEN
1601            IF (icornw(ib,2)==0) THEN
1602               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1603               CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' )
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( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjedt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1620              IF( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjeft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1621            END DO
1622         END DO
1623         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1624
1625         IF (ztestmask(1)==1) THEN
1626            IF (icorne(ib,1)==0) THEN
1627               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1628               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1629            ELSE
1630               ! This is a corner
1631               IF(lwp) WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1632               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1633               itest=itest+1
1634            ENDIF
1635         ENDIF
1636         IF (ztestmask(2)==1) THEN
1637            IF (icorne(ib,2)==0) THEN
1638               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1639               CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1640            ELSE
1641               ! This is a corner
1642               IF(lwp) WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1643               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1644               itest=itest+1
1645            ENDIF
1646         ENDIF
1647      END DO
1648      !
1649      ! South segments
1650      DO ib = 1, nbdysegs
1651         ! get mask at boundary extremities:
1652         ztestmask(1:2)=0.
1653         DO ji = 1, jpi
1654            DO jj = 1, jpj             
1655              IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1656              IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1657            END DO
1658         END DO
1659         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1660
1661         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1662            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1663            CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1664         ENDIF
1665         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1666            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1667            CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1668         ENDIF
1669      END DO
1670      !
1671      ! North segments
1672      DO ib = 1, nbdysegn
1673         ! get mask at boundary extremities:
1674         ztestmask(1:2)=0.
1675         DO ji = 1, jpi
1676            DO jj = 1, jpj             
1677               IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpindt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1678               IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpinft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1679            END DO
1680         END DO
1681         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1682
1683         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1684            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1685            CALL ctl_stop( ctmp1, ' does not start on land' )
1686         ENDIF
1687         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1688            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1689            CALL ctl_stop( ctmp1, ' does not end on land' )
1690         ENDIF
1691      END DO
1692      !
1693      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1694      !
1695      ! Other tests TBD:
1696      ! segments completly on land
1697      ! optimized open boundary array length according to landmask
1698      ! Nudging layers that overlap with interior domain
1699      !
1700   END SUBROUTINE bdy_ctl_seg
1701
1702   
1703   SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta ) 
1704      !!----------------------------------------------------------------------
1705      !!                 ***  ROUTINE bdy_coords_seg  ***
1706      !!
1707      !! ** Purpose :  build nbidta, nbidta, nbrdta for bdy built with segments
1708      !!
1709      !! ** Method  : 
1710      !!
1711      !!----------------------------------------------------------------------
1712      INTEGER, DIMENSION(:,:,:), intent(  out)  ::   nbidta, nbjdta, nbrdta   ! Index arrays: i and j indices of bdy dta
1713      !!
1714      INTEGER  ::   ii, ij, ir, iseg
1715      INTEGER  ::   igrd         ! grid type (t=1, u=2, v=3)
1716      INTEGER  ::   icount       !
1717      INTEGER  ::   ib_bdy       ! bdy number
1718      !!----------------------------------------------------------------------
1719
1720      ! East
1721      !-----
1722      DO iseg = 1, nbdysege
1723         ib_bdy = npckge(iseg)
1724         !
1725         ! ------------ T points -------------
1726         igrd=1
1727         icount=0
1728         DO ir = 1, nn_rimwidth(ib_bdy)
1729            DO ij = jpjedt(iseg), jpjeft(iseg)
1730               icount = icount + 1
1731               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir + nn_hls
1732               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1733               nbrdta(icount, igrd, ib_bdy) = ir
1734            ENDDO
1735         ENDDO
1736         !
1737         ! ------------ U points -------------
1738         igrd=2
1739         icount=0
1740         DO ir = 1, nn_rimwidth(ib_bdy)
1741            DO ij = jpjedt(iseg), jpjeft(iseg)
1742               icount = icount + 1
1743               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir + nn_hls
1744               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1745               nbrdta(icount, igrd, ib_bdy) = ir
1746            ENDDO
1747         ENDDO
1748         !
1749         ! ------------ V points -------------
1750         igrd=3
1751         icount=0
1752         DO ir = 1, nn_rimwidth(ib_bdy)
1753            !            DO ij = jpjedt(iseg), jpjeft(iseg) - 1
1754            DO ij = jpjedt(iseg), jpjeft(iseg)
1755               icount = icount + 1
1756               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir + nn_hls
1757               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1758               nbrdta(icount, igrd, ib_bdy) = ir
1759            ENDDO
1760            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1761            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1762         ENDDO
1763      ENDDO
1764      !
1765      ! West
1766      !-----
1767      DO iseg = 1, nbdysegw
1768         ib_bdy = npckgw(iseg)
1769         !
1770         ! ------------ T points -------------
1771         igrd=1
1772         icount=0
1773         DO ir = 1, nn_rimwidth(ib_bdy)
1774            DO ij = jpjwdt(iseg), jpjwft(iseg)
1775               icount = icount + 1
1776               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1777               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1778               nbrdta(icount, igrd, ib_bdy) = ir
1779            ENDDO
1780         ENDDO
1781         !
1782         ! ------------ U points -------------
1783         igrd=2
1784         icount=0
1785         DO ir = 1, nn_rimwidth(ib_bdy)
1786            DO ij = jpjwdt(iseg), jpjwft(iseg)
1787               icount = icount + 1
1788               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1789               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1790               nbrdta(icount, igrd, ib_bdy) = ir
1791            ENDDO
1792         ENDDO
1793         !
1794         ! ------------ V points -------------
1795         igrd=3
1796         icount=0
1797         DO ir = 1, nn_rimwidth(ib_bdy)
1798            !            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
1799            DO ij = jpjwdt(iseg), jpjwft(iseg)
1800               icount = icount + 1
1801               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1802               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1803               nbrdta(icount, igrd, ib_bdy) = ir
1804            ENDDO
1805            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1806            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1807         ENDDO
1808      ENDDO
1809      !
1810      ! North
1811      !-----
1812      DO iseg = 1, nbdysegn
1813         ib_bdy = npckgn(iseg)
1814         !
1815         ! ------------ T points -------------
1816         igrd=1
1817         icount=0
1818         DO ir = 1, nn_rimwidth(ib_bdy)
1819            DO ii = jpindt(iseg), jpinft(iseg)
1820               icount = icount + 1
1821               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1822               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir + nn_hls 
1823               nbrdta(icount, igrd, ib_bdy) = ir
1824            ENDDO
1825         ENDDO
1826         !
1827         ! ------------ U points -------------
1828         igrd=2
1829         icount=0
1830         DO ir = 1, nn_rimwidth(ib_bdy)
1831            !            DO ii = jpindt(iseg), jpinft(iseg) - 1
1832            DO ii = jpindt(iseg), jpinft(iseg)
1833               icount = icount + 1
1834               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1835               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir + nn_hls
1836               nbrdta(icount, igrd, ib_bdy) = ir
1837            ENDDO
1838            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1839            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1840         ENDDO
1841         !
1842         ! ------------ V points -------------
1843         igrd=3
1844         icount=0
1845         DO ir = 1, nn_rimwidth(ib_bdy)
1846            DO ii = jpindt(iseg), jpinft(iseg)
1847               icount = icount + 1
1848               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1849               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir + nn_hls
1850               nbrdta(icount, igrd, ib_bdy) = ir
1851            ENDDO
1852         ENDDO
1853      ENDDO
1854      !
1855      ! South
1856      !-----
1857      DO iseg = 1, nbdysegs
1858         ib_bdy = npckgs(iseg)
1859         !
1860         ! ------------ T points -------------
1861         igrd=1
1862         icount=0
1863         DO ir = 1, nn_rimwidth(ib_bdy)
1864            DO ii = jpisdt(iseg), jpisft(iseg)
1865               icount = icount + 1
1866               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1867               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1868               nbrdta(icount, igrd, ib_bdy) = ir
1869            ENDDO
1870         ENDDO
1871         !
1872         ! ------------ U points -------------
1873         igrd=2
1874         icount=0
1875         DO ir = 1, nn_rimwidth(ib_bdy)
1876            !            DO ii = jpisdt(iseg), jpisft(iseg) - 1
1877            DO ii = jpisdt(iseg), jpisft(iseg)
1878               icount = icount + 1
1879               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1880               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1881               nbrdta(icount, igrd, ib_bdy) = ir
1882            ENDDO
1883            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1884            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1885         ENDDO
1886         !
1887         ! ------------ V points -------------
1888         igrd=3
1889         icount=0
1890         DO ir = 1, nn_rimwidth(ib_bdy)
1891            DO ii = jpisdt(iseg), jpisft(iseg)
1892               icount = icount + 1
1893               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1894               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1895               nbrdta(icount, igrd, ib_bdy) = ir
1896            ENDDO
1897         ENDDO
1898      ENDDO
1899
1900     
1901   END SUBROUTINE bdy_coords_seg
1902   
1903   
1904   SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1905      !!----------------------------------------------------------------------
1906      !!                 ***  ROUTINE bdy_ctl_corn  ***
1907      !!
1908      !! ** Purpose :   Check numerical schemes consistency between
1909      !!                segments having a common corner
1910      !!
1911      !! ** Method  :   
1912      !!----------------------------------------------------------------------
1913      INTEGER, INTENT(in)  ::   ib1, ib2
1914      INTEGER :: itest
1915      !!----------------------------------------------------------------------
1916      itest = 0
1917
1918      IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) )   itest = itest + 1
1919      IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) )   itest = itest + 1
1920      IF( cn_tra  (ib1) /= cn_tra  (ib2) )   itest = itest + 1
1921      !
1922      IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) )   itest = itest + 1
1923      IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) )   itest = itest + 1
1924      IF( nn_tra_dta  (ib1) /= nn_tra_dta  (ib2) )   itest = itest + 1
1925      !
1926      IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) )   itest = itest + 1   
1927      !
1928      IF( itest>0 ) THEN
1929         WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2
1930         CALL ctl_stop( ctmp1, ' have different open bdy schemes' )
1931      ENDIF
1932      !
1933   END SUBROUTINE bdy_ctl_corn
1934
1935
1936   SUBROUTINE bdy_meshwri()
1937      !!----------------------------------------------------------------------
1938      !!                 ***  ROUTINE bdy_meshwri  ***
1939      !!         
1940      !! ** Purpose :   write netcdf file with nbr, flagu, flagv, ntreat for T, U
1941      !!                and V points in 2D arrays for easier visualisation/control
1942      !!
1943      !! ** Method  :   use iom_rstput as in domwri.F
1944      !!----------------------------------------------------------------------     
1945      INTEGER  ::   ib_bdy, ii, ij, igrd, ib     ! dummy loop indices
1946      INTEGER  ::   inum                                   !   -       -
1947      REAL(wp), POINTER, DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields
1948      REAL(wp)         , DIMENSION(jpi,jpj) ::   ztmp
1949      CHARACTER(LEN=1) , DIMENSION(jpbgrd)  ::   cgrid
1950      !!----------------------------------------------------------------------     
1951      cgrid = (/'t','u','v'/)
1952      CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. )
1953      DO igrd = 1, jpbgrd
1954         SELECT CASE( igrd )
1955         CASE( 1 )   ;   zmask => tmask(:,:,1)
1956         CASE( 2 )   ;   zmask => umask(:,:,1)
1957         CASE( 3 )   ;   zmask => vmask(:,:,1)
1958         END SELECT
1959         ztmp(:,:) = zmask(:,:)
1960         DO ib_bdy = 1, nb_bdy
1961            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)      ! nbr deined for all rims
1962               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1963               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1964               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10.
1965               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1966            END DO
1967         END DO
1968         CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1969         ztmp(:,:) = zmask(:,:)
1970         DO ib_bdy = 1, nb_bdy
1971            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagu defined only for rims 0 and 1
1972               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1973               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1974               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10.
1975               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1976            END DO
1977         END DO
1978         CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1979         ztmp(:,:) = zmask(:,:)
1980         DO ib_bdy = 1, nb_bdy
1981            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagv defined only for rims 0 and 1
1982               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1983               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1984               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10.
1985               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1986            END DO
1987         END DO
1988         CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1989         ztmp(:,:) = zmask(:,:)
1990         DO ib_bdy = 1, nb_bdy
1991            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! ntreat defined only for rims 0 and 1
1992               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1993               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1994               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10.
1995               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1996            END DO
1997         END DO
1998         CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1999      END DO
2000      CALL iom_close( inum )
2001
2002   END SUBROUTINE bdy_meshwri
2003   
2004   !!=================================================================================
2005END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.