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

Last change on this file since 11258 was 11258, checked in by smasson, 15 months ago

dev_r10984_HPC-13 : minor bugfix and cleaning, see #2285

  • 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', lwp )
81      ! make sur that all elements of the namelist variables have a default definition from namelist_ref
82      ln_coords_file (2:jp_bdy) = ln_coords_file (1)
83      cn_coords_file (2:jp_bdy) = cn_coords_file (1)
84      cn_dyn2d       (2:jp_bdy) = cn_dyn2d       (1)
85      nn_dyn2d_dta   (2:jp_bdy) = nn_dyn2d_dta   (1)
86      cn_dyn3d       (2:jp_bdy) = cn_dyn3d       (1)
87      nn_dyn3d_dta   (2:jp_bdy) = nn_dyn3d_dta   (1)
88      cn_tra         (2:jp_bdy) = cn_tra         (1)
89      nn_tra_dta     (2:jp_bdy) = nn_tra_dta     (1)   
90      ln_tra_dmp     (2:jp_bdy) = ln_tra_dmp     (1)
91      ln_dyn3d_dmp   (2:jp_bdy) = ln_dyn3d_dmp   (1)
92      rn_time_dmp    (2:jp_bdy) = rn_time_dmp    (1)
93      rn_time_dmp_out(2:jp_bdy) = rn_time_dmp_out(1)
94      cn_ice         (2:jp_bdy) = cn_ice         (1)
95      nn_ice_dta     (2:jp_bdy) = nn_ice_dta     (1)
96      rn_ice_tem     (2:jp_bdy) = rn_ice_tem     (1)
97      rn_ice_sal     (2:jp_bdy) = rn_ice_sal     (1)
98      rn_ice_age     (2:jp_bdy) = rn_ice_age     (1)
99      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
100      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
101902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
102      IF(lwm) WRITE ( numond, nambdy )
103
104      IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE.   ! forced for Agrif children
105
106      IF( nb_bdy == 0 ) ln_bdy = .FALSE.
107     
108      ! -----------------------------------------
109      ! unstructured open boundaries use control
110      ! -----------------------------------------
111      IF ( ln_bdy ) THEN
112         IF(lwp) WRITE(numout,*)
113         IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
114         IF(lwp) WRITE(numout,*) '~~~~~~~~'
115         !
116         ! Open boundaries definition (arrays and masks)
117         CALL bdy_def
118         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            ! Now look for crossings in user (namelist) defined open boundary segments:
389            IF( nbdysege > 0 .OR. nbdysegw > 0 .OR. nbdysegn > 0 .OR. nbdysegs > 0)   CALL bdy_ctl_seg
390
391         ELSE                                        ! Read size of arrays in boundary coordinates file.
392           
393            CALL iom_open( cn_coords_file(ib_bdy), inum )
394            DO igrd = 1, jpbgrd
395               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 
396               nblendta(igrd,ib_bdy) = MAXVAL(kdimsz)
397            END DO
398            CALL iom_close( inum )
399         ENDIF
400         !
401      END DO ! ib_bdy
402     
403      ! Allocate arrays
404      !---------------
405      jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
406      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), nbrdta(jpbdta, jpbgrd, nb_bdy) )
407   
408      ! Calculate global boundary index arrays or read in from file
409      !------------------------------------------------------------               
410      ! 1. Read global index arrays from boundary coordinates file.
411      DO ib_bdy = 1, nb_bdy
412         !
413         IF( ln_coords_file(ib_bdy) ) THEN
414            !
415            ALLOCATE( zz_read( MAXVAL(nblendta), 1 ) )         
416            CALL iom_open( cn_coords_file(ib_bdy), inum )
417            !
418            DO igrd = 1, jpbgrd
419               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
420               DO ii = 1,nblendta(igrd,ib_bdy)
421                  nbidta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
422               END DO
423               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
424               DO ii = 1,nblendta(igrd,ib_bdy)
425                  nbjdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
426               END DO
427               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
428               DO ii = 1,nblendta(igrd,ib_bdy)
429                  nbrdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
430               END DO
431               !
432               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
433               IF(lwp) WRITE(numout,*)
434               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
435               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
436               IF (ibr_max < nn_rimwidth(ib_bdy))   &
437                  CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
438            END DO
439            !
440            CALL iom_close( inum )
441            DEALLOCATE( zz_read )
442            !
443         ENDIF
444         !
445      END DO
446
447      ! 2. Now fill indices corresponding to straight open boundary arrays:
448      CALL bdy_coords_seg( nbidta, nbjdta, nbrdta )
449
450      !  Deal with duplicated points
451      !-----------------------------
452      ! We assign negative indices to duplicated points (to remove them from bdy points to be updated)
453      ! if their distance to the bdy is greater than the other
454      ! If their distance are the same, just keep only one to avoid updating a point twice
455      DO igrd = 1, jpbgrd
456         DO ib_bdy1 = 1, nb_bdy
457            DO ib_bdy2 = 1, nb_bdy
458               IF (ib_bdy1/=ib_bdy2) THEN
459                  DO ib1 = 1, nblendta(igrd,ib_bdy1)
460                     DO ib2 = 1, nblendta(igrd,ib_bdy2)
461                        IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
462                           &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
463                           !                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &
464                           !                                                       &              nbidta(ib1, igrd, ib_bdy1),      &
465                           !                                                       &              nbjdta(ib2, igrd, ib_bdy2)
466                           ! keep only points with the lowest distance to boundary:
467                           IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN
468                              nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
469                              nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
470                           ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN
471                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
472                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
473                              ! Arbitrary choice if distances are the same:
474                           ELSE
475                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
476                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
477                           ENDIF
478                        END IF
479                     END DO
480                  END DO
481               ENDIF
482            END DO
483         END DO
484      END DO
485      !
486      ! Find lenght of boundaries and rim on local mpi domain
487      !------------------------------------------------------
488      !
489      iwe = mig(1)
490      ies = mig(jpi)
491      iso = mjg(1) 
492      ino = mjg(jpj) 
493      !
494      DO ib_bdy = 1, nb_bdy
495         DO igrd = 1, jpbgrd
496            icount   = 0   ! initialization of local bdy length
497            icountr  = 0   ! initialization of local rim 0 and rim 1 bdy length
498            icountr0 = 0   ! initialization of local rim 0 bdy length
499            idx_bdy(ib_bdy)%nblen(igrd)     = 0
500            idx_bdy(ib_bdy)%nblenrim(igrd)  = 0
501            idx_bdy(ib_bdy)%nblenrim0(igrd) = 0
502            DO ib = 1, nblendta(igrd,ib_bdy)
503               ! check that data is in correct order in file
504               IF( ib > 1 ) THEN
505                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ib-1,igrd,ib_bdy) ) THEN
506                     CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', &
507                        &        ' in order of distance from edge nbr A utility for re-ordering ', &
508                        &        ' boundary coordinates and data files exists in the TOOLS/OBC directory')
509                  ENDIF
510               ENDIF
511               ! check if point is in local domain
512               IF(  nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND.   &
513                  & nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino      ) THEN
514                  !
515                  icount = icount + 1
516                  IF( nbrdta(ib,igrd,ib_bdy) == 1 .OR. nbrdta(ib,igrd,ib_bdy) == 0 )   icountr = icountr + 1
517                  IF( nbrdta(ib,igrd,ib_bdy) == 0 )   icountr0 = icountr0 + 1
518               ENDIF
519            END DO
520            idx_bdy(ib_bdy)%nblen    (igrd) = icount   !: length of boundary data on each proc
521            idx_bdy(ib_bdy)%nblenrim (igrd) = icountr  !: length of rim 0 and rim 1 boundary data on each proc   
522            idx_bdy(ib_bdy)%nblenrim0(igrd) = icountr0 !: length of rim 0 boundary data on each proc     
523         END DO   ! igrd
524
525         ! Allocate index arrays for this boundary set
526         !--------------------------------------------
527         ilen1 = MAXVAL( idx_bdy(ib_bdy)%nblen(:) )
528         ALLOCATE( idx_bdy(ib_bdy)%nbi   (ilen1,jpbgrd) ,   &
529            &      idx_bdy(ib_bdy)%nbj   (ilen1,jpbgrd) ,   &
530            &      idx_bdy(ib_bdy)%nbr   (ilen1,jpbgrd) ,   &
531            &      idx_bdy(ib_bdy)%nbd   (ilen1,jpbgrd) ,   &
532            &      idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) ,   &
533            &      idx_bdy(ib_bdy)%ntreat(ilen1,jpbgrd) ,   &
534            &      idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) ,   &
535            &      idx_bdy(ib_bdy)%nbw   (ilen1,jpbgrd) ,   &
536            &      idx_bdy(ib_bdy)%flagu (ilen1,jpbgrd) ,   &
537            &      idx_bdy(ib_bdy)%flagv (ilen1,jpbgrd) )
538
539         ! Dispatch mapping indices and discrete distances on each processor
540         ! -----------------------------------------------------------------
541         DO igrd = 1, jpbgrd
542            icount  = 0
543            ! Outer loop on rimwidth to ensure outermost points come first in the local arrays.
544            DO ir = 0, nn_rimwidth(ib_bdy)
545               DO ib = 1, nblendta(igrd,ib_bdy)
546                  ! check if point is in local domain and equals ir
547                  IF(  nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND.   &
548                     & nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino .AND.   &
549                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
550                     !
551                     icount = icount  + 1
552                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1   ! global to local indexes
553                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1   ! global to local indexes
554                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
555                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
556                  ENDIF
557               END DO
558            END DO
559         END DO   ! igrd
560
561      END DO   ! ib_bdy
562
563      ! Initialize array indicating communications in bdy
564      ! -------------------------------------------------
565      ALLOCATE( lsend_bdy(nb_bdy,jpbgrd,4,0:1), lrecv_bdy(nb_bdy,jpbgrd,4,0:1) )
566      lsend_bdy(:,:,:,:) = .false.
567      lrecv_bdy(:,:,:,:) = .false. 
568
569      DO ib_bdy = 1, nb_bdy
570         DO igrd = 1, jpbgrd
571            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! only the rim triggers communications, see bdy routines
572               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
573               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
574               IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN   ;   ir = 0
575               ELSE                                                 ;   ir = 1
576               END IF
577               !
578               ! check if point has to be sent     to   a neighbour
579               ! W neighbour and on the inner left  side
580               IF( ii == 2     .and. (nbondi == 0 .or. nbondi ==  1) )   lsend_bdy(ib_bdy,igrd,1,ir) = .true.
581               ! E neighbour and on the inner right side
582               IF( ii == jpi-1 .and. (nbondi == 0 .or. nbondi == -1) )   lsend_bdy(ib_bdy,igrd,2,ir) = .true.
583               ! S neighbour and on the inner down side
584               IF( ij == 2     .and. (nbondj == 0 .or. nbondj ==  1) )   lsend_bdy(ib_bdy,igrd,3,ir) = .true.
585               ! N neighbour and on the inner up   side
586               IF( ij == jpj-1 .and. (nbondj == 0 .or. nbondj == -1) )   lsend_bdy(ib_bdy,igrd,4,ir) = .true.
587               !
588               ! check if point has to be received from a neighbour
589               ! W neighbour and on the outter left  side
590               IF( ii == 1     .and. (nbondi == 0 .or. nbondi ==  1) )   lrecv_bdy(ib_bdy,igrd,1,ir) = .true.
591               ! E neighbour and on the outter right side
592               IF( ii == jpi   .and. (nbondi == 0 .or. nbondi == -1) )   lrecv_bdy(ib_bdy,igrd,2,ir) = .true.
593               ! S neighbour and on the outter down side
594               IF( ij == 1     .and. (nbondj == 0 .or. nbondj ==  1) )   lrecv_bdy(ib_bdy,igrd,3,ir) = .true.
595               ! N neighbour and on the outter up   side
596               IF( ij == jpj   .and. (nbondj == 0 .or. nbondj == -1) )   lrecv_bdy(ib_bdy,igrd,4,ir) = .true.
597               !
598            END DO
599         END DO  ! igrd
600
601         ! Compute rim weights for FRS scheme
602         ! ----------------------------------
603         DO igrd = 1, jpbgrd
604            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
605               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same weights
606               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 )      ! tanh formulation
607               !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic
608               !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy))       ! linear
609            END DO
610         END DO
611
612         ! Compute damping coefficients
613         ! ----------------------------
614         DO igrd = 1, jpbgrd
615            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
616               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same damping coefficients
617               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 
618                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
619               idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) & 
620                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
621            END DO
622         END DO
623
624      END DO ! ib_bdy
625
626      ! ------------------------------------------------------
627      ! Initialise masks and find normal/tangential directions
628      ! ------------------------------------------------------
629
630      ! ------------------------------------------
631      ! handle rim0, do as if rim 1 was free ocean
632      ! ------------------------------------------
633
634      ztmask(:,:) = tmask(:,:,1)   ;   zumask(:,:) = umask(:,:,1)   ;   zvmask(:,:) = vmask(:,:,1)
635      ! For the flagu/flagv calculation below we require a version of fmask without
636      ! the land boundary condition (shlat) included:
637      DO ij = 1, jpjm1
638         DO ii = 1, jpim1
639            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
640               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
641         END DO
642      END DO
643      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )
644
645      ! Read global 2D mask at T-points: bdytmask
646      ! -----------------------------------------
647      ! bdytmask = 1  on the computational domain AND on open boundaries
648      !          = 0  elsewhere   
649
650      bdytmask(:,:) = ssmask(:,:)
651
652      ! Derive mask on U and V grid from mask on T grid
653      DO ij = 1, jpjm1
654         DO ii = 1, jpim1
655            bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1,ij  )
656            bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii  ,ij+1) 
657         END DO
658      END DO
659      CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1., bdyvmask, 'V', 1. )   ! Lateral boundary cond.
660
661      ! bdy masks are now set to zero on rim 0 points:
662      DO ib_bdy = 1, nb_bdy
663         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0
664            bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
665         END DO
666         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2)   ! extent of rim 0
667            bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
668         END DO
669         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3)   ! extent of rim 0
670            bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
671         END DO
672      END DO
673
674      CALL bdy_rim_treat( zumask, zvmask, zfmask, .true. )   ! compute flagu, flagv, ntreat on rim 0
675
676      ! ------------------------------------
677      ! handle rim1, do as if rim 0 was land
678      ! ------------------------------------
679     
680      ! z[tuv]mask are now set to zero on rim 0 points:
681      DO ib_bdy = 1, nb_bdy
682         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0
683            ztmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
684         END DO
685         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2)   ! extent of rim 0
686            zumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
687         END DO
688         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3)   ! extent of rim 0
689            zvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
690         END DO
691      END DO
692
693      ! Recompute zfmask
694      DO ij = 1, jpjm1
695         DO ii = 1, jpim1
696            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
697               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
698         END DO
699      END DO
700      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )
701
702      ! bdy masks are now set to zero on rim1 points:
703      DO ib_bdy = 1, nb_bdy
704         DO ib = idx_bdy(ib_bdy)%nblenrim0(1) + 1,  idx_bdy(ib_bdy)%nblenrim(1)   ! extent of rim 1
705            bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
706         END DO
707         DO ib = idx_bdy(ib_bdy)%nblenrim0(2) + 1,  idx_bdy(ib_bdy)%nblenrim(2)   ! extent of rim 1
708            bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
709         END DO
710         DO ib = idx_bdy(ib_bdy)%nblenrim0(3) + 1,  idx_bdy(ib_bdy)%nblenrim(3)   ! extent of rim 1
711            bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
712         END DO
713      END DO
714
715      CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. )   ! compute flagu, flagv, ntreat on rim 1
716      !
717      ! Check which boundaries might need communication
718      ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,4,0:1), lrecv_bdyint(nb_bdy,jpbgrd,4,0:1) )
719      lsend_bdyint(:,:,:,:) = .false.
720      lrecv_bdyint(:,:,:,:) = .false. 
721      ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,4,0:1), lrecv_bdyext(nb_bdy,jpbgrd,4,0:1) )
722      lsend_bdyext(:,:,:,:) = .false.
723      lrecv_bdyext(:,:,:,:) = .false.
724      !
725      DO igrd = 1, jpbgrd
726         DO ib_bdy = 1, nb_bdy
727            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
728               IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE
729               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
730               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
731               ir = idx_bdy(ib_bdy)%nbr(ib,igrd)
732               flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd))
733               flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd))
734               iibe = ii - flagu   ! neighbouring point towards the exterior of the computational domain
735               ijbe = ij - flagv
736               iibi = ii + flagu   ! neighbouring point towards the interior of the computational domain
737               ijbi = ij + flagv
738               CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 )   ! free ocean neighbours
739               !
740               ! search neighbour in the  west/east  direction
741               ! Rim is on the halo and computed ocean is towards exterior of mpi domain 
742               !      <--    (o exterior)     --> 
743               ! (1)  o|x         OR    (2)   x|o
744               !       |___                 ___|
745               IF( iibi == 0     .OR. ii1 == 0     .OR. ii2 == 0     .OR. ii3 == 0     )   lrecv_bdyint(ib_bdy,igrd,1,ir) = .true.
746               IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 )   lrecv_bdyint(ib_bdy,igrd,2,ir) = .true. 
747               IF( iibe == 0                                                           )   lrecv_bdyext(ib_bdy,igrd,1,ir) = .true.
748               IF( iibe == jpi+1                                                       )   lrecv_bdyext(ib_bdy,igrd,2,ir) = .true. 
749               ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo
750               ! :¨¨¨¨¨|¨¨-->    |                                             |    <--¨¨|¨¨¨¨¨:
751               ! :     |  x:o    |    neighbour limited by ... would need o    |    o:x  |     :
752               ! :.....|_._:_____|   (1) W neighbour         E neighbour (2)   |_____:_._|.....:
753               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. &
754                  & ( iibi == 3     .OR. ii1 == 3     .OR. ii2 == 3     .OR. ii3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,1,ir)=.true.
755               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. &
756                  & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) )   lsend_bdyint(ib_bdy,igrd,2,ir)=.true.
757               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. iibe == 3     )   lsend_bdyext(ib_bdy,igrd,1,ir)=.true.
758               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 )   lsend_bdyext(ib_bdy,igrd,2,ir)=.true.
759               !
760               ! search neighbour in the north/south direction   
761               ! Rim is on the halo and computed ocean is towards exterior of mpi domain
762               !(3)   |       |         ^   ___o___     
763               !  |   |___x___|   OR    |  |   x   |
764               !  v       o           (4)  |       |
765               IF( ijbi == 0     .OR. ij1 == 0     .OR. ij2 == 0     .OR. ij3 == 0     )   lrecv_bdyint(ib_bdy,igrd,3,ir) = .true.
766               IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 )   lrecv_bdyint(ib_bdy,igrd,4,ir) = .true.
767               IF( ijbe == 0                                                           )   lrecv_bdyext(ib_bdy,igrd,3,ir) = .true.
768               IF( ijbe == jpj+1                                                       )   lrecv_bdyext(ib_bdy,igrd,4,ir) = .true.
769               ! Check if neighbour has its rim parallel to its mpi subdomain     _________  border and next to its halo
770               !   ^  |    o    |                                                :         :
771               !   |  |¨¨¨¨x¨¨¨¨|   neighbour limited by ... would need o     |  |....x....|
772               !      :_________:  (3) S neighbour          N neighbour (4)   v  |    o    |   
773               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. &
774                  & ( ijbi == 3     .OR. ij1 == 3     .OR. ij2 == 3     .OR. ij3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,3,ir)=.true.
775               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. &
776                  & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) )   lsend_bdyint(ib_bdy,igrd,4,ir)=.true.
777               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. ijbe == 3     )   lsend_bdyext(ib_bdy,igrd,3,ir)=.true.
778               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 )   lsend_bdyext(ib_bdy,igrd,4,ir)=.true.
779            END DO
780         END DO
781      END DO
782
783      DO ib_bdy = 1,nb_bdy
784         IF(  cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &
785            & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &
786            & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN
787            DO igrd = 1, jpbgrd
788               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
789                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
790                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
791                  IF(  mig(ii) > 2 .AND. mig(ii) < jpiglo-2 .AND. mjg(ij) > 2 .AND. mjg(ij) < jpjglo-2  ) THEN
792                     WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain'
793                     CALL ctl_stop( ctmp1 )
794                  END IF
795               END DO
796            END DO
797         END IF
798      END DO
799      !
800      DEALLOCATE( nbidta, nbjdta, nbrdta )
801      !
802   END SUBROUTINE bdy_def
803
804
805   SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 )
806      !!----------------------------------------------------------------------
807      !!                 ***  ROUTINE bdy_rim_treat  ***
808      !!
809      !! ** Purpose :   Initialize structures ( flagu, flagv, ntreat ) indicating how rim points
810      !!                  are to be handled in the boundary condition treatment
811      !!
812      !! ** Method  :   - to handle rim 0 zmasks must indicate ocean points      (set at one on rim 0 and rim 1 and interior)
813      !!                            and bdymasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
814      !!                    (as if rim 1 was free ocean)
815      !!                - to handle rim 1 zmasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
816      !!                            and bdymasks must indicate free ocean points (set at one on interior)
817      !!                    (as if rim 0 was land)
818      !!                - we can then check in which direction the interior of the computational domain is with the difference
819      !!                         mask array values on both sides to compute flagu and flagv
820      !!                - and look at the ocean neighbours to compute ntreat
821      !!----------------------------------------------------------------------
822      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pfmask   ! temporary fmask excluding coastal boundary condition (shlat)
823      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pumask, pvmask   ! temporary t/u/v mask array
824      LOGICAL                             , INTENT (in   ) :: lrim0    ! .true. -> rim 0   .false. -> rim 1
825      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, icount       ! dummy loop indices
826      INTEGER  ::   i_offset, j_offset, inn                ! local integer
827      INTEGER  ::   ibeg, iend                             ! local integer
828      LOGICAL  ::   llnon, llson, llean, llwen             ! local logicals indicating the presence of a ocean neighbour
829      REAL(wp), POINTER, DIMENSION(:,:)       ::   zmask   ! pointer to 2D mask fields
830      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
831      CHARACTER(LEN=1), DIMENSION(jpbgrd)     ::   cgrid
832      REAL(wp)        , DIMENSION(jpi,jpj)    ::   ztmp
833      !!----------------------------------------------------------------------
834
835      cgrid = (/'t','u','v'/)
836
837      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
838
839         ! Calculate relationship of U direction to the local orientation of the boundary
840         ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
841         ! flagu =  0 : u is tangential
842         ! flagu =  1 : u is normal to the boundary and is direction is inward
843         DO igrd = 1, jpbgrd 
844            SELECT CASE( igrd )
845               CASE( 1 )   ;   zmask => pumask     ;   i_offset = 0
846               CASE( 2 )   ;   zmask => bdytmask   ;   i_offset = 1
847               CASE( 3 )   ;   zmask => pfmask     ;   i_offset = 0
848            END SELECT
849            icount = 0
850            ztmp(:,:) = -999._wp
851            IF( lrim0 ) THEN   ! extent of rim 0
852               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
853            ELSE               ! extent of rim 1
854               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
855            END IF
856            DO ib = ibeg, iend 
857               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
858               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
859               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE
860               zwfl = zmask(ii+i_offset-1,ij)
861               zefl = zmask(ii+i_offset  ,ij)
862               ! This error check only works if you are using the bdyXmask arrays
863               IF( i_offset == 1 .and. zefl + zwfl == 2. ) THEN
864                  icount = icount + 1
865                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
866               ELSE
867                  ztmp(ii,ij) = -zwfl + zefl
868               ENDIF
869            END DO
870            IF( icount /= 0 ) THEN
871               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
872                  ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
873               CALL ctl_stop( ctmp1 )
874            ENDIF
875            SELECT CASE( igrd )
876               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
877               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
878               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
879            END SELECT
880            DO ib = ibeg, iend
881               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
882               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
883               idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij)
884            END DO
885         END DO
886
887         ! Calculate relationship of V direction to the local orientation of the boundary
888         ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
889         ! flagv =  0 : v is tangential
890         ! flagv =  1 : v is normal to the boundary and is direction is inward
891         DO igrd = 1, jpbgrd 
892            SELECT CASE( igrd )
893               CASE( 1 )   ;   zmask => pvmask     ;   j_offset = 0
894               CASE( 2 )   ;   zmask => pfmask     ;   j_offset = 0
895               CASE( 3 )   ;   zmask => bdytmask   ;   j_offset = 1
896            END SELECT
897            icount = 0
898            ztmp(:,:) = -999._wp
899            IF( lrim0 ) THEN   ! extent of rim 0
900               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
901            ELSE               ! extent of rim 1
902               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
903            END IF
904            DO ib = ibeg, iend
905               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
906               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
907               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE
908               zsfl = zmask(ii,ij+j_offset-1)
909               znfl = zmask(ii,ij+j_offset  )
910               ! This error check only works if you are using the bdyXmask arrays
911               IF( j_offset == 1 .and. znfl + zsfl == 2. ) THEN
912                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
913                  icount = icount + 1
914               ELSE
915                  ztmp(ii,ij) = -zsfl + znfl
916               END IF
917            END DO
918            IF( icount /= 0 ) THEN
919               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
920                  ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
921               CALL ctl_stop( ctmp1 )
922            ENDIF
923            SELECT CASE( igrd )
924               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
925               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
926               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
927            END SELECT
928            DO ib = ibeg, iend
929               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
930               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
931               idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij)
932            END DO
933         END DO
934         !
935      END DO ! ib_bdy
936     
937      DO ib_bdy = 1, nb_bdy
938         DO igrd = 1, jpbgrd
939            SELECT CASE( igrd )
940               CASE( 1 )   ;   zmask => bdytmask 
941               CASE( 2 )   ;   zmask => bdyumask 
942               CASE( 3 )   ;   zmask => bdyvmask 
943            END SELECT
944            ztmp(:,:) = -999._wp
945            IF( lrim0 ) THEN   ! extent of rim 0
946               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
947            ELSE               ! extent of rim 1
948               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
949            END IF
950            DO ib = ibeg, iend
951               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
952               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
953               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE
954               llnon = zmask(ii  ,ij+1) == 1. 
955               llson = zmask(ii  ,ij-1) == 1. 
956               llean = zmask(ii+1,ij  ) == 1. 
957               llwen = zmask(ii-1,ij  ) == 1. 
958               inn  = COUNT( (/ llnon, llson, llean, llwen /) )
959               IF( inn == 0 ) THEN   ! no neighbours -> interior of a corner  or  cluster of rim points
960                  !               !              !     _____     !     _____    !    __     __
961                  !  1 |   o      !  2  o   |    !  3 | x        !  4     x |   !      |   |   -> error
962                  !    |_x_ _     !    _ _x_|    !    |   o      !      o   |   !      |x_x|
963                  IF(     zmask(ii+1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 1.
964                  ELSEIF( zmask(ii-1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 2.
965                  ELSEIF( zmask(ii+1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 3.
966                  ELSEIF( zmask(ii-1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 4.
967                  ELSE                                    ;   ztmp(ii,ij) = -1.
968                     WRITE(ctmp1,*) 'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
969                       ' on boundary set ', ib_bdy, ' has no free ocean neighbour'
970                     IF( lrim0 ) THEN
971                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 0 points.'
972                     ELSE
973                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.'
974                     END IF
975                     CALL ctl_warn( ctmp1, ctmp2 )
976                  END IF
977               END IF
978               IF( inn == 1 ) THEN   ! middle of linear bdy  or incomplete corner  ! ___ o
979                  !    |         !         |   !      o     !    ______            !    |x___
980                  ! 5  | x o     ! 6   o x |   ! 7  __x__   ! 8    x
981                  !    |         !         |   !            !      o
982                  IF( llean )   ztmp(ii,ij) = 5.
983                  IF( llwen )   ztmp(ii,ij) = 6.
984                  IF( llnon )   ztmp(ii,ij) = 7.
985                  IF( llson )   ztmp(ii,ij) = 8.
986               END IF
987               IF( inn == 2 ) THEN   ! exterior of a corner
988                  !        o      !        o      !    _____|       !       |_____ 
989                  !  9 ____x o    ! 10   o x___   ! 11     x o      ! 12   o x     
990                  !         |     !       |       !        o        !        o
991                  IF( llnon .AND. llean )   ztmp(ii,ij) =  9.
992                  IF( llnon .AND. llwen )   ztmp(ii,ij) = 10.
993                  IF( llson .AND. llean )   ztmp(ii,ij) = 11.
994                  IF( llson .AND. llwen )   ztmp(ii,ij) = 12.
995               END IF
996               IF( inn == 3 ) THEN   ! 3 neighbours     __   __
997                  !    |_  o      !        o  _|  !       |_|     !       o         
998                  ! 13  _| x o    ! 14   o x |_   ! 15   o x o    ! 16  o x o       
999                  !    |   o      !        o   |  !        o      !    __|¨|__   
1000                  IF( llnon .AND. llean .AND. llson )   ztmp(ii,ij) = 13.
1001                  IF( llnon .AND. llwen .AND. llson )   ztmp(ii,ij) = 14.
1002                  IF( llwen .AND. llson .AND. llean )   ztmp(ii,ij) = 15.
1003                  IF( llwen .AND. llnon .AND. llean )   ztmp(ii,ij) = 16.
1004               END IF
1005               IF( inn == 4 ) THEN
1006                  WRITE(ctmp1,*)  'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
1007                       ' on boundary set ', ib_bdy, ' have 4 neighbours'
1008                  CALL ctl_stop( ctmp1 )
1009               END IF
1010            END DO
1011            SELECT CASE( igrd )
1012               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
1013               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
1014               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
1015            END SELECT
1016            DO ib = ibeg, iend
1017               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1018               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1019               idx_bdy(ib_bdy)%ntreat(ib,igrd) = NINT(ztmp(ii,ij))
1020            END DO
1021         END DO
1022      END DO
1023
1024    END SUBROUTINE bdy_rim_treat
1025
1026   
1027    SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )
1028      !!----------------------------------------------------------------------
1029      !!                 ***  ROUTINE find_neib  ***
1030      !!
1031      !! ** Purpose :   get ii1, ij1, ii2, ij2, ii3, ij3, the indices of
1032      !!               the free ocean neighbours of (ii,ij) for bdy treatment
1033      !!
1034      !! ** Method  :  use itreat input to select a case
1035      !!               N.B. ntreat is defined for all bdy points in routine bdy_rim_treat
1036      !!
1037      !!----------------------------------------------------------------------
1038      INTEGER, INTENT(in   )      ::   ii, ij, itreat
1039      INTEGER, INTENT(  out)      ::   ii1, ij1, ii2, ij2, ii3, ij3
1040      !!----------------------------------------------------------------------
1041      SELECT CASE( itreat )   ! points that will be used by bdy routines, -1 will be discarded
1042         !               !               !     _____     !     _____     
1043         !  1 |   o      !  2  o   |     !  3 | x        !  4     x |   
1044         !    |_x_ _     !    _ _x_|     !    |   o      !      o   |
1045      CASE( 1 )    ;   ii1 = ii+1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1046      CASE( 2 )    ;   ii1 = ii-1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1047      CASE( 3 )    ;   ii1 = ii+1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1048      CASE( 4 )    ;   ii1 = ii-1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1049         !    |          !         |     !      o        !    ______                   ! or incomplete corner
1050         ! 5  | x o      ! 6   o x |     ! 7  __x__      ! 8    x                      !  7  ____ o
1051         !    |          !         |     !               !      o                      !         |x___
1052      CASE( 5 )    ;   ii1 = ii+1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1053      CASE( 6 )    ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1054      CASE( 7 )    ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1055      CASE( 8 )    ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1056         !        o      !        o      !    _____|     !       |_____ 
1057         !  9 ____x o    ! 10   o x___   ! 11     x o    ! 12   o x     
1058         !         |     !       |       !        o      !        o     
1059      CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
1060      CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1061      CASE( 11 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1062      CASE( 12 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1063         !    |_  o      !        o  _|  !     ¨¨|_|¨¨   !       o         
1064         ! 13  _| x o    !  14  o x |_   !  15  o x o    ! 16  o x o       
1065         !    |   o      !        o   |  !        o      !    __|¨|__
1066      CASE( 13 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1   
1067      CASE( 14 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1 
1068      CASE( 15 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij   
1069      CASE( 16 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij
1070      END SELECT
1071   END SUBROUTINE find_neib
1072   
1073
1074   SUBROUTINE bdy_read_seg( kb_bdy, knblendta ) 
1075      !!----------------------------------------------------------------------
1076      !!                 ***  ROUTINE bdy_coords_seg  ***
1077      !!
1078      !! ** Purpose :  build bdy coordinates with segments defined in namelist
1079      !!
1080      !! ** Method  :  read namelist nambdy_index blocks
1081      !!
1082      !!----------------------------------------------------------------------
1083      INTEGER                   , INTENT (in   ) ::   kb_bdy           ! bdy number
1084      INTEGER, DIMENSION(jpbgrd), INTENT (  out) ::   knblendta        ! length of index arrays
1085      !!
1086      INTEGER          ::   ios                 ! Local integer output status for namelist read
1087      INTEGER          ::   nbdyind, nbdybeg, nbdyend
1088      CHARACTER(LEN=1) ::   ctypebdy   !     -        -
1089      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
1090      !!----------------------------------------------------------------------
1091
1092      ! No REWIND here because may need to read more than one nambdy_index namelist.
1093      ! Read only namelist_cfg to avoid unseccessfull overwrite
1094      ! keep full control of the configuration namelist
1095      READ  ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 )
1096904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp )
1097      IF(lwm) WRITE ( numond, nambdy_index )
1098     
1099      SELECT CASE ( TRIM(ctypebdy) )
1100      CASE( 'N' )
1101         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1102            nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain.
1103            nbdybeg  = 2
1104            nbdyend  = jpiglo - 1
1105         ENDIF
1106         nbdysegn = nbdysegn + 1
1107         npckgn(nbdysegn) = kb_bdy ! Save bdy package number
1108         jpjnob(nbdysegn) = nbdyind
1109         jpindt(nbdysegn) = nbdybeg
1110         jpinft(nbdysegn) = nbdyend
1111         !
1112      CASE( 'S' )
1113         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1114            nbdyind  = 2           ! set boundary to whole side of model domain.
1115            nbdybeg  = 2
1116            nbdyend  = jpiglo - 1
1117         ENDIF
1118         nbdysegs = nbdysegs + 1
1119         npckgs(nbdysegs) = kb_bdy ! Save bdy package number
1120         jpjsob(nbdysegs) = nbdyind
1121         jpisdt(nbdysegs) = nbdybeg
1122         jpisft(nbdysegs) = nbdyend
1123         !
1124      CASE( 'E' )
1125         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1126            nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain.
1127            nbdybeg  = 2
1128            nbdyend  = jpjglo - 1
1129         ENDIF
1130         nbdysege = nbdysege + 1 
1131         npckge(nbdysege) = kb_bdy ! Save bdy package number
1132         jpieob(nbdysege) = nbdyind
1133         jpjedt(nbdysege) = nbdybeg
1134         jpjeft(nbdysege) = nbdyend
1135         !
1136      CASE( 'W' )
1137         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1138            nbdyind  = 2           ! set boundary to whole side of model domain.
1139            nbdybeg  = 2
1140            nbdyend  = jpjglo - 1
1141         ENDIF
1142         nbdysegw = nbdysegw + 1
1143         npckgw(nbdysegw) = kb_bdy ! Save bdy package number
1144         jpiwob(nbdysegw) = nbdyind
1145         jpjwdt(nbdysegw) = nbdybeg
1146         jpjwft(nbdysegw) = nbdyend
1147         !
1148      CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
1149      END SELECT
1150     
1151      ! For simplicity we assume that in case of straight bdy, arrays have the same length
1152      ! (even if it is true that last tangential velocity points
1153      ! are useless). This simplifies a little bit boundary data format (and agrees with format
1154      ! used so far in obc package)
1155     
1156      knblendta(1:jpbgrd) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy)
1157     
1158   END SUBROUTINE bdy_read_seg
1159
1160   
1161   SUBROUTINE bdy_ctl_seg
1162      !!----------------------------------------------------------------------
1163      !!                 ***  ROUTINE bdy_ctl_seg  ***
1164      !!
1165      !! ** Purpose :   Check straight open boundary segments location
1166      !!
1167      !! ** Method  :   - Look for open boundary corners
1168      !!                - Check that segments start or end on land
1169      !!----------------------------------------------------------------------
1170      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest 
1171      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns 
1172      REAL(wp), DIMENSION(2) ::   ztestmask
1173      !!----------------------------------------------------------------------
1174      !
1175      IF (lwp) WRITE(numout,*) ' '
1176      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1177      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1178      !
1179      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege
1180      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw
1181      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn
1182      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs
1183      ! 1. Check bounds
1184      !----------------
1185      DO ib = 1, nbdysegn
1186         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
1187         IF ((jpjnob(ib).ge.jpjglo-1).or.& 
1188            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1189         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1190         IF (jpindt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1191         IF (jpinft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
1192      END DO
1193      !
1194      DO ib = 1, nbdysegs
1195         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
1196         IF ((jpjsob(ib).ge.jpjglo-1).or.& 
1197            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1198         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1199         IF (jpisdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1200         IF (jpisft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
1201      END DO
1202      !
1203      DO ib = 1, nbdysege
1204         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib)
1205         IF ((jpieob(ib).ge.jpiglo-1).or.& 
1206            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1207         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1208         IF (jpjedt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1209         IF (jpjeft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
1210      END DO
1211      !
1212      DO ib = 1, nbdysegw
1213         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib)
1214         IF ((jpiwob(ib).ge.jpiglo-1).or.& 
1215            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1216         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1217         IF (jpjwdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1218         IF (jpjwft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
1219      ENDDO
1220      !
1221      !     
1222      ! 2. Look for segment crossings
1223      !------------------------------
1224      IF (lwp) WRITE(numout,*) '**Look for segments corners  :'
1225      !
1226      itest = 0 ! corner number
1227      !
1228      ! flag to detect if start or end of open boundary belongs to a corner
1229      ! if not (=0), it must be on land.
1230      ! if a corner is detected, save bdy package number for further tests
1231      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1232      ! South/West crossings
1233      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1234         DO ib1 = 1, nbdysegw       
1235            DO ib2 = 1, nbdysegs
1236               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1237                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1238                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1239                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1240                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN 
1241                     ! We have a possible South-West corner                     
1242!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1243!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1244                     icornw(ib1,1) = npckgs(ib2)
1245                     icorns(ib2,1) = npckgw(ib1)
1246                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1247                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1248                        &                                     jpisft(ib2), jpjwft(ib1)
1249                     WRITE(ctmp2,*) ' Not allowed yet'
1250                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 
1251                        &                            ' and South segment: ',npckgs(ib2)
1252                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1253                  ELSE
1254                     WRITE(ctmp1,*) ' Check South and West Open boundary indices'
1255                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , &
1256                        &                            ' and South segment: ',npckgs(ib2)
1257                     CALL ctl_stop( ctmp1, ctmp2 )
1258                  END IF
1259               END IF
1260            END DO
1261         END DO
1262      END IF
1263      !
1264      ! South/East crossings
1265      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1266         DO ib1 = 1, nbdysege
1267            DO ib2 = 1, nbdysegs
1268               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1269                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1270                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. &
1271                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN
1272                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1273                     ! We have a possible South-East corner
1274!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1275!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1276                     icorne(ib1,1) = npckgs(ib2)
1277                     icorns(ib2,2) = npckge(ib1)
1278                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1279                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1280                        &                                     jpisdt(ib2), jpjeft(ib1)
1281                     WRITE(ctmp2,*) ' Not allowed yet'
1282                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1283                        &                            ' and South segment: ',npckgs(ib2)
1284                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1285                  ELSE
1286                     WRITE(ctmp1,*) ' Check South and East Open boundary indices'
1287                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1288                     &                               ' and South segment: ',npckgs(ib2)
1289                     CALL ctl_stop( ctmp1, ctmp2 )
1290                  END IF
1291               END IF
1292            END DO
1293         END DO
1294      END IF
1295      !
1296      ! North/West crossings
1297      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1298         DO ib1 = 1, nbdysegw       
1299            DO ib2 = 1, nbdysegn
1300               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. &
1301                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. &
1302                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1303                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1304                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1305                     ! We have a possible North-West corner
1306!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1307!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1308                     icornw(ib1,2) = npckgn(ib2)
1309                     icornn(ib2,1) = npckgw(ib1)
1310                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1311                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1312                     &                                     jpinft(ib2), jpjwdt(ib1)
1313                     WRITE(ctmp2,*) ' Not allowed yet'
1314                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1315                     &                               ' and North segment: ',npckgn(ib2)
1316                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1317                  ELSE
1318                     WRITE(ctmp1,*) ' Check North and West Open boundary indices'
1319                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1320                     &                               ' and North segment: ',npckgn(ib2)
1321                     CALL ctl_stop( ctmp1, ctmp2 )
1322                  END IF
1323               END IF
1324            END DO
1325         END DO
1326      END IF
1327      !
1328      ! North/East crossings
1329      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1330         DO ib1 = 1, nbdysege       
1331            DO ib2 = 1, nbdysegn
1332               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1333                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1334                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1335                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1336                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1337                     ! We have a possible North-East corner
1338!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1339!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1340                     icorne(ib1,2) = npckgn(ib2)
1341                     icornn(ib2,2) = npckge(ib1)
1342                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1343                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1344                     &                                     jpindt(ib2), jpjedt(ib1)
1345                     WRITE(ctmp2,*) ' Not allowed yet'
1346                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1347                     &                               ' and North segment: ',npckgn(ib2)
1348                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1349                  ELSE
1350                     WRITE(ctmp1,*) ' Check North and East Open boundary indices'
1351                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1352                     &                               ' and North segment: ',npckgn(ib2)
1353                     CALL ctl_stop( ctmp1, ctmp2 )
1354                  END IF
1355               END IF
1356            END DO
1357         END DO
1358      END IF
1359      !
1360      ! 3. Check if segment extremities are on land
1361      !--------------------------------------------
1362      !
1363      ! West segments
1364      DO ib = 1, nbdysegw
1365         ! get mask at boundary extremities:
1366         ztestmask(1:2)=0.
1367         DO ji = 1, jpi
1368            DO jj = 1, jpj             
1369              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1370               &  ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1371              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1372               &  ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1373            END DO
1374         END DO
1375         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1376
1377         IF (ztestmask(1)==1) THEN
1378            IF (icornw(ib,1)==0) THEN
1379               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1380               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1381            ELSE
1382               ! This is a corner
1383               IF(lwp) WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1384               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1385               itest=itest+1
1386            ENDIF
1387         ENDIF
1388         IF (ztestmask(2)==1) THEN
1389            IF (icornw(ib,2)==0) THEN
1390               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1391               CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' )
1392            ELSE
1393               ! This is a corner
1394               IF(lwp) WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1395               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1396               itest=itest+1
1397            ENDIF
1398         ENDIF
1399      END DO
1400      !
1401      ! East segments
1402      DO ib = 1, nbdysege
1403         ! get mask at boundary extremities:
1404         ztestmask(1:2)=0.
1405         DO ji = 1, jpi
1406            DO jj = 1, jpj             
1407              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1408               &  ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
1409              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1410               &  ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1411            END DO
1412         END DO
1413         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1414
1415         IF (ztestmask(1)==1) THEN
1416            IF (icorne(ib,1)==0) THEN
1417               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1418               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1419            ELSE
1420               ! This is a corner
1421               IF(lwp) WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1422               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1423               itest=itest+1
1424            ENDIF
1425         ENDIF
1426         IF (ztestmask(2)==1) THEN
1427            IF (icorne(ib,2)==0) THEN
1428               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1429               CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1430            ELSE
1431               ! This is a corner
1432               IF(lwp) WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1433               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1434               itest=itest+1
1435            ENDIF
1436         ENDIF
1437      END DO
1438      !
1439      ! South segments
1440      DO ib = 1, nbdysegs
1441         ! get mask at boundary extremities:
1442         ztestmask(1:2)=0.
1443         DO ji = 1, jpi
1444            DO jj = 1, jpj             
1445              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1446               &  ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1447              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1448               &  ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1449            END DO
1450         END DO
1451         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1452
1453         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1454            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1455            CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1456         ENDIF
1457         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1458            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1459            CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1460         ENDIF
1461      END DO
1462      !
1463      ! North segments
1464      DO ib = 1, nbdysegn
1465         ! get mask at boundary extremities:
1466         ztestmask(1:2)=0.
1467         DO ji = 1, jpi
1468            DO jj = 1, jpj             
1469              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1470               &  ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
1471              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1472               &  ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1473            END DO
1474         END DO
1475         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1476
1477         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1478            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1479            CALL ctl_stop( ctmp1, ' does not start on land' )
1480         ENDIF
1481         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1482            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1483            CALL ctl_stop( ctmp1, ' does not end on land' )
1484         ENDIF
1485      END DO
1486      !
1487      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1488      !
1489      ! Other tests TBD:
1490      ! segments completly on land
1491      ! optimized open boundary array length according to landmask
1492      ! Nudging layers that overlap with interior domain
1493      !
1494   END SUBROUTINE bdy_ctl_seg
1495
1496   
1497   SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta ) 
1498      !!----------------------------------------------------------------------
1499      !!                 ***  ROUTINE bdy_coords_seg  ***
1500      !!
1501      !! ** Purpose :  build nbidta, nbidta, nbrdta for bdy built with segments
1502      !!
1503      !! ** Method  : 
1504      !!
1505      !!----------------------------------------------------------------------
1506      INTEGER, DIMENSION(:,:,:), intent(  out)  ::   nbidta, nbjdta, nbrdta   ! Index arrays: i and j indices of bdy dta
1507      !!
1508      INTEGER  ::   ii, ij, ir, iseg
1509      INTEGER  ::   igrd         ! grid type (t=1, u=2, v=3)
1510      INTEGER  ::   icount       !
1511      INTEGER  ::   ib_bdy       ! bdy number
1512      !!----------------------------------------------------------------------
1513
1514      ! East
1515      !-----
1516      DO iseg = 1, nbdysege
1517         ib_bdy = npckge(iseg)
1518         !
1519         ! ------------ T points -------------
1520         igrd=1
1521         icount=0
1522         DO ir = 1, nn_rimwidth(ib_bdy)
1523            DO ij = jpjedt(iseg), jpjeft(iseg)
1524               icount = icount + 1
1525               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
1526               nbjdta(icount, igrd, ib_bdy) = ij
1527               nbrdta(icount, igrd, ib_bdy) = ir
1528            ENDDO
1529         ENDDO
1530         !
1531         ! ------------ U points -------------
1532         igrd=2
1533         icount=0
1534         DO ir = 1, nn_rimwidth(ib_bdy)
1535            DO ij = jpjedt(iseg), jpjeft(iseg)
1536               icount = icount + 1
1537               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir
1538               nbjdta(icount, igrd, ib_bdy) = ij
1539               nbrdta(icount, igrd, ib_bdy) = ir
1540            ENDDO
1541         ENDDO
1542         !
1543         ! ------------ V points -------------
1544         igrd=3
1545         icount=0
1546         DO ir = 1, nn_rimwidth(ib_bdy)
1547            !            DO ij = jpjedt(iseg), jpjeft(iseg) - 1
1548            DO ij = jpjedt(iseg), jpjeft(iseg)
1549               icount = icount + 1
1550               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
1551               nbjdta(icount, igrd, ib_bdy) = ij
1552               nbrdta(icount, igrd, ib_bdy) = ir
1553            ENDDO
1554            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1555            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1556         ENDDO
1557      ENDDO
1558      !
1559      ! West
1560      !-----
1561      DO iseg = 1, nbdysegw
1562         ib_bdy = npckgw(iseg)
1563         !
1564         ! ------------ T points -------------
1565         igrd=1
1566         icount=0
1567         DO ir = 1, nn_rimwidth(ib_bdy)
1568            DO ij = jpjwdt(iseg), jpjwft(iseg)
1569               icount = icount + 1
1570               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
1571               nbjdta(icount, igrd, ib_bdy) = ij
1572               nbrdta(icount, igrd, ib_bdy) = ir
1573            ENDDO
1574         ENDDO
1575         !
1576         ! ------------ U points -------------
1577         igrd=2
1578         icount=0
1579         DO ir = 1, nn_rimwidth(ib_bdy)
1580            DO ij = jpjwdt(iseg), jpjwft(iseg)
1581               icount = icount + 1
1582               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
1583               nbjdta(icount, igrd, ib_bdy) = ij
1584               nbrdta(icount, igrd, ib_bdy) = ir
1585            ENDDO
1586         ENDDO
1587         !
1588         ! ------------ V points -------------
1589         igrd=3
1590         icount=0
1591         DO ir = 1, nn_rimwidth(ib_bdy)
1592            !            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
1593            DO ij = jpjwdt(iseg), jpjwft(iseg)
1594               icount = icount + 1
1595               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
1596               nbjdta(icount, igrd, ib_bdy) = ij
1597               nbrdta(icount, igrd, ib_bdy) = ir
1598            ENDDO
1599            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1600            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1601         ENDDO
1602      ENDDO
1603      !
1604      ! North
1605      !-----
1606      DO iseg = 1, nbdysegn
1607         ib_bdy = npckgn(iseg)
1608         !
1609         ! ------------ T points -------------
1610         igrd=1
1611         icount=0
1612         DO ir = 1, nn_rimwidth(ib_bdy)
1613            DO ii = jpindt(iseg), jpinft(iseg)
1614               icount = icount + 1
1615               nbidta(icount, igrd, ib_bdy) = ii
1616               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
1617               nbrdta(icount, igrd, ib_bdy) = ir
1618            ENDDO
1619         ENDDO
1620         !
1621         ! ------------ U points -------------
1622         igrd=2
1623         icount=0
1624         DO ir = 1, nn_rimwidth(ib_bdy)
1625            !            DO ii = jpindt(iseg), jpinft(iseg) - 1
1626            DO ii = jpindt(iseg), jpinft(iseg)
1627               icount = icount + 1
1628               nbidta(icount, igrd, ib_bdy) = ii
1629               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
1630               nbrdta(icount, igrd, ib_bdy) = ir
1631            ENDDO
1632            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1633            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1634         ENDDO
1635         !
1636         ! ------------ V points -------------
1637         igrd=3
1638         icount=0
1639         DO ir = 1, nn_rimwidth(ib_bdy)
1640            DO ii = jpindt(iseg), jpinft(iseg)
1641               icount = icount + 1
1642               nbidta(icount, igrd, ib_bdy) = ii
1643               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir
1644               nbrdta(icount, igrd, ib_bdy) = ir
1645            ENDDO
1646         ENDDO
1647      ENDDO
1648      !
1649      ! South
1650      !-----
1651      DO iseg = 1, nbdysegs
1652         ib_bdy = npckgs(iseg)
1653         !
1654         ! ------------ T points -------------
1655         igrd=1
1656         icount=0
1657         DO ir = 1, nn_rimwidth(ib_bdy)
1658            DO ii = jpisdt(iseg), jpisft(iseg)
1659               icount = icount + 1
1660               nbidta(icount, igrd, ib_bdy) = ii
1661               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
1662               nbrdta(icount, igrd, ib_bdy) = ir
1663            ENDDO
1664         ENDDO
1665         !
1666         ! ------------ U points -------------
1667         igrd=2
1668         icount=0
1669         DO ir = 1, nn_rimwidth(ib_bdy)
1670            !            DO ii = jpisdt(iseg), jpisft(iseg) - 1
1671            DO ii = jpisdt(iseg), jpisft(iseg)
1672               icount = icount + 1
1673               nbidta(icount, igrd, ib_bdy) = ii
1674               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
1675               nbrdta(icount, igrd, ib_bdy) = ir
1676            ENDDO
1677            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1678            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1679         ENDDO
1680         !
1681         ! ------------ V points -------------
1682         igrd=3
1683         icount=0
1684         DO ir = 1, nn_rimwidth(ib_bdy)
1685            DO ii = jpisdt(iseg), jpisft(iseg)
1686               icount = icount + 1
1687               nbidta(icount, igrd, ib_bdy) = ii
1688               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
1689               nbrdta(icount, igrd, ib_bdy) = ir
1690            ENDDO
1691         ENDDO
1692      ENDDO
1693
1694     
1695   END SUBROUTINE bdy_coords_seg
1696   
1697   
1698   SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1699      !!----------------------------------------------------------------------
1700      !!                 ***  ROUTINE bdy_ctl_corn  ***
1701      !!
1702      !! ** Purpose :   Check numerical schemes consistency between
1703      !!                segments having a common corner
1704      !!
1705      !! ** Method  :   
1706      !!----------------------------------------------------------------------
1707      INTEGER, INTENT(in)  ::   ib1, ib2
1708      INTEGER :: itest
1709      !!----------------------------------------------------------------------
1710      itest = 0
1711
1712      IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) )   itest = itest + 1
1713      IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) )   itest = itest + 1
1714      IF( cn_tra  (ib1) /= cn_tra  (ib2) )   itest = itest + 1
1715      !
1716      IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) )   itest = itest + 1
1717      IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) )   itest = itest + 1
1718      IF( nn_tra_dta  (ib1) /= nn_tra_dta  (ib2) )   itest = itest + 1
1719      !
1720      IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) )   itest = itest + 1   
1721      !
1722      IF( itest>0 ) THEN
1723         WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2
1724         CALL ctl_stop( ctmp1, ' have different open bdy schemes' )
1725      ENDIF
1726      !
1727   END SUBROUTINE bdy_ctl_corn
1728
1729
1730   SUBROUTINE bdy_meshwri()
1731      !!----------------------------------------------------------------------
1732      !!                 ***  ROUTINE bdy_meshwri  ***
1733      !!         
1734      !! ** Purpose :   write netcdf file with nbr, flagu, flagv, ntreat for T, U
1735      !!                and V points in 2D arrays for easier visualisation/control
1736      !!
1737      !! ** Method  :   use iom_rstput as in domwri.F
1738      !!----------------------------------------------------------------------     
1739      INTEGER  ::   ib_bdy, ii, ij, igrd, ib     ! dummy loop indices
1740      INTEGER  ::   inum                                   !   -       -
1741      REAL(wp), POINTER, DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields
1742      REAL(wp)         , DIMENSION(jpi,jpj) ::   ztmp
1743      CHARACTER(LEN=1) , DIMENSION(jpbgrd)  ::   cgrid
1744      !!----------------------------------------------------------------------     
1745      cgrid = (/'t','u','v'/)
1746      CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. )
1747      DO igrd = 1, jpbgrd
1748         SELECT CASE( igrd )
1749         CASE( 1 )   ;   zmask => tmask(:,:,1)
1750         CASE( 2 )   ;   zmask => umask(:,:,1)
1751         CASE( 3 )   ;   zmask => vmask(:,:,1)
1752         END SELECT
1753         ztmp(:,:) = zmask(:,:)
1754         DO ib_bdy = 1, nb_bdy
1755            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)      ! nbr deined for all rims
1756               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1757               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1758               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10.
1759               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1760            END DO
1761         END DO
1762         CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1763         ztmp(:,:) = zmask(:,:)
1764         DO ib_bdy = 1, nb_bdy
1765            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagu defined only for rims 0 and 1
1766               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1767               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1768               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10.
1769               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1770            END DO
1771         END DO
1772         CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1773         ztmp(:,:) = zmask(:,:)
1774         DO ib_bdy = 1, nb_bdy
1775            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagv defined only for rims 0 and 1
1776               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1777               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1778               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10.
1779               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1780            END DO
1781         END DO
1782         CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1783         ztmp(:,:) = zmask(:,:)
1784         DO ib_bdy = 1, nb_bdy
1785            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! ntreat defined only for rims 0 and 1
1786               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1787               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1788               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10.
1789               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1790            END DO
1791         END DO
1792         CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1793      END DO
1794      CALL iom_close( inum )
1795
1796   END SUBROUTINE bdy_meshwri
1797   
1798   !!=================================================================================
1799END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.