source: NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_CO8_package/src/OCE/BDY/bdyini.F90 @ 12094

Last change on this file since 12094 was 12094, checked in by deazer, 18 months ago

Allow ts_init to work properly for already interpolated scoords
Allow bdy to be spread over arbitary rim points

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