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

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

trunk: bugfix following [15354], see #2731

  • Property svn:keywords set to Id
File size: 111.7 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         .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! so inner side
596                  IF( mpiSnei(nn_hls,jpso) > -1 ) THEN   ;   lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.   ! send to so neighbourg
597                  ELSE                                   ;   CALL ctl_stop( 'bdyini send olr so-side' )
598                  ENDIF
599               ENDIF
600               IF( ii  < Nis0                          .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! so side we-halo
601                  IF( nn_comm == 1 .AND. mpiSnei(nn_hls,jpso) > -1 )   lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
602               ENDIF
603               IF( ii  > Nie0                          .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! so side ea-halo
604                  IF( nn_comm == 1 .AND. mpiSnei(nn_hls,jpso) > -1 )   lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
605               ENDIF
606               IF( ii >= Nis0 .AND. ii <= Nie0         .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! no inner side
607                  IF( mpiSnei(nn_hls,jpno) > -1 ) THEN   ;   lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.   ! send to no neighbourg
608                  ELSE                                   ;   CALL ctl_stop( 'bdyini send olr no-side' )
609                  ENDIF
610               ENDIF
611               IF( ii  < Nis0                          .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! no side we-halo
612                  IF( nn_comm == 1 .AND. mpiSnei(nn_hls,jpno) > -1 )   lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
613               ENDIF
614               IF( ii  > Nie0                          .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! no side ea-halo
615                  IF( nn_comm == 1 .AND. mpiSnei(nn_hls,jpno) > -1 )   lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
616               ENDIF
617               IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! sw inner corner
618                  IF( mpiSnei(nn_hls,jpsw) > -1 ) THEN   ;   lsend_bdyolr(ib_bdy,igrd,jpsw,ir) = .TRUE.   ! send to sw neighbourg
619                  ELSEIF( nn_comm /= 1          ) THEN   ;   CALL ctl_stop( 'bdyini send olr sw-corner' )
620                  ENDIF
621               ENDIF
622               IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN   ! se inner corner
623                  IF( mpiSnei(nn_hls,jpse) > -1 ) THEN   ;   lsend_bdyolr(ib_bdy,igrd,jpse,ir) = .TRUE.   ! send to se neighbourg
624                  ELSEIF( nn_comm /= 1          ) THEN   ;   CALL ctl_stop( 'bdyini send olr se-corner' )
625                  ENDIF
626               ENDIF
627               IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! nw inner corner
628                  IF( mpiSnei(nn_hls,jpnw) > -1 ) THEN   ;   lsend_bdyolr(ib_bdy,igrd,jpnw,ir) = .TRUE.   ! send to nw neighbourg
629                  ELSEIF( nn_comm /= 1          ) THEN   ;   CALL ctl_stop( 'bdyini send olr nw-corner' )
630                  ENDIF
631               ENDIF
632               IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN   ! ne inner corner
633                  IF( mpiSnei(nn_hls,jpne) > -1 ) THEN   ;   lsend_bdyolr(ib_bdy,igrd,jpne,ir) = .TRUE.   ! send to ne neighbourg
634                  ELSEIF( nn_comm /= 1          ) THEN   ;   CALL ctl_stop( 'bdyini send olr ne-corner' )
635                  ENDIF
636               ENDIF
637               !
638               ! check if point has to be received from a neighbour
639               IF( ii  < Nis0                  .AND. ij >= Njs0 .AND. ij <= Nje0 ) THEN   ! we side
640                  IF( mpiRnei(nn_hls,jpwe) > -1 ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpwe,ir) = .TRUE.   ! rcv from we nei
641                  ELSE                                   ;   CALL ctl_stop( 'bdyini recv olr we-side ' )
642                  ENDIF
643               ENDIF
644               IF( ii  > Nie0                  .AND. ij >= Njs0 .AND. ij <= Nje0 ) THEN   ! ea side
645                  IF( mpiRnei(nn_hls,jpea) > -1 ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpea,ir) = .TRUE.   ! rcv from ea nei
646                  ELSE                                   ;   CALL ctl_stop( 'bdyini recv olr ea-side ' )
647                  ENDIF
648               ENDIF
649               IF( ii >= Nis0 .AND. ii <= Nie0 .AND. ij  < Njs0                  ) THEN   ! so side
650                  IF( mpiRnei(nn_hls,jpso) > -1 ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.   ! rcv from so nei
651                  ELSE                                   ;   CALL ctl_stop( 'bdyini recv olr so-side ' )
652                  ENDIF
653               ENDIF
654               IF( ii >= Nis0 .AND. ii <= Nie0 .AND. ij  > Nje0                  ) THEN   ! no side
655                  IF( mpiRnei(nn_hls,jpno) > -1 ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.   ! rcv from no nei
656                  ELSE                                   ;   CALL ctl_stop( 'bdyini recv olr no-side ' )
657                  ENDIF
658               ENDIF
659               IF( ii  < Nis0                  .AND. ij  < Njs0                  ) THEN   ! sw corner
660                  IF(     mpiRnei(nn_hls,jpsw) > -1                    ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpsw,ir) = .TRUE.
661                  ELSEIF( mpiRnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
662                  ELSE                                                          ;   CALL ctl_stop( 'bdyini recv olr sw-corner' )
663                  ENDIF
664               ENDIF
665               IF( ii  > Nie0                  .AND. ij  < Njs0                  ) THEN   ! se corner
666                  IF(     mpiRnei(nn_hls,jpse) > -1                    ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpse,ir) = .TRUE.
667                  ELSEIF( mpiRnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
668                  ELSE                                                          ;   CALL ctl_stop( 'bdyini recv olr se-corner' )
669                  ENDIF
670               ENDIF
671               IF( ii  < Nis0                  .AND. ij  > Nje0                  ) THEN   ! nw corner
672                  IF(     mpiRnei(nn_hls,jpnw) > -1                    ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpnw,ir) = .TRUE.
673                  ELSEIF( mpiRnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
674                  ELSE                                                          ;   CALL ctl_stop( 'bdyini recv olr nw-corner' )
675                  ENDIF
676               ENDIF
677               IF( ii  > Nie0                  .AND. ij  > Nje0                  ) THEN   ! ne corner
678                  IF(     mpiRnei(nn_hls,jpne) > -1                    ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpne,ir) = .TRUE.
679                  ELSEIF( mpiRnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 ) THEN   ;   lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
680                  ELSE                                                          ;   CALL ctl_stop( 'bdyini recv olr ne-corner' )
681                  ENDIF
682               ENDIF
683               !
684            END DO
685         END DO   !   igrd
686         
687         ! Comment out for debug
688!!$         DO ir = 0,1
689!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing,   &
690!!$               &                              lsend = lsend_bdyolr(ib_bdy,1,:,ir), lrecv = lrecv_bdyolr(ib_bdy,1,:,ir) )
691!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug olr T', ir ; CALL FLUSH(numout)
692!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing,   &
693!!$               &                              lsend = lsend_bdyolr(ib_bdy,2,:,ir), lrecv = lrecv_bdyolr(ib_bdy,2,:,ir) )
694!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug olr U', ir ; CALL FLUSH(numout)
695!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing,   &
696!!$               &                              lsend = lsend_bdyolr(ib_bdy,3,:,ir), lrecv = lrecv_bdyolr(ib_bdy,3,:,ir) )   
697!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug olr V', ir ; CALL FLUSH(numout)
698!!$         END DO
699         
700         ! Compute rim weights for FRS scheme
701         ! ----------------------------------
702         DO igrd = 1, jpbgrd
703            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
704               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same weights
705               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 )      ! tanh formulation
706               !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic
707               !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy))       ! linear
708            END DO
709         END DO
710
711         ! Compute damping coefficients
712         ! ----------------------------
713         DO igrd = 1, jpbgrd
714            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
715               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same damping coefficients
716               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 
717                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
718               idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) & 
719                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
720            END DO
721         END DO
722
723      END DO ! ib_bdy
724
725      ! ------------------------------------------------------
726      ! Initialise masks and find normal/tangential directions
727      ! ------------------------------------------------------
728
729      ! ------------------------------------------
730      ! handle rim0, do as if rim 1 was free ocean
731      ! ------------------------------------------
732
733      ztmask(:,:) = tmask(:,:,1)   ;   zumask(:,:) = umask(:,:,1)   ;   zvmask(:,:) = vmask(:,:,1)
734      ! For the flagu/flagv calculation below we require a version of fmask without
735      ! the land boundary condition (shlat) included:
736      DO_2D( 0, 0, 0, 0 )
737         zfmask(ji,jj) =  ztmask(ji,jj  ) * ztmask(ji+1,jj  )   &
738            &           * ztmask(ji,jj+1) * ztmask(ji+1,jj+1)
739      END_2D
740      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp )
741
742      ! Read global 2D mask at T-points: bdytmask
743      ! -----------------------------------------
744      ! bdytmask = 1  on the computational domain but not on open boundaries
745      !          = 0  elsewhere   
746
747      bdytmask(:,:) = ssmask(:,:)
748
749      ! Derive mask on U and V grid from mask on T grid
750      DO_2D( 0, 0, 0, 0 )
751            bdyumask(ji,jj) = bdytmask(ji,jj) * bdytmask(ji+1,jj  )
752            bdyvmask(ji,jj) = bdytmask(ji,jj) * bdytmask(ji  ,jj+1) 
753      END_2D
754      CALL lbc_lnk( 'bdyini', bdyumask, 'U', 1.0_wp , bdyvmask, 'V', 1.0_wp )   ! Lateral boundary cond.
755
756      ! bdy masks are now set to zero on rim 0 points:
757      DO ib_bdy = 1, nb_bdy
758         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0
759            bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
760         END DO
761         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2)   ! extent of rim 0
762            bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
763         END DO
764         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3)   ! extent of rim 0
765            bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
766         END DO
767      END DO
768
769      CALL bdy_rim_treat( zumask, zvmask, zfmask, .true. )   ! compute flagu, flagv, ntreat on rim 0
770
771      ! ------------------------------------
772      ! handle rim1, do as if rim 0 was land
773      ! ------------------------------------
774     
775      ! z[tuv]mask are now set to zero on rim 0 points:
776      DO ib_bdy = 1, nb_bdy
777         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0
778            ztmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
779         END DO
780         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2)   ! extent of rim 0
781            zumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
782         END DO
783         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3)   ! extent of rim 0
784            zvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
785         END DO
786      END DO
787
788      ! Recompute zfmask
789      DO_2D( 0, 0, 0, 0 )
790         zfmask(ji,jj) =  ztmask(ji,jj  ) * ztmask(ji+1,jj  )   &
791            &           * ztmask(ji,jj+1) * ztmask(ji+1,jj+1)
792      END_2D
793      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp )
794
795      ! bdy masks are now set to zero on rim1 points:
796      DO ib_bdy = 1, nb_bdy
797         DO ib = idx_bdy(ib_bdy)%nblenrim0(1) + 1,  idx_bdy(ib_bdy)%nblenrim(1)   ! extent of rim 1
798            bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
799         END DO
800         DO ib = idx_bdy(ib_bdy)%nblenrim0(2) + 1,  idx_bdy(ib_bdy)%nblenrim(2)   ! extent of rim 1
801            bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
802         END DO
803         DO ib = idx_bdy(ib_bdy)%nblenrim0(3) + 1,  idx_bdy(ib_bdy)%nblenrim(3)   ! extent of rim 1
804            bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
805         END DO
806      END DO
807
808      CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. )   ! compute flagu, flagv, ntreat on rim 1
809      !
810      ! Check which boundaries might need communication
811      ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,8,0:1), lrecv_bdyint(nb_bdy,jpbgrd,8,0:1) )
812      lsend_bdyint(:,:,:,:) = .false.
813      lrecv_bdyint(:,:,:,:) = .false. 
814      ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,8,0:1), lrecv_bdyext(nb_bdy,jpbgrd,8,0:1) )
815      lsend_bdyext(:,:,:,:) = .false.
816      lrecv_bdyext(:,:,:,:) = .false.
817      !
818      DO ib_bdy = 1, nb_bdy
819         DO igrd = 1, jpbgrd
820            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
821               IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE
822               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
823               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
824               ir = idx_bdy(ib_bdy)%nbr(ib,igrd)
825               flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd))
826               flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd))
827               iibe = ii - flagu   ! neighbouring point towards the exterior of the computational domain
828               ijbe = ij - flagv
829               iibi = ii + flagu   ! neighbouring point towards the interior of the computational domain
830               ijbi = ij + flagv
831               CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 )   ! free ocean neighbours
832               !
833               !  take care of the 4 sides
834               !
835               DO icnt = 1, 4
836                  SELECT CASE( icnt )
837                     !                                           ... _____
838                  CASE( 1 )   ! x: rim on rcvwe/sndea-side         o|  :
839                     !          o: potential neighbour(s)          o|x :
840                     !             outside of the MPI domain     ..o|__:__
841                     cRdir    = 'we'             ;   cSdir    = 'ea'
842                     iRnei    = jpwe             ;   iSnei    = jpea
843                     iiRst    = 1                ;   ijRst    = 2               ! Rcv we-side starting point, excluding sw-corner
844                     iiRnd    = 1                ;   ijRnd    = jpj-1           ! Rcv we-side   ending point, excluding nw-corner
845                     iiSst    = jpi-2*nn_hls+1   ;   ijSst    = 2               ! Snd ea-side starting point, excluding se-corner
846                     iiSnd    = jpi-2*nn_hls+1   ;   ijSnd    = jpj-1           ! Snd ea-side   ending point, excluding ne-corner
847                     iioutdir = -1               ;   ijoutdir = -999            ! outside MPI domain: westward
848                     !                                           ______....
849                  CASE( 2 )   ! x: rim on rcvea/sndwe-side          :  |o
850                     !          o: potential neighbour(s)           : x|o
851                     !             outside of the MPI domain     ___:__|o..
852                     cRdir    = 'ea'             ;   cSdir    = 'we'
853                     iRnei    = jpea             ;   iSnei    = jpwe
854                     iiRst    = jpi              ;   ijRst    = 2                ! Rcv ea-side starting point, excluding se-corner
855                     iiRnd    = jpi              ;   ijRnd    = jpj-1            ! Rcv ea-side   ending point, excluding ne-corner
856                     iiSst    = 2*nn_hls         ;   ijSst    = 2                ! Snd we-side starting point, excluding sw-corner
857                     iiSnd    = 2*nn_hls         ;   ijSnd    = jpj-1            ! Snd we-side   ending point, excluding nw-corner
858                     iioutdir = 1                ;   ijoutdir = -999             ! outside MPI domain: eastward
859                     !
860                  CASE( 3 )   ! x: rim on rcvso/sndno-side       |       |
861                     !          o: potential neighbour(s)        |¨¨¨¨¨¨¨|
862                     !             outside of the MPI domain     |___x___|
863                     !                                           : o o o :
864                     !                                           :       :
865                     cRdir    = 'so'             ;   cSdir    = 'no'
866                     iRnei    = jpso             ;   iSnei    = jpno
867                     iiRst    = 2                ;   ijRst    = 1                ! Rcv so-side starting point, excluding sw-corner
868                     iiRnd    = jpi-1            ;   ijRnd    = 1                ! Rcv so-side   ending point, excluding se-corner
869                     iiSst    = 2                ;   ijSst    = jpj-2*nn_hls+1   ! Snd no-side starting point, excluding nw-corner
870                     iiSnd    = jpi-1            ;   ijSnd    = jpj-2*nn_hls+1   ! Snd no-side   ending point, excluding ne-corner
871                     iioutdir = -999             ;   ijoutdir = -1               ! outside MPI domain: southward
872                     !                                           :       :
873                  CASE( 4 )   ! x: rim on rcvno/sndso-side       :_o_o_o_:
874                     !          o: potential neighbour(s)        |   x   |
875                     !             outside of the MPI domain     |       |
876                     !                                           |¨¨¨¨¨¨¨|
877                     cRdir    = 'no'             ;   cSdir    = 'so'
878                     iRnei    = jpno             ;   iSnei    = jpso
879                     iiRst    = 2                ;   ijRst    = jpj              ! Rcv no-side starting point, excluding nw-corner
880                     iiRnd    = jpi-1            ;   ijRnd    = jpj              ! Rcv no-side   ending point, excluding ne-corner
881                     iiSst    = 2                ;   ijSst    =     2*nn_hls     ! Snd so-side starting point, excluding sw-corner
882                     iiSnd    = jpi-1            ;   ijSnd    =     2*nn_hls     ! Snd so-side   ending point, excluding se-corner
883                     iioutdir = -999             ;   ijoutdir = 1                ! outside MPI domain: northward
884                  END SELECT
885                  !
886                  IF( ii >= iiRst .AND. ii <= iiRnd .AND. ij >= ijRst .AND. ij <= ijRnd ) THEN   ! rim point in recv side
887                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the MPI domain?
888                     ! take care of neighbourg(s) in the interior of the computational domain
889                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of the MPI domain
890                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN     ! -> I cannot compute it -> recv it
891                        IF( mpiRnei(nn_hls,iRnei) > -1 ) THEN   ;   lrecv_bdyint(ib_bdy,igrd,iRnei,ir) = .TRUE.
892                        ELSE                                    ;   CALL ctl_stop( 'bdyini recv int '//cRdir//'-side ' )
893                        ENDIF
894                     ENDIF
895                     ! take care of neighbourg in the exterior of the computational domain
896                     IF(  iibe==iiout .OR. ijbe==ijout ) THEN   ! Neib outside of the MPI domain -> I cannot compute it -> recv it
897                        IF( mpiRnei(nn_hls,iRnei) > -1 ) THEN   ;   lrecv_bdyext(ib_bdy,igrd,iRnei,ir) = .TRUE.
898                        ELSE                                    ;   CALL ctl_stop( 'bdyini recv ext '//cRdir//'-side ' )
899                        ENDIF
900                     ENDIF
901                  ENDIF
902                 
903                  IF( ii >= iiSst .AND. ii <= iiSnd .AND. ij >= ijSst .AND. ij <= ijSnd ) THEN   ! rim point in send side
904                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the nei MPI domain?
905                     ! take care of neighbourg(s) in the interior of the computational domain
906                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of nei MPI domain
907                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN     ! -> nei cannot compute it
908                        IF( mpiSnei(nn_hls,iSnei) > -1 ) THEN   ;   lsend_bdyint(ib_bdy,igrd,iSnei,ir) = .TRUE.   ! -> send to nei
909                        ELSE                                    ;   CALL ctl_stop( 'bdyini send int '//cSdir//'-side ' )
910                        ENDIF
911                     ENDIF
912                     ! take care of neighbourg in the exterior of the computational domain
913                     IF( iibe == iiout .OR. ijbe == ijout ) THEN   ! Neib outside of the nei MPI domain -> nei cannot compute it
914                        IF( mpiSnei(nn_hls,iSnei) > -1 ) THEN   ;   lsend_bdyext(ib_bdy,igrd,iSnei,ir) = .TRUE.   ! -> send to nei
915                        ELSE                                    ;   CALL ctl_stop( 'bdyini send ext '//cSdir//'-side ' )
916                        ENDIF
917                     ENDIF
918                  END IF
919
920               END DO   ! 4 sides
921               !
922               ! specific treatment for the corners
923               !
924               DO icnt = 1, 4
925                  SELECT CASE( icnt )
926                     !                                       ...|....
927                  CASE( 1 )   ! x: rim on sw-corner            o|   :
928                     !          o: potential neighbour(s)      o|x__:__
929                     !             outside of the MPI domain   o o o:
930                     !                                              :
931                     cRdir    = 'sw'
932                     iRdiag   = jpsw             ;   iRsono   = jpso             ! Recv: for sw or so
933                     iSdiag   = jpne             ;   iSsono   = jpno             ! Send: to ne or no
934                     iiRcorn  = 1                ;   ijRcorn  = 1                ! receiving sw-corner
935                     iiSdiag  = jpi-2*nn_hls+1   ;   ijSdiag  = jpj-2*nn_hls+1   ! send to sw-corner of ne neighbourg
936                     iiSsono  = 1                ;   ijSsono  = jpj-2*nn_hls+1   ! send to sw-corner of no neighbourg
937                     iioutdir = -1               ;   ijoutdir = -1               ! outside MPI domain: westward or southward
938                     !                                          ....|...
939                  CASE( 2 )   ! x: rim on se-corner             :   |o
940                     !          o: potential neighbour(s)     __:__x|o
941                     !             outside of the MPI domain    :o o o
942                     !                                          :   
943                     cRdir    = 'se'
944                     iRdiag   = jpse             ;   iRsono   = jpso             ! Recv: for se or so
945                     iSdiag   = jpnw             ;   iSsono   = jpno             ! Send: to nw or no
946                     iiRcorn  = jpi              ;   ijRcorn  = 1                ! receiving se-corner
947                     iiSdiag  =     2*nn_hls     ;   ijSdiag  = jpj-2*nn_hls+1   ! send to se-corner of nw neighbourg
948                     iiSsono  = jpi              ;   ijSsono  = jpj-2*nn_hls+1   ! send to se-corner of no neighbourg
949                     iioutdir = 1                ;   ijoutdir = -1               ! outside MPI domain: eastward or southward
950                     !                                              :       
951                     !                                         o o_o:___
952                  CASE( 3 )   ! x: rim on nw-corner            o|x  :
953                     !          o: potential neighbour(s)    ..o|...:
954                     !             outside of the MPI domain    |
955                     cRdir    = 'nw'
956                     iRdiag   = jpnw             ;   iRsono   = jpno             ! Recv: for nw or no
957                     iSdiag   = jpse             ;   iSsono   = jpso             ! Send: to se or so
958                     iiRcorn  = 1                ;   ijRcorn  = jpj              ! receiving nw-corner
959                     iiSdiag  = jpi-2*nn_hls+1   ;   ijSdiag  =     2*nn_hls     ! send to nw-corner of se neighbourg
960                     iiSsono  = 1                ;   ijSsono  =     2*nn_hls     ! send to nw-corner of so neighbourg
961                     iioutdir = -1               ;   ijoutdir =  1               ! outside MPI domain: westward or northward
962                     !                                          :       
963                     !                                       ___:o_o o
964                  CASE( 4 )   ! x: rim on ne-corner             :  x|o
965                     !          o: potential neighbour(s)       :...|o...
966                     !             outside of the MPI domain        |
967                     cRdir    = 'ne'
968                     iRdiag   = jpne             ;   iRsono   = jpno             ! Recv: for ne or no
969                     iSdiag   = jpsw             ;   iSsono   = jpso             ! Send: to sw or so
970                     iiRcorn  = jpi              ;   ijRcorn  = jpj              ! receiving ne-corner
971                     iiSdiag  =     2*nn_hls     ;   ijSdiag  =     2*nn_hls     ! send to ne-corner of sw neighbourg
972                     iiSsono  = jpi              ;   ijSsono  =     2*nn_hls     ! send to ne-corner of so neighbourg
973                     iioutdir = 1                ;   ijoutdir = 1                ! outside MPI domain: eastward or southward
974                  END SELECT
975                  !
976                  ! Check if we need to receive data for this rim point
977                  IF( ii == iiRcorn .AND. ij == ijRcorn ) THEN   ! the rim point is located on the corner for the MPI domain
978                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the MPI domain?
979                     ! take care of neighbourg(s) in the interior of the computational domain
980                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of the MPI domain
981                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN     ! -> I cannot compute it -> recv it
982                        IF(     mpiRnei(nn_hls,iRdiag) > -1 ) THEN
983                           lrecv_bdyint(ib_bdy,igrd,iRdiag,ir) = .TRUE.   ! Receive directly from diagonal neighbourg
984                        ELSEIF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 ) THEN
985                           lrecv_bdyint(ib_bdy,igrd,iRsono,ir) = .TRUE.   ! Receive through the South/North neighbourg
986                        ELSE
987                           CALL ctl_stop( 'bdyini recv int '//cRdir//'-corner ' )
988                        ENDIF
989                     ENDIF
990                     ! take care of neighbourg in the exterior of the computational domain
991                     IF(  iibe==iiout .OR. ijbe==ijout ) THEN   ! Neib outside of the MPI domain -> I cannot compute it -> recv it
992                        IF(     mpiRnei(nn_hls,iRdiag) > -1 ) THEN
993                           lrecv_bdyext(ib_bdy,igrd,iRdiag,ir) = .TRUE.   ! Receive directly from diagonal neighbourg
994                        ELSEIF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 ) THEN
995                           lrecv_bdyext(ib_bdy,igrd,iRsono,ir) = .TRUE.   ! Receive through the South/North neighbourg
996                        ELSE
997                           CALL ctl_stop( 'bdyini recv ext '//cRdir//'-corner ' )
998                        ENDIF
999                     ENDIF
1000                  ENDIF
1001                  !
1002                  ! Check if this rim point corresponds to the corner of one neighbourg. if yes, do we need to send data?
1003                  ! Direct send to diag: Is this rim point the corner point of a diag neighbour with which we communicate?
1004                  IF( ii == iiSdiag .AND. ij == ijSdiag .AND. mpiSnei(nn_hls,iSdiag) > -1 ) THEN
1005                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the nei MPI domain?
1006                     ! take care of neighbourg(s) in the interior of the computational domain
1007                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of diag nei MPI
1008                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout )      &   ! domain -> nei cannot compute it
1009                        &    lsend_bdyint(ib_bdy,igrd,iSdiag,ir) = .TRUE.                        ! send rim point data to diag nei
1010                     ! take care of neighbourg in the exterior of the computational domain
1011                     IF(  iibe==iiout .OR. ijbe==ijout )   &                                 
1012                        &    lsend_bdyext(ib_bdy,igrd,iSdiag,ir) = .TRUE.
1013                  ENDIF
1014                  ! Indirect send to diag (through so/no): rim point is the corner point of a so/no nei with which we communicate
1015                  IF( ii == iiSsono .AND. ij == ijSsono .AND. mpiSnei(nn_hls,iSsono) > -1 .AND. nn_comm == 1 ) THEN
1016                     iiout = ii+iioutdir ; ijout = ij+ijoutdir        ! in which direction do we go outside of the nei MPI domain?
1017                     ! take care of neighbourg(s) in the interior of the computational domain
1018                     IF(  iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR.   &   ! Neib outside of so/no nei MPI
1019                        & ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout )      &   ! domain -> nei cannot compute it
1020                        &    lsend_bdyint(ib_bdy,igrd,iSsono,ir) = .TRUE.                        ! send rim point data to so/no nei
1021                     ! take care of neighbourg in the exterior of the computational domain
1022                     IF(  iibe==iiout .OR. ijbe==ijout)   &
1023                        &    lsend_bdyext(ib_bdy,igrd,iSsono,ir) = .TRUE.
1024                  ENDIF
1025                  !
1026               END DO   ! 4 corners
1027            END DO   ! ib
1028         END DO   ! igrd
1029
1030         ! Comment out for debug
1031!!$         DO ir = 0,1
1032!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing,   &
1033!!$               &                              lsend = lsend_bdyint(ib_bdy,1,:,ir), lrecv = lrecv_bdyint(ib_bdy,1,:,ir) )
1034!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug int T', ir ; CALL FLUSH(numout)
1035!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing,   &
1036!!$               &                              lsend = lsend_bdyint(ib_bdy,2,:,ir), lrecv = lrecv_bdyint(ib_bdy,2,:,ir) )
1037!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug int U', ir ; CALL FLUSH(numout)
1038!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing,   &
1039!!$               &                              lsend = lsend_bdyint(ib_bdy,3,:,ir), lrecv = lrecv_bdyint(ib_bdy,3,:,ir) )   
1040!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug int V', ir ; CALL FLUSH(numout)
1041!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing,   &
1042!!$               &                              lsend = lsend_bdyext(ib_bdy,1,:,ir), lrecv = lrecv_bdyext(ib_bdy,1,:,ir) )
1043!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug ext T', ir ; CALL FLUSH(numout)
1044!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing,   &
1045!!$               &                              lsend = lsend_bdyext(ib_bdy,2,:,ir), lrecv = lrecv_bdyext(ib_bdy,2,:,ir) )
1046!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug ext U', ir ; CALL FLUSH(numout)
1047!!$            zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing,   &
1048!!$               &                              lsend = lsend_bdyext(ib_bdy,3,:,ir), lrecv = lrecv_bdyext(ib_bdy,3,:,ir) )   
1049!!$            IF(lwp) WRITE(numout,*) ' seb bdy debug ext V', ir ; CALL FLUSH(numout)
1050!!$         END DO
1051         
1052      END DO   ! ib_bdy
1053
1054      DO ib_bdy = 1,nb_bdy
1055         IF(  cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &
1056            & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &
1057            & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN
1058            DO igrd = 1, jpbgrd
1059               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1060                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1061                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1062                  IF(  mig0(ii) > 2 .AND. mig0(ii) < Ni0glo-2 .AND. mjg0(ij) > 2 .AND. mjg0(ij) < Nj0glo-2  ) THEN
1063                     WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain'
1064                     CALL ctl_stop( ctmp1 )
1065                  END IF
1066               END DO
1067            END DO
1068         END IF
1069      END DO
1070      !
1071      DEALLOCATE( nbidta, nbjdta, nbrdta )
1072      !
1073   END SUBROUTINE bdy_def
1074
1075
1076   SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 )
1077      !!----------------------------------------------------------------------
1078      !!                 ***  ROUTINE bdy_rim_treat  ***
1079      !!
1080      !! ** Purpose :   Initialize structures ( flagu, flagv, ntreat ) indicating how rim points
1081      !!                  are to be handled in the boundary condition treatment
1082      !!
1083      !! ** Method  :   - to handle rim 0 zmasks must indicate ocean points      (set at one on rim 0 and rim 1 and interior)
1084      !!                            and bdymasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
1085      !!                    (as if rim 1 was free ocean)
1086      !!                - to handle rim 1 zmasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
1087      !!                            and bdymasks must indicate free ocean points (set at one on interior)
1088      !!                    (as if rim 0 was land)
1089      !!                - we can then check in which direction the interior of the computational domain is with the difference
1090      !!                         mask array values on both sides to compute flagu and flagv
1091      !!                - and look at the ocean neighbours to compute ntreat
1092      !!----------------------------------------------------------------------
1093      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pumask, pvmask   ! temporary u/v mask array
1094      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pfmask           ! temporary fmask excluding coastal boundary condition (shlat)
1095      LOGICAL                             , INTENT (in   ) :: lrim0            ! .true. -> rim 0   .false. -> rim 1
1096      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, icount       ! dummy loop indices
1097      INTEGER  ::   i_offset, j_offset, inn                ! local integer
1098      INTEGER  ::   ibeg, iend                             ! local integer
1099      LOGICAL  ::   llnon, llson, llean, llwen             ! local logicals indicating the presence of a ocean neighbour
1100      REAL(wp), POINTER, DIMENSION(:,:)       ::   zmask   ! pointer to 2D mask fields
1101      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
1102      CHARACTER(LEN=1), DIMENSION(jpbgrd)     ::   cgrid
1103      REAL(wp)        , DIMENSION(jpi,jpj)    ::   ztmp
1104      !!----------------------------------------------------------------------
1105
1106      cgrid = (/'t','u','v'/)
1107
1108      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
1109
1110         DO igrd = 1, jpbgrd
1111           
1112            IF( lrim0 ) THEN   ! extent of rim 0
1113               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
1114            ELSE               ! extent of rim 1
1115               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
1116            END IF
1117
1118            ! Calculate relationship of U direction to the local orientation of the boundary
1119            ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
1120            ! flagu =  0 : u is tangential
1121            ! flagu =  1 : u is normal to the boundary and is direction is inward
1122            SELECT CASE( igrd )
1123               CASE( 1 )   ;   zmask => pumask     ;   i_offset = 0   ! U(i-1)   T(i)   U(i  )
1124               CASE( 2 )   ;   zmask => bdytmask   ;   i_offset = 1   ! T(i  )   U(i)   T(i+1)
1125               CASE( 3 )   ;   zmask => pfmask     ;   i_offset = 0   ! F(i-1)   V(i)   F(i  )
1126            END SELECT
1127            icount = 0
1128            ztmp(:,:) = -999._wp
1129            DO ib = ibeg, iend 
1130               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1131               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1132               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
1133               zwfl = zmask(ii+i_offset-1,ij)
1134               zefl = zmask(ii+i_offset  ,ij)
1135               ! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
1136               IF( i_offset == 1 .and. zefl + zwfl == 2._wp ) THEN
1137                  icount = icount + 1
1138                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
1139               ELSE
1140                  ztmp(ii,ij) = -zwfl + zefl
1141               ENDIF
1142            END DO
1143            IF( icount /= 0 ) THEN
1144               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
1145                  ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1146               CALL ctl_stop( ctmp1 )
1147            ENDIF
1148            SELECT CASE( igrd )
1149               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1150               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1151               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1152            END SELECT
1153            DO ib = ibeg, iend
1154               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1155               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1156               idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij)
1157            END DO
1158           
1159            ! Calculate relationship of V direction to the local orientation of the boundary
1160            ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
1161            ! flagv =  0 : v is tangential
1162            ! flagv =  1 : v is normal to the boundary and is direction is inward
1163            SELECT CASE( igrd )
1164               CASE( 1 )   ;   zmask => pvmask     ;   j_offset = 0
1165               CASE( 2 )   ;   zmask => pfmask     ;   j_offset = 0
1166               CASE( 3 )   ;   zmask => bdytmask   ;   j_offset = 1
1167            END SELECT
1168            icount = 0
1169            ztmp(:,:) = -999._wp
1170            DO ib = ibeg, iend
1171               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1172               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1173               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
1174               zsfl = zmask(ii,ij+j_offset-1)
1175               znfl = zmask(ii,ij+j_offset  )
1176               ! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
1177               IF( j_offset == 1 .and. znfl + zsfl == 2._wp ) THEN
1178                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
1179                  icount = icount + 1
1180               ELSE
1181                  ztmp(ii,ij) = -zsfl + znfl
1182               END IF
1183            END DO
1184            IF( icount /= 0 ) THEN
1185               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
1186                  ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1187               CALL ctl_stop( ctmp1 )
1188            ENDIF
1189            SELECT CASE( igrd )
1190               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1191               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1192               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1193            END SELECT
1194            DO ib = ibeg, iend
1195               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1196               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1197               idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij)
1198            END DO
1199     
1200            ! Calculate ntreat
1201            SELECT CASE( igrd )
1202               CASE( 1 )   ;   zmask => bdytmask 
1203               CASE( 2 )   ;   zmask => bdyumask 
1204               CASE( 3 )   ;   zmask => bdyvmask 
1205            END SELECT
1206            ztmp(:,:) = -999._wp
1207            DO ib = ibeg, iend
1208               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1209               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1210               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
1211               llnon = zmask(ii  ,ij+1) == 1._wp 
1212               llson = zmask(ii  ,ij-1) == 1._wp 
1213               llean = zmask(ii+1,ij  ) == 1._wp 
1214               llwen = zmask(ii-1,ij  ) == 1._wp 
1215               inn  = COUNT( (/ llnon, llson, llean, llwen /) )
1216               IF( inn == 0 ) THEN   ! no neighbours -> interior of a corner  or  cluster of rim points
1217                  !               !              !     _____     !     _____    !    __     __
1218                  !  1 |   o      !  2  o   |    !  3 | x        !  4     x |   !      |   |   -> error
1219                  !    |_x_ _     !    _ _x_|    !    |   o      !      o   |   !      |x_x|
1220                  IF(     zmask(ii+1,ij+1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 1._wp
1221                  ELSEIF( zmask(ii-1,ij+1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 2._wp
1222                  ELSEIF( zmask(ii+1,ij-1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 3._wp
1223                  ELSEIF( zmask(ii-1,ij-1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 4._wp
1224                  ELSE                                       ;   ztmp(ii,ij) = -1._wp
1225                     WRITE(ctmp1,*) 'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
1226                       ' on boundary set ', ib_bdy, ' has no free ocean neighbour'
1227                     IF( lrim0 ) THEN
1228                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 0 points.'
1229                     ELSE
1230                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.'
1231                     END IF
1232                     CALL ctl_warn( ctmp1, ctmp2 )
1233                  END IF
1234               END IF
1235               IF( inn == 1 ) THEN   ! middle of linear bdy  or incomplete corner  ! ___ o
1236                  !    |         !         |   !      o     !    ______            !    |x___
1237                  ! 5  | x o     ! 6   o x |   ! 7  __x__   ! 8    x
1238                  !    |         !         |   !            !      o
1239                  IF( llean )   ztmp(ii,ij) = 5._wp
1240                  IF( llwen )   ztmp(ii,ij) = 6._wp
1241                  IF( llnon )   ztmp(ii,ij) = 7._wp
1242                  IF( llson )   ztmp(ii,ij) = 8._wp
1243               END IF
1244               IF( inn == 2 ) THEN   ! exterior of a corner
1245                  !        o      !        o      !    _____|       !       |_____ 
1246                  !  9 ____x o    ! 10   o x___   ! 11     x o      ! 12   o x     
1247                  !         |     !       |       !        o        !        o
1248                  IF( llnon .AND. llean )   ztmp(ii,ij) =  9._wp
1249                  IF( llnon .AND. llwen )   ztmp(ii,ij) = 10._wp
1250                  IF( llson .AND. llean )   ztmp(ii,ij) = 11._wp
1251                  IF( llson .AND. llwen )   ztmp(ii,ij) = 12._wp
1252               END IF
1253               IF( inn == 3 ) THEN   ! 3 neighbours     __   __
1254                  !    |_  o      !        o  _|  !       |_|     !       o         
1255                  ! 13  _| x o    ! 14   o x |_   ! 15   o x o    ! 16  o x o       
1256                  !    |   o      !        o   |  !        o      !    __|¨|__   
1257                  IF( llnon .AND. llean .AND. llson )   ztmp(ii,ij) = 13._wp
1258                  IF( llnon .AND. llwen .AND. llson )   ztmp(ii,ij) = 14._wp
1259                  IF( llwen .AND. llson .AND. llean )   ztmp(ii,ij) = 15._wp
1260                  IF( llwen .AND. llnon .AND. llean )   ztmp(ii,ij) = 16._wp
1261               END IF
1262               IF( inn == 4 ) THEN
1263                  WRITE(ctmp1,*)  'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
1264                       ' on boundary set ', ib_bdy, ' have 4 neighbours'
1265                  CALL ctl_stop( ctmp1 )
1266               END IF
1267            END DO
1268            SELECT CASE( igrd )
1269               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1270               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1271               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1272            END SELECT
1273            DO ib = ibeg, iend
1274               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1275               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1276               idx_bdy(ib_bdy)%ntreat(ib,igrd) = NINT(ztmp(ii,ij))
1277            END DO
1278            !
1279         END DO   ! jpbgrd
1280         !
1281      END DO   ! ib_bdy
1282
1283    END SUBROUTINE bdy_rim_treat
1284
1285   
1286    SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )
1287      !!----------------------------------------------------------------------
1288      !!                 ***  ROUTINE find_neib  ***
1289      !!
1290      !! ** Purpose :   get ii1, ij1, ii2, ij2, ii3, ij3, the indices of
1291      !!               the free ocean neighbours of (ii,ij) for bdy treatment
1292      !!
1293      !! ** Method  :  use itreat input to select a case
1294      !!               N.B. ntreat is defined for all bdy points in routine bdy_rim_treat
1295      !!
1296      !!----------------------------------------------------------------------
1297      INTEGER, INTENT(in   )      ::   ii, ij, itreat
1298      INTEGER, INTENT(  out)      ::   ii1, ij1, ii2, ij2, ii3, ij3
1299      !!----------------------------------------------------------------------
1300      SELECT CASE( itreat )   ! points that will be used by bdy routines, -1 will be discarded
1301         !               !               !     _____     !     _____     
1302         !  1 |   o      !  2  o   |     !  3 | x        !  4     x |   
1303         !    |_x_ _     !    _ _x_|     !    |   o      !      o   |
1304      CASE( 1 )    ;   ii1 = ii+1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1305      CASE( 2 )    ;   ii1 = ii-1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1306      CASE( 3 )    ;   ii1 = ii+1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1307      CASE( 4 )    ;   ii1 = ii-1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1308         !    |          !         |     !      o        !    ______                   ! or incomplete corner
1309         ! 5  | x o      ! 6   o x |     ! 7  __x__      ! 8    x                      !  7  ____ o
1310         !    |          !         |     !               !      o                      !         |x___
1311      CASE( 5 )    ;   ii1 = ii+1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1312      CASE( 6 )    ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1313      CASE( 7 )    ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1314      CASE( 8 )    ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1315         !        o      !        o      !    _____|     !       |_____ 
1316         !  9 ____x o    ! 10   o x___   ! 11     x o    ! 12   o x     
1317         !         |     !       |       !        o      !        o     
1318      CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
1319      CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1320      CASE( 11 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1321      CASE( 12 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1322         !    |_  o      !        o  _|  !     ¨¨|_|¨¨   !       o         
1323         ! 13  _| x o    !  14  o x |_   !  15  o x o    ! 16  o x o       
1324         !    |   o      !        o   |  !        o      !    __|¨|__
1325      CASE( 13 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1   
1326      CASE( 14 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1 
1327      CASE( 15 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij   
1328      CASE( 16 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij
1329      END SELECT
1330   END SUBROUTINE find_neib
1331   
1332
1333   SUBROUTINE bdy_read_seg( kb_bdy, knblendta ) 
1334      !!----------------------------------------------------------------------
1335      !!                 ***  ROUTINE bdy_read_seg  ***
1336      !!
1337      !! ** Purpose :  build bdy coordinates with segments defined in namelist
1338      !!
1339      !! ** Method  :  read namelist nambdy_index blocks
1340      !!
1341      !!----------------------------------------------------------------------
1342      INTEGER                   , INTENT (in   ) ::   kb_bdy           ! bdy number
1343      INTEGER, DIMENSION(jpbgrd), INTENT (  out) ::   knblendta        ! length of index arrays
1344      !!
1345      INTEGER          ::   ios                 ! Local integer output status for namelist read
1346      INTEGER          ::   nbdyind, nbdybeg, nbdyend
1347      INTEGER          ::   nbdy_count, nbdy_rdstart, nbdy_loc
1348      CHARACTER(LEN=1) ::   ctypebdy   !     -        -
1349      CHARACTER(LEN=50)::   cerrmsg    !     -        -
1350      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
1351      !!----------------------------------------------------------------------
1352      ! Need to support possibility of reading more than one nambdy_index from
1353      ! the namelist_cfg internal file.
1354      ! Do this by finding the kb_bdy'th occurence of nambdy_index in the
1355      ! character buffer as the starting point.
1356      nbdy_rdstart = 1
1357      DO nbdy_count = 1, kb_bdy
1358       nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_index' )
1359       IF( nbdy_loc .GT. 0 ) THEN
1360          nbdy_rdstart = nbdy_rdstart + nbdy_loc
1361       ELSE
1362          WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',kb_bdy,' of nambdy_index not found'
1363          ios = -1
1364          CALL ctl_nam ( ios , cerrmsg )
1365       ENDIF
1366      END DO
1367      nbdy_rdstart = MAX( 1, nbdy_rdstart - 2 )
1368      READ  ( numnam_cfg( nbdy_rdstart: ), nambdy_index, IOSTAT = ios, ERR = 904)
1369904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' )
1370      IF(lwm) WRITE ( numond, nambdy_index )
1371     
1372      SELECT CASE ( TRIM(ctypebdy) )
1373      CASE( 'N' )
1374         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1375            nbdyind  = Nj0glo - 2  ! set boundary to whole side of model domain.
1376            nbdybeg  = 2
1377            nbdyend  = Ni0glo - 1
1378         ENDIF
1379         nbdysegn = nbdysegn + 1
1380         npckgn(nbdysegn) = kb_bdy ! Save bdy package number
1381         jpjnob(nbdysegn) = nbdyind 
1382         jpindt(nbdysegn) = nbdybeg
1383         jpinft(nbdysegn) = nbdyend
1384         !
1385      CASE( 'S' )
1386         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1387            nbdyind  = 2           ! set boundary to whole side of model domain.
1388            nbdybeg  = 2
1389            nbdyend  = Ni0glo - 1
1390         ENDIF
1391         nbdysegs = nbdysegs + 1
1392         npckgs(nbdysegs) = kb_bdy ! Save bdy package number
1393         jpjsob(nbdysegs) = nbdyind
1394         jpisdt(nbdysegs) = nbdybeg
1395         jpisft(nbdysegs) = nbdyend
1396         !
1397      CASE( 'E' )
1398         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1399            nbdyind  = Ni0glo - 2  ! set boundary to whole side of model domain.
1400            nbdybeg  = 2
1401            nbdyend  = Nj0glo - 1
1402         ENDIF
1403         nbdysege = nbdysege + 1 
1404         npckge(nbdysege) = kb_bdy ! Save bdy package number
1405         jpieob(nbdysege) = nbdyind
1406         jpjedt(nbdysege) = nbdybeg
1407         jpjeft(nbdysege) = nbdyend
1408         !
1409      CASE( 'W' )
1410         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1411            nbdyind  = 2           ! set boundary to whole side of model domain.
1412            nbdybeg  = 2
1413            nbdyend  = Nj0glo - 1
1414         ENDIF
1415         nbdysegw = nbdysegw + 1
1416         npckgw(nbdysegw) = kb_bdy ! Save bdy package number
1417         jpiwob(nbdysegw) = nbdyind
1418         jpjwdt(nbdysegw) = nbdybeg
1419         jpjwft(nbdysegw) = nbdyend
1420         !
1421      CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
1422      END SELECT
1423     
1424      ! For simplicity we assume that in case of straight bdy, arrays have the same length
1425      ! (even if it is true that last tangential velocity points
1426      ! are useless). This simplifies a little bit boundary data format (and agrees with format
1427      ! used so far in obc package)
1428     
1429      knblendta(1:jpbgrd) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy)
1430     
1431   END SUBROUTINE bdy_read_seg
1432
1433   
1434   SUBROUTINE bdy_ctl_seg
1435      !!----------------------------------------------------------------------
1436      !!                 ***  ROUTINE bdy_ctl_seg  ***
1437      !!
1438      !! ** Purpose :   Check straight open boundary segments location
1439      !!
1440      !! ** Method  :   - Look for open boundary corners
1441      !!                - Check that segments start or end on land
1442      !!----------------------------------------------------------------------
1443      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest 
1444      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns 
1445      REAL(wp), DIMENSION(2) ::   ztestmask
1446      !!----------------------------------------------------------------------
1447      !
1448      IF (lwp) WRITE(numout,*) ' '
1449      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1450      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1451      !
1452      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege
1453      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw
1454      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn
1455      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs
1456      !
1457      ! 1. Check bounds
1458      !----------------
1459      DO ib = 1, nbdysegn
1460         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
1461         IF ((jpjnob(ib).ge.Nj0glo-1).or.& 
1462            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1463         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1464         IF (jpindt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1465         IF (jpinft(ib).gt.Ni0glo)     CALL ctl_stop( 'End index out of domain' )
1466      END DO
1467      !
1468      DO ib = 1, nbdysegs
1469         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
1470         IF ((jpjsob(ib).ge.Nj0glo-1).or.& 
1471            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1472         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1473         IF (jpisdt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1474         IF (jpisft(ib).gt.Ni0glo)     CALL ctl_stop( 'End index out of domain' )
1475      END DO
1476      !
1477      DO ib = 1, nbdysege
1478         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib)
1479         IF ((jpieob(ib).ge.Ni0glo-1).or.& 
1480            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1481         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1482         IF (jpjedt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1483         IF (jpjeft(ib).gt.Nj0glo)     CALL ctl_stop( 'End index out of domain' )
1484      END DO
1485      !
1486      DO ib = 1, nbdysegw
1487         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib)
1488         IF ((jpiwob(ib).ge.Ni0glo-1).or.& 
1489            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1490         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1491         IF (jpjwdt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1492         IF (jpjwft(ib).gt.Nj0glo)     CALL ctl_stop( 'End index out of domain' )
1493      ENDDO
1494      !     
1495      ! 2. Look for segment crossings
1496      !------------------------------
1497      IF (lwp) WRITE(numout,*) '**Look for segments corners  :'
1498      !
1499      itest = 0 ! corner number
1500      !
1501      ! flag to detect if start or end of open boundary belongs to a corner
1502      ! if not (=0), it must be on land.
1503      ! if a corner is detected, save bdy package number for further tests
1504      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1505      ! South/West crossings
1506      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1507         DO ib1 = 1, nbdysegw       
1508            DO ib2 = 1, nbdysegs
1509               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1510                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1511                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1512                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1513                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN 
1514                     ! We have a possible South-West corner                     
1515!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1516!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1517                     icornw(ib1,1) = npckgs(ib2)
1518                     icorns(ib2,1) = npckgw(ib1)
1519                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1520                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1521                        &                                     jpisft(ib2), jpjwft(ib1)
1522                     WRITE(ctmp2,*) ' Not allowed yet'
1523                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 
1524                        &                            ' and South segment: ',npckgs(ib2)
1525                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1526                  ELSE
1527                     WRITE(ctmp1,*) ' Check South and West Open boundary indices'
1528                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , &
1529                        &                            ' and South segment: ',npckgs(ib2)
1530                     CALL ctl_stop( ctmp1, ctmp2 )
1531                  END IF
1532               END IF
1533            END DO
1534         END DO
1535      END IF
1536      !
1537      ! South/East crossings
1538      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1539         DO ib1 = 1, nbdysege
1540            DO ib2 = 1, nbdysegs
1541               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1542                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1543                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. &
1544                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN
1545                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1546                     ! We have a possible South-East corner
1547!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1548!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1549                     icorne(ib1,1) = npckgs(ib2)
1550                     icorns(ib2,2) = npckge(ib1)
1551                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1552                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1553                        &                                     jpisdt(ib2), jpjeft(ib1)
1554                     WRITE(ctmp2,*) ' Not allowed yet'
1555                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1556                        &                            ' and South segment: ',npckgs(ib2)
1557                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1558                  ELSE
1559                     WRITE(ctmp1,*) ' Check South and East Open boundary indices'
1560                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1561                     &                               ' and South segment: ',npckgs(ib2)
1562                     CALL ctl_stop( ctmp1, ctmp2 )
1563                  END IF
1564               END IF
1565            END DO
1566         END DO
1567      END IF
1568      !
1569      ! North/West crossings
1570      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1571         DO ib1 = 1, nbdysegw       
1572            DO ib2 = 1, nbdysegn
1573               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. &
1574                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. &
1575                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1576                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1577                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1578                     ! We have a possible North-West corner
1579!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1580!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1581                     icornw(ib1,2) = npckgn(ib2)
1582                     icornn(ib2,1) = npckgw(ib1)
1583                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1584                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1585                     &                                     jpinft(ib2), jpjwdt(ib1)
1586                     WRITE(ctmp2,*) ' Not allowed yet'
1587                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1588                     &                               ' and North segment: ',npckgn(ib2)
1589                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1590                  ELSE
1591                     WRITE(ctmp1,*) ' Check North and West Open boundary indices'
1592                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1593                     &                               ' and North segment: ',npckgn(ib2)
1594                     CALL ctl_stop( ctmp1, ctmp2 )
1595                  END IF
1596               END IF
1597            END DO
1598         END DO
1599      END IF
1600      !
1601      ! North/East crossings
1602      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1603         DO ib1 = 1, nbdysege       
1604            DO ib2 = 1, nbdysegn
1605               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1606                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1607                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1608                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1609                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1610                     ! We have a possible North-East corner
1611!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1612!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1613                     icorne(ib1,2) = npckgn(ib2)
1614                     icornn(ib2,2) = npckge(ib1)
1615                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1616                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1617                     &                                     jpindt(ib2), jpjedt(ib1)
1618                     WRITE(ctmp2,*) ' Not allowed yet'
1619                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1620                     &                               ' and North segment: ',npckgn(ib2)
1621                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1622                  ELSE
1623                     WRITE(ctmp1,*) ' Check North and East Open boundary indices'
1624                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1625                     &                               ' and North segment: ',npckgn(ib2)
1626                     CALL ctl_stop( ctmp1, ctmp2 )
1627                  END IF
1628               END IF
1629            END DO
1630         END DO
1631      END IF
1632      !
1633      ! 3. Check if segment extremities are on land
1634      !--------------------------------------------
1635      !
1636      ! West segments
1637      DO ib = 1, nbdysegw
1638         ! get mask at boundary extremities:
1639         ztestmask(1:2)=0.
1640         DO ji = 1, jpi
1641            DO jj = 1, jpj             
1642              IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1643              IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1644            END DO
1645         END DO
1646         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1647
1648         IF (ztestmask(1)==1) THEN
1649            IF (icornw(ib,1)==0) THEN
1650               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1651               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1652            ELSE
1653               ! This is a corner
1654               IF(lwp) WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1655               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1656               itest=itest+1
1657            ENDIF
1658         ENDIF
1659         IF (ztestmask(2)==1) THEN
1660            IF (icornw(ib,2)==0) THEN
1661               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1662               CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' )
1663            ELSE
1664               ! This is a corner
1665               IF(lwp) WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1666               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1667               itest=itest+1
1668            ENDIF
1669         ENDIF
1670      END DO
1671      !
1672      ! East segments
1673      DO ib = 1, nbdysege
1674         ! get mask at boundary extremities:
1675         ztestmask(1:2)=0.
1676         DO ji = 1, jpi
1677            DO jj = 1, jpj             
1678              IF( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjedt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1679              IF( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjeft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1680            END DO
1681         END DO
1682         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1683
1684         IF (ztestmask(1)==1) THEN
1685            IF (icorne(ib,1)==0) THEN
1686               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1687               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1688            ELSE
1689               ! This is a corner
1690               IF(lwp) WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1691               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1692               itest=itest+1
1693            ENDIF
1694         ENDIF
1695         IF (ztestmask(2)==1) THEN
1696            IF (icorne(ib,2)==0) THEN
1697               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1698               CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1699            ELSE
1700               ! This is a corner
1701               IF(lwp) WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1702               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1703               itest=itest+1
1704            ENDIF
1705         ENDIF
1706      END DO
1707      !
1708      ! South segments
1709      DO ib = 1, nbdysegs
1710         ! get mask at boundary extremities:
1711         ztestmask(1:2)=0.
1712         DO ji = 1, jpi
1713            DO jj = 1, jpj             
1714              IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1715              IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1716            END DO
1717         END DO
1718         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1719
1720         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1721            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1722            CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1723         ENDIF
1724         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1725            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1726            CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1727         ENDIF
1728      END DO
1729      !
1730      ! North segments
1731      DO ib = 1, nbdysegn
1732         ! get mask at boundary extremities:
1733         ztestmask(1:2)=0.
1734         DO ji = 1, jpi
1735            DO jj = 1, jpj             
1736               IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpindt(ib) )   ztestmask(1) = tmask(ji,jj,1)
1737               IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpinft(ib) )   ztestmask(2) = tmask(ji,jj,1) 
1738            END DO
1739         END DO
1740         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1741
1742         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1743            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1744            CALL ctl_stop( ctmp1, ' does not start on land' )
1745         ENDIF
1746         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1747            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1748            CALL ctl_stop( ctmp1, ' does not end on land' )
1749         ENDIF
1750      END DO
1751      !
1752      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1753      !
1754      ! Other tests TBD:
1755      ! segments completly on land
1756      ! optimized open boundary array length according to landmask
1757      ! Nudging layers that overlap with interior domain
1758      !
1759   END SUBROUTINE bdy_ctl_seg
1760
1761   
1762   SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta ) 
1763      !!----------------------------------------------------------------------
1764      !!                 ***  ROUTINE bdy_coords_seg  ***
1765      !!
1766      !! ** Purpose :  build nbidta, nbidta, nbrdta for bdy built with segments
1767      !!
1768      !! ** Method  : 
1769      !!
1770      !!----------------------------------------------------------------------
1771      INTEGER, DIMENSION(:,:,:), intent(  out)  ::   nbidta, nbjdta, nbrdta   ! Index arrays: i and j indices of bdy dta
1772      !!
1773      INTEGER  ::   ii, ij, ir, iseg
1774      INTEGER  ::   igrd         ! grid type (t=1, u=2, v=3)
1775      INTEGER  ::   icount       !
1776      INTEGER  ::   ib_bdy       ! bdy number
1777      !!----------------------------------------------------------------------
1778
1779      ! East
1780      !-----
1781      DO iseg = 1, nbdysege
1782         ib_bdy = npckge(iseg)
1783         !
1784         ! ------------ T points -------------
1785         igrd=1
1786         icount=0
1787         DO ir = 1, nn_rimwidth(ib_bdy)
1788            DO ij = jpjedt(iseg), jpjeft(iseg)
1789               icount = icount + 1
1790               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir + nn_hls
1791               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1792               nbrdta(icount, igrd, ib_bdy) = ir
1793            ENDDO
1794         ENDDO
1795         !
1796         ! ------------ U points -------------
1797         igrd=2
1798         icount=0
1799         DO ir = 1, nn_rimwidth(ib_bdy)
1800            DO ij = jpjedt(iseg), jpjeft(iseg)
1801               icount = icount + 1
1802               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir + nn_hls
1803               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1804               nbrdta(icount, igrd, ib_bdy) = ir
1805            ENDDO
1806         ENDDO
1807         !
1808         ! ------------ V points -------------
1809         igrd=3
1810         icount=0
1811         DO ir = 1, nn_rimwidth(ib_bdy)
1812            !            DO ij = jpjedt(iseg), jpjeft(iseg) - 1
1813            DO ij = jpjedt(iseg), jpjeft(iseg)
1814               icount = icount + 1
1815               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir + nn_hls
1816               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1817               nbrdta(icount, igrd, ib_bdy) = ir
1818            ENDDO
1819            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1820            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1821         ENDDO
1822      ENDDO
1823      !
1824      ! West
1825      !-----
1826      DO iseg = 1, nbdysegw
1827         ib_bdy = npckgw(iseg)
1828         !
1829         ! ------------ T points -------------
1830         igrd=1
1831         icount=0
1832         DO ir = 1, nn_rimwidth(ib_bdy)
1833            DO ij = jpjwdt(iseg), jpjwft(iseg)
1834               icount = icount + 1
1835               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1836               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1837               nbrdta(icount, igrd, ib_bdy) = ir
1838            ENDDO
1839         ENDDO
1840         !
1841         ! ------------ U points -------------
1842         igrd=2
1843         icount=0
1844         DO ir = 1, nn_rimwidth(ib_bdy)
1845            DO ij = jpjwdt(iseg), jpjwft(iseg)
1846               icount = icount + 1
1847               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1848               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1849               nbrdta(icount, igrd, ib_bdy) = ir
1850            ENDDO
1851         ENDDO
1852         !
1853         ! ------------ V points -------------
1854         igrd=3
1855         icount=0
1856         DO ir = 1, nn_rimwidth(ib_bdy)
1857            !            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
1858            DO ij = jpjwdt(iseg), jpjwft(iseg)
1859               icount = icount + 1
1860               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1861               nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1862               nbrdta(icount, igrd, ib_bdy) = ir
1863            ENDDO
1864            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1865            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1866         ENDDO
1867      ENDDO
1868      !
1869      ! North
1870      !-----
1871      DO iseg = 1, nbdysegn
1872         ib_bdy = npckgn(iseg)
1873         !
1874         ! ------------ T points -------------
1875         igrd=1
1876         icount=0
1877         DO ir = 1, nn_rimwidth(ib_bdy)
1878            DO ii = jpindt(iseg), jpinft(iseg)
1879               icount = icount + 1
1880               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1881               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir + nn_hls 
1882               nbrdta(icount, igrd, ib_bdy) = ir
1883            ENDDO
1884         ENDDO
1885         !
1886         ! ------------ U points -------------
1887         igrd=2
1888         icount=0
1889         DO ir = 1, nn_rimwidth(ib_bdy)
1890            !            DO ii = jpindt(iseg), jpinft(iseg) - 1
1891            DO ii = jpindt(iseg), jpinft(iseg)
1892               icount = icount + 1
1893               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1894               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir + nn_hls
1895               nbrdta(icount, igrd, ib_bdy) = ir
1896            ENDDO
1897            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1898            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1899         ENDDO
1900         !
1901         ! ------------ V points -------------
1902         igrd=3
1903         icount=0
1904         DO ir = 1, nn_rimwidth(ib_bdy)
1905            DO ii = jpindt(iseg), jpinft(iseg)
1906               icount = icount + 1
1907               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1908               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir + nn_hls
1909               nbrdta(icount, igrd, ib_bdy) = ir
1910            ENDDO
1911         ENDDO
1912      ENDDO
1913      !
1914      ! South
1915      !-----
1916      DO iseg = 1, nbdysegs
1917         ib_bdy = npckgs(iseg)
1918         !
1919         ! ------------ T points -------------
1920         igrd=1
1921         icount=0
1922         DO ir = 1, nn_rimwidth(ib_bdy)
1923            DO ii = jpisdt(iseg), jpisft(iseg)
1924               icount = icount + 1
1925               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1926               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1927               nbrdta(icount, igrd, ib_bdy) = ir
1928            ENDDO
1929         ENDDO
1930         !
1931         ! ------------ U points -------------
1932         igrd=2
1933         icount=0
1934         DO ir = 1, nn_rimwidth(ib_bdy)
1935            !            DO ii = jpisdt(iseg), jpisft(iseg) - 1
1936            DO ii = jpisdt(iseg), jpisft(iseg)
1937               icount = icount + 1
1938               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1939               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1940               nbrdta(icount, igrd, ib_bdy) = ir
1941            ENDDO
1942            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1943            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1944         ENDDO
1945         !
1946         ! ------------ V points -------------
1947         igrd=3
1948         icount=0
1949         DO ir = 1, nn_rimwidth(ib_bdy)
1950            DO ii = jpisdt(iseg), jpisft(iseg)
1951               icount = icount + 1
1952               nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1953               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1954               nbrdta(icount, igrd, ib_bdy) = ir
1955            ENDDO
1956         ENDDO
1957      ENDDO
1958
1959     
1960   END SUBROUTINE bdy_coords_seg
1961   
1962   
1963   SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1964      !!----------------------------------------------------------------------
1965      !!                 ***  ROUTINE bdy_ctl_corn  ***
1966      !!
1967      !! ** Purpose :   Check numerical schemes consistency between
1968      !!                segments having a common corner
1969      !!
1970      !! ** Method  :   
1971      !!----------------------------------------------------------------------
1972      INTEGER, INTENT(in)  ::   ib1, ib2
1973      INTEGER :: itest
1974      !!----------------------------------------------------------------------
1975      itest = 0
1976
1977      IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) )   itest = itest + 1
1978      IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) )   itest = itest + 1
1979      IF( cn_tra  (ib1) /= cn_tra  (ib2) )   itest = itest + 1
1980      !
1981      IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) )   itest = itest + 1
1982      IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) )   itest = itest + 1
1983      IF( nn_tra_dta  (ib1) /= nn_tra_dta  (ib2) )   itest = itest + 1
1984      !
1985      IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) )   itest = itest + 1   
1986      !
1987      IF( itest>0 ) THEN
1988         WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2
1989         CALL ctl_stop( ctmp1, ' have different open bdy schemes' )
1990      ENDIF
1991      !
1992   END SUBROUTINE bdy_ctl_corn
1993
1994
1995   SUBROUTINE bdy_meshwri()
1996      !!----------------------------------------------------------------------
1997      !!                 ***  ROUTINE bdy_meshwri  ***
1998      !!         
1999      !! ** Purpose :   write netcdf file with nbr, flagu, flagv, ntreat for T, U
2000      !!                and V points in 2D arrays for easier visualisation/control
2001      !!
2002      !! ** Method  :   use iom_rstput as in domwri.F
2003      !!----------------------------------------------------------------------     
2004      INTEGER  ::   ib_bdy, ii, ij, igrd, ib     ! dummy loop indices
2005      INTEGER  ::   inum                                   !   -       -
2006      REAL(wp), POINTER, DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields
2007      REAL(wp)         , DIMENSION(jpi,jpj) ::   ztmp
2008      CHARACTER(LEN=1) , DIMENSION(jpbgrd)  ::   cgrid
2009      !!----------------------------------------------------------------------     
2010      cgrid = (/'t','u','v'/)
2011      CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. )
2012      DO igrd = 1, jpbgrd
2013         SELECT CASE( igrd )
2014         CASE( 1 )   ;   zmask => tmask(:,:,1)
2015         CASE( 2 )   ;   zmask => umask(:,:,1)
2016         CASE( 3 )   ;   zmask => vmask(:,:,1)
2017         END SELECT
2018         ztmp(:,:) = zmask(:,:)
2019         DO ib_bdy = 1, nb_bdy
2020            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)      ! nbr deined for all rims
2021               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
2022               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
2023               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10.
2024               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
2025            END DO
2026         END DO
2027         CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
2028         ztmp(:,:) = zmask(:,:)
2029         DO ib_bdy = 1, nb_bdy
2030            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagu defined only for rims 0 and 1
2031               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
2032               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
2033               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10.
2034               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
2035            END DO
2036         END DO
2037         CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
2038         ztmp(:,:) = zmask(:,:)
2039         DO ib_bdy = 1, nb_bdy
2040            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagv defined only for rims 0 and 1
2041               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
2042               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
2043               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10.
2044               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
2045            END DO
2046         END DO
2047         CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
2048         ztmp(:,:) = zmask(:,:)
2049         DO ib_bdy = 1, nb_bdy
2050            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! ntreat defined only for rims 0 and 1
2051               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
2052               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
2053               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10.
2054               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
2055            END DO
2056         END DO
2057         CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
2058      END DO
2059      CALL iom_close( inum )
2060
2061   END SUBROUTINE bdy_meshwri
2062   
2063   !!=================================================================================
2064END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.