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

Last change on this file since 11352 was 11352, checked in by smasson, 2 years ago

dev_r10984_HPC-13 : small bugfixes

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