source: NEMO/trunk/src/OCE/BDY/bdyini.F90

Last change on this file was 13286, checked in by smasson, 5 weeks ago

trunk: merge extra halos branch in trunk, see #2366

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