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

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

trunk: BDY compliant with corner communications, see #2731

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