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 @ 15368

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

trunk: final version (hopefully) for ticket #2731

  • Property svn:keywords set to Id
File size: 109.4 KB
Line 
1MODULE bdyini
2   !!======================================================================
3   !!                       ***  MODULE  bdyini  ***
4   !! Unstructured open boundaries : initialisation
5   !!======================================================================
6   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
7   !!             -   !  2007-01  (D. Storkey) Update to use IOM module
8   !!             -   !  2007-01  (D. Storkey) Tidal forcing
9   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
10   !!            3.3  !  2010-09  (E.O'Dea) updates for Shelf configurations
11   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
12   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
13   !!            3.4  !  2012     (J. Chanut) straight open boundary case update
14   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) optimization of BDY communications
15   !!            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, iiSstdiag, iiSnddiag, iiSstsono, iiSndsono
150      INTEGER  ::   ijRst, ijRnd, ijSst, ijSnd, ijSstdiag, ijSnddiag, ijSstsono, ijSndsono
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    = Njs0            ! Rcv we-side starting point, excluding sw-corner
810                     iiRnd    = nn_hls           ;   ijRnd    = Nje0            ! Rcv we-side   ending point, excluding nw-corner
811                     iiSst    = Nie0-nn_hls+1    ;   ijSst    = Njs0            ! Snd ea-side starting point, excluding se-corner
812                     iiSnd    = Nie0             ;   ijSnd    = Nje0            ! 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    = Nie0+1           ;   ijRst    =  Njs0            ! Rcv ea-side starting point, excluding se-corner
820                     iiRnd    = jpi              ;   ijRnd    =  Nje0            ! Rcv ea-side   ending point, excluding ne-corner
821                     iiSst    = Nis0             ;   ijSst    =  Njs0            ! Snd we-side starting point, excluding sw-corner
822                     iiSnd    = Nis0+nn_hls-1    ;   ijSnd    =  Nje0            ! 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    = Nis0             ;   ijRst    = 1                ! Rcv so-side starting point, excluding sw-corner
832                     iiRnd    = Nie0             ;   ijRnd    = nn_hls           ! Rcv so-side   ending point, excluding se-corner
833                     iiSst    = Nis0             ;   ijSst    = Nje0-nn_hls+1    ! Snd no-side starting point, excluding nw-corner
834                     iiSnd    = Nie0             ;   ijSnd    = Nje0             ! 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    = Nis0             ;   ijRst    = Nje0+1           ! Rcv no-side starting point, excluding nw-corner
843                     iiRnd    = Nie0             ;   ijRnd    = jpj              ! Rcv no-side   ending point, excluding ne-corner
844                     iiSst    = Nis0             ;   ijSst    = Njs0             ! Snd so-side starting point, excluding sw-corner
845                     iiSnd    = Nie0             ;   ijSnd    = Njs0+nn_hls-1    ! 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                     iiRst     = 1               ;   ijRst     = 1               ! Rcv sw-corner starting point
889                     iiRnd     = nn_hls          ;   ijRnd     = nn_hls          ! Rcv sw-corner   ending point
890                     iiSstdiag = Nie0-nn_hls+1   ;   ijSstdiag = Nje0-nn_hls+1   ! send to sw-corner of ne neighbourg
891                     iiSnddiag = Nie0            ;   ijSnddiag = Nje0            ! send to sw-corner of ne neighbourg
892                     iiSstsono = 1               ;   ijSstsono = Nje0-nn_hls+1   ! send to sw-corner of no neighbourg
893                     iiSndsono = nn_hls          ;   ijSndsono = Nje0            ! send to sw-corner of no neighbourg
894                     iioutdir  = -1              ;   ijoutdir  = -1              ! outside MPI domain: westward or southward
895                     !                                          ....|...
896                  CASE( 2 )   ! x: rim on se-corner             :   |o
897                     !          o: potential neighbour(s)     __:__x|o
898                     !             outside of the MPI domain    :o o o
899                     !                                          :   
900                     iRdiag    = jpse            ;   iRsono    = jpso            ! Recv: for se or so
901                     iSdiag    = jpnw            ;   iSsono    = jpno            ! Send: to nw or no
902                     iiRst     = Nie0+1          ;   ijRst     = 1               ! Rcv se-corner starting point
903                     iiRnd     = jpi             ;   ijRnd     = nn_hls          ! Rcv se-corner   ending point
904                     iiSstdiag = Nis0            ;   ijSstdiag = Nje0-nn_hls+1   ! send to se-corner of nw neighbourg
905                     iiSnddiag = Nis0+nn_hls-1   ;   ijSnddiag = Nje0            ! send to se-corner of nw neighbourg
906                     iiSstsono = Nie0+1          ;   ijSstsono = Nje0-nn_hls+1   ! send to se-corner of no neighbourg
907                     iiSndsono = jpi             ;   ijSndsono = Nje0            ! send to se-corner of no neighbourg
908                     iioutdir  = 1               ;   ijoutdir  = -1              ! outside MPI domain: eastward or southward
909                     !                                              :       
910                     !                                         o o_o:___
911                  CASE( 3 )   ! x: rim on nw-corner            o|x  :
912                     !          o: potential neighbour(s)    ..o|...:
913                     !             outside of the MPI domain    |
914                     iRdiag    = jpnw            ;   iRsono    = jpno            ! Recv: for nw or no
915                     iSdiag    = jpse            ;   iSsono    = jpso            ! Send: to se or so
916                     iiRst     = 1               ;   ijRst     = Nje0+1          ! Rcv nw-corner starting point
917                     iiRnd     = nn_hls          ;   ijRnd     = jpj             ! Rcv nw-corner   ending point
918                     iiSstdiag = Nie0-nn_hls+1   ;   ijSstdiag = Njs0            ! send to nw-corner of se neighbourg
919                     iiSnddiag = Nie0            ;   ijSnddiag = Njs0+nn_hls-1   ! send to nw-corner of se neighbourg
920                     iiSstsono = 1               ;   ijSstsono = Njs0            ! send to nw-corner of so neighbourg
921                     iiSndsono = nn_hls          ;   ijSndsono = Njs0+nn_hls-1   ! send to nw-corner of so neighbourg
922                     iioutdir  = -1              ;   ijoutdir  =  1              ! outside MPI domain: westward or northward
923                     !                                          :       
924                     !                                       ___:o_o o
925                  CASE( 4 )   ! x: rim on ne-corner             :  x|o
926                     !          o: potential neighbour(s)       :...|o...
927                     !             outside of the MPI domain        |
928                     iRdiag    = jpne            ;   iRsono    = jpno            ! Recv: for ne or no
929                     iSdiag    = jpsw            ;   iSsono    = jpso            ! Send: to sw or so
930                     iiRst     = Nie0+1          ;   ijRst     = Nje0+1          ! Rcv ne-corner starting point
931                     iiRnd     = jpi             ;   ijRnd     = jpj             ! Rcv ne-corner   ending point
932                     iiSstdiag = Nis0            ;   ijSstdiag = Njs0            ! send to ne-corner of sw neighbourg
933                     iiSnddiag = Nis0+nn_hls-1   ;   ijSnddiag = Njs0+nn_hls-1   ! send to ne-corner of sw neighbourg
934                     iiSstsono = Nie0+1          ;   ijSstsono = Njs0            ! send to ne-corner of so neighbourg
935                     iiSndsono = jpi             ;   ijSndsono = Njs0+nn_hls-1   ! send to ne-corner of so neighbourg
936                     iioutdir  = 1               ;   ijoutdir  = 1               ! outside MPI domain: eastward or southward
937                  END SELECT
938                  !
939                  ! Check if we need to receive data for this rim point
940                  IF( ii >= iiRst .AND. ii <= iiRnd .AND. ij >= ijRst .AND. ij <= ijRnd ) THEN   ! rim point on the corner
941                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the MPI domain?
942                     ! take care of neighbourg(s) in the interior of the computational domain
943                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of the MPI domain
944                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN     ! -> I cannot compute it -> recv it
945                        IF( mpiRnei(nn_hls,iRdiag) > -1                    )   lrecv_bdyint(ib_bdy,igrd,iRdiag,ir) = .TRUE.   ! Receive directly from diagonal neighbourg
946                        IF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 )   lrecv_bdyint(ib_bdy,igrd,iRsono,ir) = .TRUE.   ! Receive through the South/North neighbourg
947                     ENDIF
948                     ! take care of neighbourg in the exterior of the computational domain
949                     IF(  iibe==iiout .OR. ijbe==ijout ) THEN   ! Neib outside of the MPI domain -> I cannot compute it -> recv it
950                        IF( mpiRnei(nn_hls,iRdiag) > -1                    )   lrecv_bdyext(ib_bdy,igrd,iRdiag,ir) = .TRUE.   ! Receive directly from diagonal neighbourg
951                        IF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 )   lrecv_bdyext(ib_bdy,igrd,iRsono,ir) = .TRUE.   ! Receive through the South/North neighbourg
952                     ENDIF
953                  ENDIF
954                  !
955                  ! Check if this rim point corresponds to the corner of one neighbourg. if yes, do we need to send data?
956                  ! Direct send to diag: Is this rim point the corner point of a diag neighbour with which we communicate?
957                  IF( ii >= iiSstdiag .AND. ii <= iiSnddiag .AND. ij >= ijSstdiag .AND. ij <= ijSnddiag   &
958                     &                .AND. mpiSnei(nn_hls,iSdiag) > -1 ) THEN
959                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the nei MPI domain?
960                     ! take care of neighbourg(s) in the interior of the computational domain
961                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of diag nei MPI
962                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout )      &   ! domain -> nei cannot compute it
963                        &    lsend_bdyint(ib_bdy,igrd,iSdiag,ir) = .TRUE.                        ! send rim point data to diag nei
964                     ! take care of neighbourg in the exterior of the computational domain
965                     IF(  iibe==iiout .OR. ijbe==ijout )   &                                 
966                        &    lsend_bdyext(ib_bdy,igrd,iSdiag,ir) = .TRUE.
967                  ENDIF
968                  ! Indirect send to diag (through so/no): rim point is the corner point of a so/no nei with which we communicate
969                  IF( ii >= iiSstsono .AND. ii <= iiSndsono .AND. ij >= ijSstsono .AND. ij <= ijSndsono   &
970                     &                .AND. mpiSnei(nn_hls,iSsono) > -1 .AND. nn_comm == 1 ) THEN
971                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the nei MPI domain?
972                     ! take care of neighbourg(s) in the interior of the computational domain
973                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of so/no nei MPI
974                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout )      &   ! domain -> nei cannot compute it
975                        &    lsend_bdyint(ib_bdy,igrd,iSsono,ir) = .TRUE.                        ! send rim point data to so/no nei
976                     ! take care of neighbourg in the exterior of the computational domain
977                     IF(  iibe==iiout .OR. ijbe==ijout )   &
978                        &    lsend_bdyext(ib_bdy,igrd,iSsono,ir) = .TRUE.
979                  ENDIF
980                  !
981               END DO   ! 4 corners
982            END DO   ! ib
983         END DO   ! igrd
984
985         ! Comment out for debug
986!!$         DO ir = 0,1
987!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing,   &
988!!$               &                              lsend = lsend_bdyint(ib_bdy,1,:,ir), lrecv = lrecv_bdyint(ib_bdy,1,:,ir) )
989!!$            IF(lwp) WRITE(numout,*) ' bdy debug int T', ir ; CALL FLUSH(numout)
990!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing,   &
991!!$               &                              lsend = lsend_bdyint(ib_bdy,2,:,ir), lrecv = lrecv_bdyint(ib_bdy,2,:,ir) )
992!!$            IF(lwp) WRITE(numout,*) ' bdy debug int U', ir ; CALL FLUSH(numout)
993!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing,   &
994!!$               &                              lsend = lsend_bdyint(ib_bdy,3,:,ir), lrecv = lrecv_bdyint(ib_bdy,3,:,ir) )   
995!!$            IF(lwp) WRITE(numout,*) ' bdy debug int V', ir ; CALL FLUSH(numout)
996!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing,   &
997!!$               &                              lsend = lsend_bdyext(ib_bdy,1,:,ir), lrecv = lrecv_bdyext(ib_bdy,1,:,ir) )
998!!$            IF(lwp) WRITE(numout,*) ' bdy debug ext T', ir ; CALL FLUSH(numout)
999!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing,   &
1000!!$               &                              lsend = lsend_bdyext(ib_bdy,2,:,ir), lrecv = lrecv_bdyext(ib_bdy,2,:,ir) )
1001!!$            IF(lwp) WRITE(numout,*) ' bdy debug ext U', ir ; CALL FLUSH(numout)
1002!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing,   &
1003!!$               &                              lsend = lsend_bdyext(ib_bdy,3,:,ir), lrecv = lrecv_bdyext(ib_bdy,3,:,ir) )   
1004!!$            IF(lwp) WRITE(numout,*) ' bdy debug ext V', ir ; CALL FLUSH(numout)
1005!!$         END DO
1006         
1007      END DO   ! ib_bdy
1008
1009      DO ib_bdy = 1,nb_bdy
1010         IF(  cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &
1011            & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &
1012            & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN
1013            DO igrd = 1, jpbgrd
1014               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1015                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1016                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1017                  IF(  mig0(ii) > 2 .AND. mig0(ii) < Ni0glo-2 .AND. mjg0(ij) > 2 .AND. mjg0(ij) < Nj0glo-2  ) THEN
1018                     WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain'
1019                     CALL ctl_stop( ctmp1 )
1020                  END IF
1021               END DO
1022            END DO
1023         END IF
1024      END DO
1025      !
1026      DEALLOCATE( nbidta, nbjdta, nbrdta )
1027      !
1028   END SUBROUTINE bdy_def
1029
1030
1031   SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 )
1032      !!----------------------------------------------------------------------
1033      !!                 ***  ROUTINE bdy_rim_treat  ***
1034      !!
1035      !! ** Purpose :   Initialize structures ( flagu, flagv, ntreat ) indicating how rim points
1036      !!                  are to be handled in the boundary condition treatment
1037      !!
1038      !! ** Method  :   - to handle rim 0 zmasks must indicate ocean points      (set at one on rim 0 and rim 1 and interior)
1039      !!                            and bdymasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
1040      !!                    (as if rim 1 was free ocean)
1041      !!                - to handle rim 1 zmasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
1042      !!                            and bdymasks must indicate free ocean points (set at one on interior)
1043      !!                    (as if rim 0 was land)
1044      !!                - we can then check in which direction the interior of the computational domain is with the difference
1045      !!                         mask array values on both sides to compute flagu and flagv
1046      !!                - and look at the ocean neighbours to compute ntreat
1047      !!----------------------------------------------------------------------
1048      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pumask, pvmask   ! temporary u/v mask array
1049      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pfmask           ! temporary fmask excluding coastal boundary condition (shlat)
1050      LOGICAL                             , INTENT (in   ) :: lrim0            ! .true. -> rim 0   .false. -> rim 1
1051      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, icount       ! dummy loop indices
1052      INTEGER  ::   i_offset, j_offset, inn                ! local integer
1053      INTEGER  ::   ibeg, iend                             ! local integer
1054      LOGICAL  ::   llnon, llson, llean, llwen             ! local logicals indicating the presence of a ocean neighbour
1055      REAL(wp), POINTER, DIMENSION(:,:)       ::   zmask   ! pointer to 2D mask fields
1056      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
1057      CHARACTER(LEN=1), DIMENSION(jpbgrd)     ::   cgrid
1058      REAL(wp)        , DIMENSION(jpi,jpj)    ::   ztmp
1059      !!----------------------------------------------------------------------
1060
1061      cgrid = (/'t','u','v'/)
1062
1063      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
1064
1065         DO igrd = 1, jpbgrd
1066           
1067            IF( lrim0 ) THEN   ! extent of rim 0
1068               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
1069            ELSE               ! extent of rim 1
1070               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
1071            END IF
1072
1073            ! Calculate relationship of U direction to the local orientation of the boundary
1074            ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
1075            ! flagu =  0 : u is tangential
1076            ! flagu =  1 : u is normal to the boundary and is direction is inward
1077            SELECT CASE( igrd )
1078               CASE( 1 )   ;   zmask => pumask     ;   i_offset = 0   ! U(i-1)   T(i)   U(i  )
1079               CASE( 2 )   ;   zmask => bdytmask   ;   i_offset = 1   ! T(i  )   U(i)   T(i+1)
1080               CASE( 3 )   ;   zmask => pfmask     ;   i_offset = 0   ! F(i-1)   V(i)   F(i  )
1081            END SELECT
1082            icount = 0
1083            ztmp(:,:) = -999._wp
1084            DO ib = ibeg, iend 
1085               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1086               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1087               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
1088               zwfl = zmask(ii+i_offset-1,ij)
1089               zefl = zmask(ii+i_offset  ,ij)
1090               ! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
1091               IF( i_offset == 1 .and. zefl + zwfl == 2._wp ) THEN
1092                  icount = icount + 1
1093                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
1094               ELSE
1095                  ztmp(ii,ij) = -zwfl + zefl
1096               ENDIF
1097            END DO
1098            IF( icount /= 0 ) THEN
1099               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
1100                  ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1101               CALL ctl_stop( ctmp1 )
1102            ENDIF
1103            SELECT CASE( igrd )
1104               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1105               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1106               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1107            END SELECT
1108            DO ib = ibeg, iend
1109               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1110               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1111               idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij)
1112            END DO
1113           
1114            ! Calculate relationship of V direction to the local orientation of the boundary
1115            ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
1116            ! flagv =  0 : v is tangential
1117            ! flagv =  1 : v is normal to the boundary and is direction is inward
1118            SELECT CASE( igrd )
1119               CASE( 1 )   ;   zmask => pvmask     ;   j_offset = 0
1120               CASE( 2 )   ;   zmask => pfmask     ;   j_offset = 0
1121               CASE( 3 )   ;   zmask => bdytmask   ;   j_offset = 1
1122            END SELECT
1123            icount = 0
1124            ztmp(:,:) = -999._wp
1125            DO ib = ibeg, iend
1126               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1127               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1128               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
1129               zsfl = zmask(ii,ij+j_offset-1)
1130               znfl = zmask(ii,ij+j_offset  )
1131               ! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
1132               IF( j_offset == 1 .and. znfl + zsfl == 2._wp ) THEN
1133                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
1134                  icount = icount + 1
1135               ELSE
1136                  ztmp(ii,ij) = -zsfl + znfl
1137               END IF
1138            END DO
1139            IF( icount /= 0 ) THEN
1140               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
1141                  ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1142               CALL ctl_stop( ctmp1 )
1143            ENDIF
1144            SELECT CASE( igrd )
1145               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1146               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1147               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1148            END SELECT
1149            DO ib = ibeg, iend
1150               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1151               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1152               idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij)
1153            END DO
1154     
1155            ! Calculate ntreat
1156            SELECT CASE( igrd )
1157               CASE( 1 )   ;   zmask => bdytmask 
1158               CASE( 2 )   ;   zmask => bdyumask 
1159               CASE( 3 )   ;   zmask => bdyvmask 
1160            END SELECT
1161            ztmp(:,:) = -999._wp
1162            DO ib = ibeg, iend
1163               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1164               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1165               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
1166               llnon = zmask(ii  ,ij+1) == 1._wp 
1167               llson = zmask(ii  ,ij-1) == 1._wp 
1168               llean = zmask(ii+1,ij  ) == 1._wp 
1169               llwen = zmask(ii-1,ij  ) == 1._wp 
1170               inn  = COUNT( (/ llnon, llson, llean, llwen /) )
1171               IF( inn == 0 ) THEN   ! no neighbours -> interior of a corner  or  cluster of rim points
1172                  !               !              !     _____     !     _____    !    __     __
1173                  !  1 |   o      !  2  o   |    !  3 | x        !  4     x |   !      |   |   -> error
1174                  !    |_x_ _     !    _ _x_|    !    |   o      !      o   |   !      |x_x|
1175                  IF(     zmask(ii+1,ij+1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 1._wp
1176                  ELSEIF( zmask(ii-1,ij+1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 2._wp
1177                  ELSEIF( zmask(ii+1,ij-1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 3._wp
1178                  ELSEIF( zmask(ii-1,ij-1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 4._wp
1179                  ELSE                                       ;   ztmp(ii,ij) = -1._wp
1180                     WRITE(ctmp1,*) 'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
1181                       ' on boundary set ', ib_bdy, ' has no free ocean neighbour'
1182                     IF( lrim0 ) THEN
1183                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 0 points.'
1184                     ELSE
1185                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.'
1186                     END IF
1187                     CALL ctl_warn( ctmp1, ctmp2 )
1188                  END IF
1189               END IF
1190               IF( inn == 1 ) THEN   ! middle of linear bdy  or incomplete corner  ! ___ o
1191                  !    |         !         |   !      o     !    ______            !    |x___
1192                  ! 5  | x o     ! 6   o x |   ! 7  __x__   ! 8    x
1193                  !    |         !         |   !            !      o
1194                  IF( llean )   ztmp(ii,ij) = 5._wp
1195                  IF( llwen )   ztmp(ii,ij) = 6._wp
1196                  IF( llnon )   ztmp(ii,ij) = 7._wp
1197                  IF( llson )   ztmp(ii,ij) = 8._wp
1198               END IF
1199               IF( inn == 2 ) THEN   ! exterior of a corner
1200                  !        o      !        o      !    _____|       !       |_____ 
1201                  !  9 ____x o    ! 10   o x___   ! 11     x o      ! 12   o x     
1202                  !         |     !       |       !        o        !        o
1203                  IF( llnon .AND. llean )   ztmp(ii,ij) =  9._wp
1204                  IF( llnon .AND. llwen )   ztmp(ii,ij) = 10._wp
1205                  IF( llson .AND. llean )   ztmp(ii,ij) = 11._wp
1206                  IF( llson .AND. llwen )   ztmp(ii,ij) = 12._wp
1207               END IF
1208               IF( inn == 3 ) THEN   ! 3 neighbours     __   __
1209                  !    |_  o      !        o  _|  !       |_|     !       o         
1210                  ! 13  _| x o    ! 14   o x |_   ! 15   o x o    ! 16  o x o       
1211                  !    |   o      !        o   |  !        o      !    __|¨|__   
1212                  IF( llnon .AND. llean .AND. llson )   ztmp(ii,ij) = 13._wp
1213                  IF( llnon .AND. llwen .AND. llson )   ztmp(ii,ij) = 14._wp
1214                  IF( llwen .AND. llson .AND. llean )   ztmp(ii,ij) = 15._wp
1215                  IF( llwen .AND. llnon .AND. llean )   ztmp(ii,ij) = 16._wp
1216               END IF
1217               IF( inn == 4 ) THEN
1218                  WRITE(ctmp1,*)  'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
1219                       ' on boundary set ', ib_bdy, ' have 4 neighbours'
1220                  CALL ctl_stop( ctmp1 )
1221               END IF
1222            END DO
1223            SELECT CASE( igrd )
1224               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1225               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1226               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1227            END SELECT
1228            DO ib = ibeg, iend
1229               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1230               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1231               idx_bdy(ib_bdy)%ntreat(ib,igrd) = NINT(ztmp(ii,ij))
1232            END DO
1233            !
1234         END DO   ! jpbgrd
1235         !
1236      END DO   ! ib_bdy
1237
1238    END SUBROUTINE bdy_rim_treat
1239
1240   
1241    SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )
1242      !!----------------------------------------------------------------------
1243      !!                 ***  ROUTINE find_neib  ***
1244      !!
1245      !! ** Purpose :   get ii1, ij1, ii2, ij2, ii3, ij3, the indices of
1246      !!               the free ocean neighbours of (ii,ij) for bdy treatment
1247      !!
1248      !! ** Method  :  use itreat input to select a case
1249      !!               N.B. ntreat is defined for all bdy points in routine bdy_rim_treat
1250      !!
1251      !!----------------------------------------------------------------------
1252      INTEGER, INTENT(in   )      ::   ii, ij, itreat
1253      INTEGER, INTENT(  out)      ::   ii1, ij1, ii2, ij2, ii3, ij3
1254      !!----------------------------------------------------------------------
1255      SELECT CASE( itreat )   ! points that will be used by bdy routines, -1 will be discarded
1256         !               !               !     _____     !     _____     
1257         !  1 |   o      !  2  o   |     !  3 | x        !  4     x |   
1258         !    |_x_ _     !    _ _x_|     !    |   o      !      o   |
1259      CASE( 1 )    ;   ii1 = ii+1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1260      CASE( 2 )    ;   ii1 = ii-1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1261      CASE( 3 )    ;   ii1 = ii+1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1262      CASE( 4 )    ;   ii1 = ii-1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1263         !    |          !         |     !      o        !    ______                   ! or incomplete corner
1264         ! 5  | x o      ! 6   o x |     ! 7  __x__      ! 8    x                      !  7  ____ o
1265         !    |          !         |     !               !      o                      !         |x___
1266      CASE( 5 )    ;   ii1 = ii+1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1267      CASE( 6 )    ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1268      CASE( 7 )    ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1269      CASE( 8 )    ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1270         !        o      !        o      !    _____|     !       |_____ 
1271         !  9 ____x o    ! 10   o x___   ! 11     x o    ! 12   o x     
1272         !         |     !       |       !        o      !        o     
1273      CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
1274      CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1275      CASE( 11 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1276      CASE( 12 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1277         !    |_  o      !        o  _|  !     ¨¨|_|¨¨   !       o         
1278         ! 13  _| x o    !  14  o x |_   !  15  o x o    ! 16  o x o       
1279         !    |   o      !        o   |  !        o      !    __|¨|__
1280      CASE( 13 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1   
1281      CASE( 14 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1 
1282      CASE( 15 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij   
1283      CASE( 16 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij
1284      END SELECT
1285   END SUBROUTINE find_neib
1286   
1287
1288   SUBROUTINE bdy_read_seg( kb_bdy, knblendta ) 
1289      !!----------------------------------------------------------------------
1290      !!                 ***  ROUTINE bdy_read_seg  ***
1291      !!
1292      !! ** Purpose :  build bdy coordinates with segments defined in namelist
1293      !!
1294      !! ** Method  :  read namelist nambdy_index blocks
1295      !!
1296      !!----------------------------------------------------------------------
1297      INTEGER                   , INTENT (in   ) ::   kb_bdy           ! bdy number
1298      INTEGER, DIMENSION(jpbgrd), INTENT (  out) ::   knblendta        ! length of index arrays
1299      !!
1300      INTEGER          ::   ios                 ! Local integer output status for namelist read
1301      INTEGER          ::   nbdyind, nbdybeg, nbdyend
1302      INTEGER          ::   nbdy_count, nbdy_rdstart, nbdy_loc
1303      CHARACTER(LEN=1) ::   ctypebdy   !     -        -
1304      CHARACTER(LEN=50)::   cerrmsg    !     -        -
1305      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
1306      !!----------------------------------------------------------------------
1307      ! Need to support possibility of reading more than one nambdy_index from
1308      ! the namelist_cfg internal file.
1309      ! Do this by finding the kb_bdy'th occurence of nambdy_index in the
1310      ! character buffer as the starting point.
1311      nbdy_rdstart = 1
1312      DO nbdy_count = 1, kb_bdy
1313       nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_index' )
1314       IF( nbdy_loc .GT. 0 ) THEN
1315          nbdy_rdstart = nbdy_rdstart + nbdy_loc
1316       ELSE
1317          WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',kb_bdy,' of nambdy_index not found'
1318          ios = -1
1319          CALL ctl_nam ( ios , cerrmsg )
1320       ENDIF
1321      END DO
1322      nbdy_rdstart = MAX( 1, nbdy_rdstart - 2 )
1323      READ  ( numnam_cfg( nbdy_rdstart: ), nambdy_index, IOSTAT = ios, ERR = 904)
1324904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' )
1325      IF(lwm) WRITE ( numond, nambdy_index )
1326     
1327      SELECT CASE ( TRIM(ctypebdy) )
1328      CASE( 'N' )
1329         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1330            nbdyind  = Nj0glo - 2  ! set boundary to whole side of model domain.
1331            nbdybeg  = 2
1332            nbdyend  = Ni0glo - 1
1333         ENDIF
1334         nbdysegn = nbdysegn + 1
1335         npckgn(nbdysegn) = kb_bdy ! Save bdy package number
1336         jpjnob(nbdysegn) = nbdyind 
1337         jpindt(nbdysegn) = nbdybeg
1338         jpinft(nbdysegn) = nbdyend
1339         !
1340      CASE( 'S' )
1341         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1342            nbdyind  = 2           ! set boundary to whole side of model domain.
1343            nbdybeg  = 2
1344            nbdyend  = Ni0glo - 1
1345         ENDIF
1346         nbdysegs = nbdysegs + 1
1347         npckgs(nbdysegs) = kb_bdy ! Save bdy package number
1348         jpjsob(nbdysegs) = nbdyind
1349         jpisdt(nbdysegs) = nbdybeg
1350         jpisft(nbdysegs) = nbdyend
1351         !
1352      CASE( 'E' )
1353         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1354            nbdyind  = Ni0glo - 2  ! set boundary to whole side of model domain.
1355            nbdybeg  = 2
1356            nbdyend  = Nj0glo - 1
1357         ENDIF
1358         nbdysege = nbdysege + 1 
1359         npckge(nbdysege) = kb_bdy ! Save bdy package number
1360         jpieob(nbdysege) = nbdyind
1361         jpjedt(nbdysege) = nbdybeg
1362         jpjeft(nbdysege) = nbdyend
1363         !
1364      CASE( 'W' )
1365         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1366            nbdyind  = 2           ! set boundary to whole side of model domain.
1367            nbdybeg  = 2
1368            nbdyend  = Nj0glo - 1
1369         ENDIF
1370         nbdysegw = nbdysegw + 1
1371         npckgw(nbdysegw) = kb_bdy ! Save bdy package number
1372         jpiwob(nbdysegw) = nbdyind
1373         jpjwdt(nbdysegw) = nbdybeg
1374         jpjwft(nbdysegw) = nbdyend
1375         !
1376      CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
1377      END SELECT
1378     
1379      ! For simplicity we assume that in case of straight bdy, arrays have the same length
1380      ! (even if it is true that last tangential velocity points
1381      ! are useless). This simplifies a little bit boundary data format (and agrees with format
1382      ! used so far in obc package)
1383     
1384      knblendta(1:jpbgrd) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy)
1385     
1386   END SUBROUTINE bdy_read_seg
1387
1388   
1389   SUBROUTINE bdy_ctl_seg
1390      !!----------------------------------------------------------------------
1391      !!                 ***  ROUTINE bdy_ctl_seg  ***
1392      !!
1393      !! ** Purpose :   Check straight open boundary segments location
1394      !!
1395      !! ** Method  :   - Look for open boundary corners
1396      !!                - Check that segments start or end on land
1397      !!----------------------------------------------------------------------
1398      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest 
1399      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns 
1400      REAL(wp), DIMENSION(2) ::   ztestmask
1401      !!----------------------------------------------------------------------
1402      !
1403      IF (lwp) WRITE(numout,*) ' '
1404      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1405      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1406      !
1407      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege
1408      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw
1409      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn
1410      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs
1411      !
1412      ! 1. Check bounds
1413      !----------------
1414      DO ib = 1, nbdysegn
1415         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
1416         IF ((jpjnob(ib).ge.Nj0glo-1).or.& 
1417            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1418         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1419         IF (jpindt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1420         IF (jpinft(ib).gt.Ni0glo)     CALL ctl_stop( 'End index out of domain' )
1421      END DO
1422      !
1423      DO ib = 1, nbdysegs
1424         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
1425         IF ((jpjsob(ib).ge.Nj0glo-1).or.& 
1426            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1427         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1428         IF (jpisdt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1429         IF (jpisft(ib).gt.Ni0glo)     CALL ctl_stop( 'End index out of domain' )
1430      END DO
1431      !
1432      DO ib = 1, nbdysege
1433         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib)
1434         IF ((jpieob(ib).ge.Ni0glo-1).or.& 
1435            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1436         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1437         IF (jpjedt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1438         IF (jpjeft(ib).gt.Nj0glo)     CALL ctl_stop( 'End index out of domain' )
1439      END DO
1440      !
1441      DO ib = 1, nbdysegw
1442         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib)
1443         IF ((jpiwob(ib).ge.Ni0glo-1).or.& 
1444            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1445         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1446         IF (jpjwdt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1447         IF (jpjwft(ib).gt.Nj0glo)     CALL ctl_stop( 'End index out of domain' )
1448      ENDDO
1449      !     
1450      ! 2. Look for segment crossings
1451      !------------------------------
1452      IF (lwp) WRITE(numout,*) '**Look for segments corners  :'
1453      !
1454      itest = 0 ! corner number
1455      !
1456      ! flag to detect if start or end of open boundary belongs to a corner
1457      ! if not (=0), it must be on land.
1458      ! if a corner is detected, save bdy package number for further tests
1459      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1460      ! South/West crossings
1461      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1462         DO ib1 = 1, nbdysegw       
1463            DO ib2 = 1, nbdysegs
1464               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1465                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1466                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1467                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1468                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN 
1469                     ! We have a possible South-West corner                     
1470!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1471!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1472                     icornw(ib1,1) = npckgs(ib2)
1473                     icorns(ib2,1) = npckgw(ib1)
1474                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1475                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1476                        &                                     jpisft(ib2), jpjwft(ib1)
1477                     WRITE(ctmp2,*) ' Not allowed yet'
1478                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 
1479                        &                            ' and South segment: ',npckgs(ib2)
1480                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1481                  ELSE
1482                     WRITE(ctmp1,*) ' Check South and West Open boundary indices'
1483                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , &
1484                        &                            ' and South segment: ',npckgs(ib2)
1485                     CALL ctl_stop( ctmp1, ctmp2 )
1486                  END IF
1487               END IF
1488            END DO
1489         END DO
1490      END IF
1491      !
1492      ! South/East crossings
1493      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1494         DO ib1 = 1, nbdysege
1495            DO ib2 = 1, nbdysegs
1496               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1497                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1498                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. &
1499                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN
1500                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1501                     ! We have a possible South-East corner
1502!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1503!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1504                     icorne(ib1,1) = npckgs(ib2)
1505                     icorns(ib2,2) = npckge(ib1)
1506                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1507                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1508                        &                                     jpisdt(ib2), jpjeft(ib1)
1509                     WRITE(ctmp2,*) ' Not allowed yet'
1510                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1511                        &                            ' and South segment: ',npckgs(ib2)
1512                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1513                  ELSE
1514                     WRITE(ctmp1,*) ' Check South and East Open boundary indices'
1515                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1516                     &                               ' and South segment: ',npckgs(ib2)
1517                     CALL ctl_stop( ctmp1, ctmp2 )
1518                  END IF
1519               END IF
1520            END DO
1521         END DO
1522      END IF
1523      !
1524      ! North/West crossings
1525      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1526         DO ib1 = 1, nbdysegw       
1527            DO ib2 = 1, nbdysegn
1528               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. &
1529                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. &
1530                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1531                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1532                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1533                     ! We have a possible North-West corner
1534!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1535!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1536                     icornw(ib1,2) = npckgn(ib2)
1537                     icornn(ib2,1) = npckgw(ib1)
1538                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1539                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1540                     &                                     jpinft(ib2), jpjwdt(ib1)
1541                     WRITE(ctmp2,*) ' Not allowed yet'
1542                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1543                     &                               ' and North segment: ',npckgn(ib2)
1544                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1545                  ELSE
1546                     WRITE(ctmp1,*) ' Check North and West Open boundary indices'
1547                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1548                     &                               ' and North segment: ',npckgn(ib2)
1549                     CALL ctl_stop( ctmp1, ctmp2 )
1550                  END IF
1551               END IF
1552            END DO
1553         END DO
1554      END IF
1555      !
1556      ! North/East crossings
1557      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1558         DO ib1 = 1, nbdysege       
1559            DO ib2 = 1, nbdysegn
1560               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1561                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1562                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1563                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1564                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1565                     ! We have a possible North-East corner
1566!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1567!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1568                     icorne(ib1,2) = npckgn(ib2)
1569                     icornn(ib2,2) = npckge(ib1)
1570                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1571                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1572                     &                                     jpindt(ib2), jpjedt(ib1)
1573                     WRITE(ctmp2,*) ' Not allowed yet'
1574                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1575                     &                               ' and North segment: ',npckgn(ib2)
1576                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1577                  ELSE
1578                     WRITE(ctmp1,*) ' Check North and East Open boundary indices'
1579                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1580                     &                               ' and North segment: ',npckgn(ib2)
1581                     CALL ctl_stop( ctmp1, ctmp2 )
1582                  END IF
1583               END IF
1584            END DO
1585         END DO
1586      END IF
1587      !
1588      ! 3. Check if segment extremities are on land
1589      !--------------------------------------------
1590      !
1591      ! West segments
1592      DO ib = 1, nbdysegw
1593         ! get mask at boundary extremities:
1594         ztestmask(1:2)=0.
1595         DO ji = 1, jpi
1596            DO jj = 1, jpj             
1597              IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1598              IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1599            END DO
1600         END DO
1601         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1602
1603         IF (ztestmask(1)==1) THEN
1604            IF (icornw(ib,1)==0) THEN
1605               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1606               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1607            ELSE
1608               ! This is a corner
1609               IF(lwp) WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1610               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1611               itest=itest+1
1612            ENDIF
1613         ENDIF
1614         IF (ztestmask(2)==1) THEN
1615            IF (icornw(ib,2)==0) THEN
1616               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1617               CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' )
1618            ELSE
1619               ! This is a corner
1620               IF(lwp) WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1621               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1622               itest=itest+1
1623            ENDIF
1624         ENDIF
1625      END DO
1626      !
1627      ! East segments
1628      DO ib = 1, nbdysege
1629         ! get mask at boundary extremities:
1630         ztestmask(1:2)=0.
1631         DO ji = 1, jpi
1632            DO jj = 1, jpj             
1633              IF( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjedt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1634              IF( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjeft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1635            END DO
1636         END DO
1637         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1638
1639         IF (ztestmask(1)==1) THEN
1640            IF (icorne(ib,1)==0) THEN
1641               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1642               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1643            ELSE
1644               ! This is a corner
1645               IF(lwp) WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1646               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1647               itest=itest+1
1648            ENDIF
1649         ENDIF
1650         IF (ztestmask(2)==1) THEN
1651            IF (icorne(ib,2)==0) THEN
1652               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1653               CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1654            ELSE
1655               ! This is a corner
1656               IF(lwp) WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1657               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1658               itest=itest+1
1659            ENDIF
1660         ENDIF
1661      END DO
1662      !
1663      ! South segments
1664      DO ib = 1, nbdysegs
1665         ! get mask at boundary extremities:
1666         ztestmask(1:2)=0.
1667         DO ji = 1, jpi
1668            DO jj = 1, jpj             
1669              IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1670              IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1671            END DO
1672         END DO
1673         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1674
1675         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1676            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1677            CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1678         ENDIF
1679         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1680            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1681            CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1682         ENDIF
1683      END DO
1684      !
1685      ! North segments
1686      DO ib = 1, nbdysegn
1687         ! get mask at boundary extremities:
1688         ztestmask(1:2)=0.
1689         DO ji = 1, jpi
1690            DO jj = 1, jpj             
1691               IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpindt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1692               IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpinft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1693            END DO
1694         END DO
1695         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1696
1697         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1698            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1699            CALL ctl_stop( ctmp1, ' does not start on land' )
1700         ENDIF
1701         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1702            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1703            CALL ctl_stop( ctmp1, ' does not end on land' )
1704         ENDIF
1705      END DO
1706      !
1707      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1708      !
1709      ! Other tests TBD:
1710      ! segments completly on land
1711      ! optimized open boundary array length according to landmask
1712      ! Nudging layers that overlap with interior domain
1713      !
1714   END SUBROUTINE bdy_ctl_seg
1715
1716   
1717   SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta ) 
1718      !!----------------------------------------------------------------------
1719      !!                 ***  ROUTINE bdy_coords_seg  ***
1720      !!
1721      !! ** Purpose :  build nbidta, nbidta, nbrdta for bdy built with segments
1722      !!
1723      !! ** Method  : 
1724      !!
1725      !!----------------------------------------------------------------------
1726      INTEGER, DIMENSION(:,:,:), intent(  out)  ::   nbidta, nbjdta, nbrdta   ! Index arrays: i and j indices of bdy dta
1727      !!
1728      INTEGER  ::   ii, ij, ir, iseg
1729      INTEGER  ::   igrd         ! grid type (t=1, u=2, v=3)
1730      INTEGER  ::   icount       !
1731      INTEGER  ::   ib_bdy       ! bdy number
1732      !!----------------------------------------------------------------------
1733
1734      ! East
1735      !-----
1736      DO iseg = 1, nbdysege
1737         ib_bdy = npckge(iseg)
1738         !
1739         ! ------------ T points -------------
1740         igrd=1
1741         icount=0
1742         DO ir = 1, nn_rimwidth(ib_bdy)
1743            DO ij = jpjedt(iseg), jpjeft(iseg)
1744               icount = icount + 1
1745               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir + nn_hls
1746               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1747               nbrdta(icount, igrd, ib_bdy) = ir
1748            ENDDO
1749         ENDDO
1750         !
1751         ! ------------ U points -------------
1752         igrd=2
1753         icount=0
1754         DO ir = 1, nn_rimwidth(ib_bdy)
1755            DO ij = jpjedt(iseg), jpjeft(iseg)
1756               icount = icount + 1
1757               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir + nn_hls
1758               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1759               nbrdta(icount, igrd, ib_bdy) = ir
1760            ENDDO
1761         ENDDO
1762         !
1763         ! ------------ V points -------------
1764         igrd=3
1765         icount=0
1766         DO ir = 1, nn_rimwidth(ib_bdy)
1767            !            DO ij = jpjedt(iseg), jpjeft(iseg) - 1
1768            DO ij = jpjedt(iseg), jpjeft(iseg)
1769               icount = icount + 1
1770               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir + nn_hls
1771               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1772               nbrdta(icount, igrd, ib_bdy) = ir
1773            ENDDO
1774            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1775            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1776         ENDDO
1777      ENDDO
1778      !
1779      ! West
1780      !-----
1781      DO iseg = 1, nbdysegw
1782         ib_bdy = npckgw(iseg)
1783         !
1784         ! ------------ T points -------------
1785         igrd=1
1786         icount=0
1787         DO ir = 1, nn_rimwidth(ib_bdy)
1788            DO ij = jpjwdt(iseg), jpjwft(iseg)
1789               icount = icount + 1
1790               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1791               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1792               nbrdta(icount, igrd, ib_bdy) = ir
1793            ENDDO
1794         ENDDO
1795         !
1796         ! ------------ U points -------------
1797         igrd=2
1798         icount=0
1799         DO ir = 1, nn_rimwidth(ib_bdy)
1800            DO ij = jpjwdt(iseg), jpjwft(iseg)
1801               icount = icount + 1
1802               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1803               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1804               nbrdta(icount, igrd, ib_bdy) = ir
1805            ENDDO
1806         ENDDO
1807         !
1808         ! ------------ V points -------------
1809         igrd=3
1810         icount=0
1811         DO ir = 1, nn_rimwidth(ib_bdy)
1812            !            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
1813            DO ij = jpjwdt(iseg), jpjwft(iseg)
1814               icount = icount + 1
1815               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1816               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1817               nbrdta(icount, igrd, ib_bdy) = ir
1818            ENDDO
1819            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1820            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1821         ENDDO
1822      ENDDO
1823      !
1824      ! North
1825      !-----
1826      DO iseg = 1, nbdysegn
1827         ib_bdy = npckgn(iseg)
1828         !
1829         ! ------------ T points -------------
1830         igrd=1
1831         icount=0
1832         DO ir = 1, nn_rimwidth(ib_bdy)
1833            DO ii = jpindt(iseg), jpinft(iseg)
1834               icount = icount + 1
1835               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1836               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir + nn_hls 
1837               nbrdta(icount, igrd, ib_bdy) = ir
1838            ENDDO
1839         ENDDO
1840         !
1841         ! ------------ U points -------------
1842         igrd=2
1843         icount=0
1844         DO ir = 1, nn_rimwidth(ib_bdy)
1845            !            DO ii = jpindt(iseg), jpinft(iseg) - 1
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) + 2 - ir + nn_hls
1850               nbrdta(icount, igrd, ib_bdy) = ir
1851            ENDDO
1852            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1853            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1854         ENDDO
1855         !
1856         ! ------------ V points -------------
1857         igrd=3
1858         icount=0
1859         DO ir = 1, nn_rimwidth(ib_bdy)
1860            DO ii = jpindt(iseg), jpinft(iseg)
1861               icount = icount + 1
1862               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1863               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir + nn_hls
1864               nbrdta(icount, igrd, ib_bdy) = ir
1865            ENDDO
1866         ENDDO
1867      ENDDO
1868      !
1869      ! South
1870      !-----
1871      DO iseg = 1, nbdysegs
1872         ib_bdy = npckgs(iseg)
1873         !
1874         ! ------------ T points -------------
1875         igrd=1
1876         icount=0
1877         DO ir = 1, nn_rimwidth(ib_bdy)
1878            DO ii = jpisdt(iseg), jpisft(iseg)
1879               icount = icount + 1
1880               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1881               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1882               nbrdta(icount, igrd, ib_bdy) = ir
1883            ENDDO
1884         ENDDO
1885         !
1886         ! ------------ U points -------------
1887         igrd=2
1888         icount=0
1889         DO ir = 1, nn_rimwidth(ib_bdy)
1890            !            DO ii = jpisdt(iseg), jpisft(iseg) - 1
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            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1898            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1899         ENDDO
1900         !
1901         ! ------------ V points -------------
1902         igrd=3
1903         icount=0
1904         DO ir = 1, nn_rimwidth(ib_bdy)
1905            DO ii = jpisdt(iseg), jpisft(iseg)
1906               icount = icount + 1
1907               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1908               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1909               nbrdta(icount, igrd, ib_bdy) = ir
1910            ENDDO
1911         ENDDO
1912      ENDDO
1913
1914     
1915   END SUBROUTINE bdy_coords_seg
1916   
1917   
1918   SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1919      !!----------------------------------------------------------------------
1920      !!                 ***  ROUTINE bdy_ctl_corn  ***
1921      !!
1922      !! ** Purpose :   Check numerical schemes consistency between
1923      !!                segments having a common corner
1924      !!
1925      !! ** Method  :   
1926      !!----------------------------------------------------------------------
1927      INTEGER, INTENT(in)  ::   ib1, ib2
1928      INTEGER :: itest
1929      !!----------------------------------------------------------------------
1930      itest = 0
1931
1932      IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) )   itest = itest + 1
1933      IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) )   itest = itest + 1
1934      IF( cn_tra  (ib1) /= cn_tra  (ib2) )   itest = itest + 1
1935      !
1936      IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) )   itest = itest + 1
1937      IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) )   itest = itest + 1
1938      IF( nn_tra_dta  (ib1) /= nn_tra_dta  (ib2) )   itest = itest + 1
1939      !
1940      IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) )   itest = itest + 1   
1941      !
1942      IF( itest>0 ) THEN
1943         WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2
1944         CALL ctl_stop( ctmp1, ' have different open bdy schemes' )
1945      ENDIF
1946      !
1947   END SUBROUTINE bdy_ctl_corn
1948
1949
1950   SUBROUTINE bdy_meshwri()
1951      !!----------------------------------------------------------------------
1952      !!                 ***  ROUTINE bdy_meshwri  ***
1953      !!         
1954      !! ** Purpose :   write netcdf file with nbr, flagu, flagv, ntreat for T, U
1955      !!                and V points in 2D arrays for easier visualisation/control
1956      !!
1957      !! ** Method  :   use iom_rstput as in domwri.F
1958      !!----------------------------------------------------------------------     
1959      INTEGER  ::   ib_bdy, ii, ij, igrd, ib     ! dummy loop indices
1960      INTEGER  ::   inum                                   !   -       -
1961      REAL(wp), POINTER, DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields
1962      REAL(wp)         , DIMENSION(jpi,jpj) ::   ztmp
1963      CHARACTER(LEN=1) , DIMENSION(jpbgrd)  ::   cgrid
1964      !!----------------------------------------------------------------------     
1965      cgrid = (/'t','u','v'/)
1966      CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. )
1967      DO igrd = 1, jpbgrd
1968         SELECT CASE( igrd )
1969         CASE( 1 )   ;   zmask => tmask(:,:,1)
1970         CASE( 2 )   ;   zmask => umask(:,:,1)
1971         CASE( 3 )   ;   zmask => vmask(:,:,1)
1972         END SELECT
1973         ztmp(:,:) = zmask(:,:)
1974         DO ib_bdy = 1, nb_bdy
1975            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)      ! nbr deined for all rims
1976               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1977               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1978               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10.
1979               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1980            END DO
1981         END DO
1982         CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1983         ztmp(:,:) = zmask(:,:)
1984         DO ib_bdy = 1, nb_bdy
1985            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagu defined only for rims 0 and 1
1986               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1987               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1988               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10.
1989               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1990            END DO
1991         END DO
1992         CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1993         ztmp(:,:) = zmask(:,:)
1994         DO ib_bdy = 1, nb_bdy
1995            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagv defined only for rims 0 and 1
1996               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1997               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1998               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10.
1999               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
2000            END DO
2001         END DO
2002         CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
2003         ztmp(:,:) = zmask(:,:)
2004         DO ib_bdy = 1, nb_bdy
2005            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! ntreat defined only for rims 0 and 1
2006               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
2007               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
2008               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10.
2009               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
2010            END DO
2011         END DO
2012         CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
2013      END DO
2014      CALL iom_close( inum )
2015
2016   END SUBROUTINE bdy_meshwri
2017   
2018   !!=================================================================================
2019END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.