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/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY – NEMO

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90 @ 11223

Last change on this file since 11223 was 11223, checked in by smasson, 5 years ago

dev_r10984_HPC-13 : cleaning of rewriting of bdydta, see #2285

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