New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
bdyini.F90 in NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/BDY – NEMO

source: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/BDY/bdyini.F90 @ 15127

Last change on this file since 15127 was 15127, checked in by cetlod, 3 years ago

dev_PISCO : merge with trunk@15119

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