source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/OCE/BDY/bdyini.F90 @ 13565

Last change on this file since 13565 was 13565, checked in by jchanut, 5 months ago

#2222, 1) Added parent bathymetry volume consistency check 2) Added velocity extrapolation in update 3) Corrected bdy issue #2519

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