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

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

source: branches/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 3901

Last change on this file since 3901 was 3901, checked in by clevy, 11 years ago

Configuration Setting/Step2, see ticket:#1074
  • Property svn:keywords set to Id
File size: 74.6 KB
Line 
1MODULE bdyini
2   !!======================================================================
3   !!                       ***  MODULE  bdyini  ***
4   !! Unstructured open boundaries : initialisation
5   !!======================================================================
6   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
7   !!             -   !  2007-01  (D. Storkey) Update to use IOM module
8   !!             -   !  2007-01  (D. Storkey) Tidal forcing
9   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
10   !!            3.3  !  2010-09  (E.O'Dea) updates for Shelf configurations
11   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
12   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
13   !!            3.4  !  2012     (J. Chanut) straight open boundary case update
14   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Updates for the
15   !!                             optimization of BDY communications
16   !!----------------------------------------------------------------------
17#if defined key_bdy
18   !!----------------------------------------------------------------------
19   !!   'key_bdy'                     Unstructured Open Boundary Conditions
20   !!----------------------------------------------------------------------
21   !!   bdy_init       : Initialization of unstructured open boundaries
22   !!----------------------------------------------------------------------
23   USE timing          ! Timing
24   USE oce             ! ocean dynamics and tracers variables
25   USE dom_oce         ! ocean space and time domain
26   USE bdy_oce         ! unstructured open boundary conditions
27   USE in_out_manager  ! I/O units
28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
29   USE lib_mpp         ! for mpp_sum 
30   USE iom             ! I/O
31   USE sbctide, ONLY: lk_tide ! Tidal forcing or not
32   USE phycst, ONLY: rday
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   bdy_init   ! routine called in nemo_init
38
39   INTEGER, PARAMETER          :: jp_nseg = 100
40   INTEGER, PARAMETER          :: nrimmax = 20 ! maximum rimwidth in structured
41                                               ! open boundary data files
42   ! Straight open boundary segment parameters:
43   INTEGER  :: nbdysege, nbdysegw, nbdysegn, nbdysegs 
44   INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge
45   INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw
46   INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn
47   INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
50   !! $Id$
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54   
55   SUBROUTINE bdy_init
56      !!----------------------------------------------------------------------
57      !!                 ***  ROUTINE bdy_init  ***
58      !!         
59      !! ** Purpose :   Initialization of the dynamics and tracer fields with
60      !!              unstructured open boundaries.
61      !!
62      !! ** Method  :   Read initialization arrays (mask, indices) to identify
63      !!              an unstructured open boundary
64      !!
65      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
66      !!----------------------------------------------------------------------     
67      ! namelist variables
68      !-------------------
69      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile
70      CHARACTER(LEN=1)   ::   ctypebdy
71      INTEGER :: nbdyind, nbdybeg, nbdyend
72
73      ! local variables
74      !-------------------
75      INTEGER  ::   ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
76      INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers
77      INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       -
78      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       -
79      INTEGER  ::   jpbdtau, jpbdtas                       !   -       -
80      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       -
81      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts
82      REAL   , POINTER  ::  flagu, flagv                   !    -   -
83      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
84      INTEGER, DIMENSION (2)                  ::   kdimsz
85      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays
86      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
87      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
88      CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid
89      INTEGER :: com_east, com_west, com_south, com_north          ! Flags for boundaries sending
90      INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b  ! Flags for boundaries receiving
91      INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4)                ! Arrays for neighbours coordinates
92
93      !!
94      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             &
95         &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, &
96         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         & 
97         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp,              &
98#if defined key_lim2
99         &             nn_ice_lim2, nn_ice_lim2_dta,                       &
100#endif
101         &             ln_vol, nn_volctl, nn_rimwidth
102      !!
103      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
104      INTEGER  ::   ios                 ! Local integer output status for namelist read
105      !!----------------------------------------------------------------------
106
107      IF( nn_timing == 1 ) CALL timing_start('bdy_init')
108
109      IF( bdy_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' )
110
111      IF(lwp) WRITE(numout,*)
112      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
113      IF(lwp) WRITE(numout,*) '~~~~~~~~'
114      !
115
116      IF( jperio /= 0 )   CALL ctl_stop( 'Cyclic or symmetric,',   &
117         &                               ' and general open boundary condition are not compatible' )
118
119      cgrid= (/'t','u','v'/)
120     
121      ! ------------------------
122      ! Read namelist parameters
123      ! ------------------------
124
125      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries 
126      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901)
127901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
128
129      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
130      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
131902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
132      WRITE ( numond, nambdy )
133
134      ! -----------------------------------------
135      ! Check and write out namelist parameters
136      ! -----------------------------------------
137      !                                   ! control prints
138      IF(lwp) WRITE(numout,*) '         nambdy'
139
140      IF( nb_bdy .eq. 0 ) THEN
141        IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
142      ELSE
143        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy
144      ENDIF
145
146      DO ib_bdy = 1,nb_bdy
147        IF(lwp) WRITE(numout,*) ' ' 
148        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' 
149
150        IF( ln_coords_file(ib_bdy) ) THEN
151           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
152        ELSE
153           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
154        ENDIF
155        IF(lwp) WRITE(numout,*)
156
157        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
158        SELECT CASE( nn_dyn2d(ib_bdy) )                 
159          CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
160          CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
161          CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
162          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
163        END SELECT
164        IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN
165           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !
166              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
167              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
168              CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file'
169              CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files'
170              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
171           END SELECT
172           IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.lk_tide)) THEN
173             CALL ctl_stop( 'You must activate key_tide to add tidal forcing at open boundaries' )
174           ENDIF
175        ENDIF
176        IF(lwp) WRITE(numout,*)
177
178        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
179        SELECT CASE( nn_dyn3d(ib_bdy) )                 
180          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
181          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
182          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
183          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
184          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
185        END SELECT
186        IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN
187           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !
188              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
189              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
190              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
191           END SELECT
192        ENDIF
193
194        IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
195           IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN
196              IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
197              ln_dyn3d_dmp(ib_bdy)=.false.
198           ELSEIF ( nn_dyn3d(ib_bdy).EQ.1 ) THEN
199              CALL ctl_stop( 'Use FRS OR relaxation' )
200           ELSE
201              IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone'
202              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
203              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
204           ENDIF
205        ELSE
206           IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities'
207        ENDIF
208        IF(lwp) WRITE(numout,*)
209
210        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
211        SELECT CASE( nn_tra(ib_bdy) )                 
212          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
213          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
214          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
215          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions'
216          CASE( 4 )      ;   IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity'
217          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
218        END SELECT
219        IF( nn_tra(ib_bdy) .gt. 0 ) THEN
220           SELECT CASE( nn_tra_dta(ib_bdy) )                   !
221              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
222              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
223              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
224           END SELECT
225        ENDIF
226
227        IF ( ln_tra_dmp(ib_bdy) ) THEN
228           IF ( nn_tra(ib_bdy).EQ.0 ) THEN
229              IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
230              ln_tra_dmp(ib_bdy)=.false.
231           ELSEIF ( nn_tra(ib_bdy).EQ.1 ) THEN
232              CALL ctl_stop( 'Use FRS OR relaxation' )
233           ELSE
234              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
235              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
236              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
237           ENDIF
238        ELSE
239           IF(lwp) WRITE(numout,*) '      NO T/S relaxation'
240        ENDIF
241        IF(lwp) WRITE(numout,*)
242
243#if defined key_lim2
244        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
245        SELECT CASE( nn_ice_lim2(ib_bdy) )                 
246          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
247          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
248          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
249        END SELECT
250        IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN
251           SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   !
252              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
253              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
254              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' )
255           END SELECT
256        ENDIF
257        IF(lwp) WRITE(numout,*)
258#endif
259
260        IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
261        IF(lwp) WRITE(numout,*)
262
263      ENDDO
264
265     IF (nb_bdy .gt. 0) THEN
266        IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
267          IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
268          IF(lwp) WRITE(numout,*)
269          SELECT CASE ( nn_volctl )
270            CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
271            CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
272            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
273          END SELECT
274          IF(lwp) WRITE(numout,*)
275        ELSE
276          IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
277          IF(lwp) WRITE(numout,*)
278        ENDIF
279     ENDIF
280
281      ! -------------------------------------------------
282      ! Initialise indices arrays for open boundaries
283      ! -------------------------------------------------
284
285      ! Work out global dimensions of boundary data
286      ! ---------------------------------------------
287      REWIND( numnam_cfg )     
288
289      !!----------------------------------------------------------------------
290
291             
292               
293      nblendta(:,:) = 0
294      nbdysege = 0
295      nbdysegw = 0
296      nbdysegn = 0
297      nbdysegs = 0
298      icount   = 0 ! count user defined segments
299      ! Dimensions below are used to allocate arrays to read external data
300      jpbdtas = 1 ! Maximum size of boundary data (structured case)
301      jpbdtau = 1 ! Maximum size of boundary data (unstructured case)
302
303      DO ib_bdy = 1, nb_bdy
304
305         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
306 
307            icount = icount + 1
308            ! No REWIND here because may need to read more than one nambdy_index namelist.
309            ! Read only namelist_cfg to avoid unseccessfull overwrite
310!!          REWIND( numnam_ref )              ! Namelist nambdy_index in reference namelist : Open boundaries indexes
311!!          READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 903)
312!!903       IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in reference namelist', lwp )
313
314!!          REWIND( numnam_cfg )              ! Namelist nambdy_index in configuration namelist : Open boundaries indexes
315            READ  ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 )
316904         IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp )
317            WRITE ( numond, nambdy_index )
318
319            SELECT CASE ( TRIM(ctypebdy) )
320              CASE( 'N' )
321                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
322                    nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain.
323                    nbdybeg  = 2
324                    nbdyend  = jpiglo - 1
325                 ENDIF
326                 nbdysegn = nbdysegn + 1
327                 npckgn(nbdysegn) = ib_bdy ! Save bdy package number
328                 jpjnob(nbdysegn) = nbdyind
329                 jpindt(nbdysegn) = nbdybeg
330                 jpinft(nbdysegn) = nbdyend
331                 !
332              CASE( 'S' )
333                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
334                    nbdyind  = 2           ! set boundary to whole side of model domain.
335                    nbdybeg  = 2
336                    nbdyend  = jpiglo - 1
337                 ENDIF
338                 nbdysegs = nbdysegs + 1
339                 npckgs(nbdysegs) = ib_bdy ! Save bdy package number
340                 jpjsob(nbdysegs) = nbdyind
341                 jpisdt(nbdysegs) = nbdybeg
342                 jpisft(nbdysegs) = nbdyend
343                 !
344              CASE( 'E' )
345                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
346                    nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain.
347                    nbdybeg  = 2
348                    nbdyend  = jpjglo - 1
349                 ENDIF
350                 nbdysege = nbdysege + 1 
351                 npckge(nbdysege) = ib_bdy ! Save bdy package number
352                 jpieob(nbdysege) = nbdyind
353                 jpjedt(nbdysege) = nbdybeg
354                 jpjeft(nbdysege) = nbdyend
355                 !
356              CASE( 'W' )
357                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
358                    nbdyind  = 2           ! set boundary to whole side of model domain.
359                    nbdybeg  = 2
360                    nbdyend  = jpjglo - 1
361                 ENDIF
362                 nbdysegw = nbdysegw + 1
363                 npckgw(nbdysegw) = ib_bdy ! Save bdy package number
364                 jpiwob(nbdysegw) = nbdyind
365                 jpjwdt(nbdysegw) = nbdybeg
366                 jpjwft(nbdysegw) = nbdyend
367                 !
368              CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
369            END SELECT
370
371            ! For simplicity we assume that in case of straight bdy, arrays have the same length
372            ! (even if it is true that last tangential velocity points
373            ! are useless). This simplifies a little bit boundary data format (and agrees with format
374            ! used so far in obc package)
375
376            nblendta(1:jpbgrd,ib_bdy) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy)
377            jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1))
378            IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) &
379            & CALL ctl_stop( 'rimwidth must be lower than nrimmax' )
380
381         ELSE            ! Read size of arrays in boundary coordinates file.
382            CALL iom_open( cn_coords_file(ib_bdy), inum )
383            DO igrd = 1, jpbgrd
384               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 
385               nblendta(igrd,ib_bdy) = kdimsz(1)
386               jpbdtau = MAX(jpbdtau, kdimsz(1))
387            ENDDO
388            CALL iom_close( inum )
389
390         ENDIF
391
392      ENDDO ! ib_bdy
393
394      IF (nb_bdy>0) THEN
395         jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
396
397         ! Allocate arrays
398         !---------------
399         ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
400            &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
401
402         ALLOCATE( dta_global(jpbdtau, 1, jpk) )
403         IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) )
404         !
405      ENDIF
406
407      ! Now look for crossings in user (namelist) defined open boundary segments:
408      !--------------------------------------------------------------------------
409      IF ( icount>0 ) CALL bdy_ctl_seg
410
411      ! Calculate global boundary index arrays or read in from file
412      !------------------------------------------------------------               
413      ! 1. Read global index arrays from boundary coordinates file.
414      DO ib_bdy = 1, nb_bdy
415
416         IF( ln_coords_file(ib_bdy) ) THEN
417
418            CALL iom_open( cn_coords_file(ib_bdy), inum )
419            DO igrd = 1, jpbgrd
420               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
421               DO ii = 1,nblendta(igrd,ib_bdy)
422                  nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
423               END DO
424               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
425               DO ii = 1,nblendta(igrd,ib_bdy)
426                  nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
427               END DO
428               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
429               DO ii = 1,nblendta(igrd,ib_bdy)
430                  nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
431               END DO
432
433               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
434               IF(lwp) WRITE(numout,*)
435               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
436               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
437               IF (ibr_max < nn_rimwidth(ib_bdy))   &
438                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
439            END DO
440            CALL iom_close( inum )
441
442         ENDIF
443
444      ENDDO     
445   
446      ! 2. Now fill indices corresponding to straight open boundary arrays:
447      ! East
448      !-----
449      DO iseg = 1, nbdysege
450         ib_bdy = npckge(iseg)
451         !
452         ! ------------ T points -------------
453         igrd=1
454         icount=0
455         DO ir = 1, nn_rimwidth(ib_bdy)
456            DO ij = jpjedt(iseg), jpjeft(iseg)
457               icount = icount + 1
458               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
459               nbjdta(icount, igrd, ib_bdy) = ij
460               nbrdta(icount, igrd, ib_bdy) = ir
461            ENDDO
462         ENDDO
463         !
464         ! ------------ U points -------------
465         igrd=2
466         icount=0
467         DO ir = 1, nn_rimwidth(ib_bdy)
468            DO ij = jpjedt(iseg), jpjeft(iseg)
469               icount = icount + 1
470               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir
471               nbjdta(icount, igrd, ib_bdy) = ij
472               nbrdta(icount, igrd, ib_bdy) = ir
473            ENDDO
474         ENDDO
475         !
476         ! ------------ V points -------------
477         igrd=3
478         icount=0
479         DO ir = 1, nn_rimwidth(ib_bdy)
480!            DO ij = jpjedt(iseg), jpjeft(iseg) - 1
481            DO ij = jpjedt(iseg), jpjeft(iseg)
482               icount = icount + 1
483               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
484               nbjdta(icount, igrd, ib_bdy) = ij
485               nbrdta(icount, igrd, ib_bdy) = ir
486            ENDDO
487            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
488            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
489         ENDDO
490      ENDDO
491      !
492      ! West
493      !-----
494      DO iseg = 1, nbdysegw
495         ib_bdy = npckgw(iseg)
496         !
497         ! ------------ T points -------------
498         igrd=1
499         icount=0
500         DO ir = 1, nn_rimwidth(ib_bdy)
501            DO ij = jpjwdt(iseg), jpjwft(iseg)
502               icount = icount + 1
503               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
504               nbjdta(icount, igrd, ib_bdy) = ij
505               nbrdta(icount, igrd, ib_bdy) = ir
506            ENDDO
507         ENDDO
508         !
509         ! ------------ U points -------------
510         igrd=2
511         icount=0
512         DO ir = 1, nn_rimwidth(ib_bdy)
513            DO ij = jpjwdt(iseg), jpjwft(iseg)
514               icount = icount + 1
515               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
516               nbjdta(icount, igrd, ib_bdy) = ij
517               nbrdta(icount, igrd, ib_bdy) = ir
518            ENDDO
519         ENDDO
520         !
521         ! ------------ V points -------------
522         igrd=3
523         icount=0
524         DO ir = 1, nn_rimwidth(ib_bdy)
525!            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
526            DO ij = jpjwdt(iseg), jpjwft(iseg)
527               icount = icount + 1
528               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
529               nbjdta(icount, igrd, ib_bdy) = ij
530               nbrdta(icount, igrd, ib_bdy) = ir
531            ENDDO
532            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
533            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
534         ENDDO
535      ENDDO
536      !
537      ! North
538      !-----
539      DO iseg = 1, nbdysegn
540         ib_bdy = npckgn(iseg)
541         !
542         ! ------------ T points -------------
543         igrd=1
544         icount=0
545         DO ir = 1, nn_rimwidth(ib_bdy)
546            DO ii = jpindt(iseg), jpinft(iseg)
547               icount = icount + 1
548               nbidta(icount, igrd, ib_bdy) = ii
549               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
550               nbrdta(icount, igrd, ib_bdy) = ir
551            ENDDO
552         ENDDO
553         !
554         ! ------------ U points -------------
555         igrd=2
556         icount=0
557         DO ir = 1, nn_rimwidth(ib_bdy)
558!            DO ii = jpindt(iseg), jpinft(iseg) - 1
559            DO ii = jpindt(iseg), jpinft(iseg)
560               icount = icount + 1
561               nbidta(icount, igrd, ib_bdy) = ii
562               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
563               nbrdta(icount, igrd, ib_bdy) = ir
564            ENDDO
565            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
566            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
567         ENDDO
568         !
569         ! ------------ V points -------------
570         igrd=3
571         icount=0
572         DO ir = 1, nn_rimwidth(ib_bdy)
573            DO ii = jpindt(iseg), jpinft(iseg)
574               icount = icount + 1
575               nbidta(icount, igrd, ib_bdy) = ii
576               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir
577               nbrdta(icount, igrd, ib_bdy) = ir
578            ENDDO
579         ENDDO
580      ENDDO
581      !
582      ! South
583      !-----
584      DO iseg = 1, nbdysegs
585         ib_bdy = npckgs(iseg)
586         !
587         ! ------------ T points -------------
588         igrd=1
589         icount=0
590         DO ir = 1, nn_rimwidth(ib_bdy)
591            DO ii = jpisdt(iseg), jpisft(iseg)
592               icount = icount + 1
593               nbidta(icount, igrd, ib_bdy) = ii
594               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
595               nbrdta(icount, igrd, ib_bdy) = ir
596            ENDDO
597         ENDDO
598         !
599         ! ------------ U points -------------
600         igrd=2
601         icount=0
602         DO ir = 1, nn_rimwidth(ib_bdy)
603!            DO ii = jpisdt(iseg), jpisft(iseg) - 1
604            DO ii = jpisdt(iseg), jpisft(iseg)
605               icount = icount + 1
606               nbidta(icount, igrd, ib_bdy) = ii
607               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
608               nbrdta(icount, igrd, ib_bdy) = ir
609            ENDDO
610            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
611            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
612         ENDDO
613         !
614         ! ------------ V points -------------
615         igrd=3
616         icount=0
617         DO ir = 1, nn_rimwidth(ib_bdy)
618            DO ii = jpisdt(iseg), jpisft(iseg)
619               icount = icount + 1
620               nbidta(icount, igrd, ib_bdy) = ii
621               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
622               nbrdta(icount, igrd, ib_bdy) = ir
623            ENDDO
624         ENDDO
625      ENDDO
626
627      !  Deal with duplicated points
628      !-----------------------------
629      ! We assign negative indices to duplicated points (to remove them from bdy points to be updated)
630      ! if their distance to the bdy is greater than the other
631      ! If their distance are the same, just keep only one to avoid updating a point twice
632      DO igrd = 1, jpbgrd
633         DO ib_bdy1 = 1, nb_bdy
634            DO ib_bdy2 = 1, nb_bdy
635               IF (ib_bdy1/=ib_bdy2) THEN
636                  DO ib1 = 1, nblendta(igrd,ib_bdy1)
637                     DO ib2 = 1, nblendta(igrd,ib_bdy2)
638                        IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
639                        &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
640!                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &
641!                                                       &              nbidta(ib1, igrd, ib_bdy1),      &
642!                                                       &              nbjdta(ib2, igrd, ib_bdy2)
643                           ! keep only points with the lowest distance to boundary:
644                           IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN
645                             nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
646                             nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
647                           ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN
648                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
649                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
650                           ! Arbitrary choice if distances are the same:
651                           ELSE
652                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
653                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
654                           ENDIF
655                        END IF
656                     END DO
657                  END DO
658               ENDIF
659            END DO
660         END DO
661      END DO
662
663      ! Work out dimensions of boundary data on each processor
664      ! ------------------------------------------------------
665
666      ! Rather assume that boundary data indices are given on global domain
667      ! TO BE DISCUSSED ?
668!      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
669!      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
670!      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
671!      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1     
672      iw = mig(1) - jpizoom + 2         ! if monotasking and no zoom, iw=2
673      ie = mig(1) + nlci - jpizoom - 1  ! if monotasking and no zoom, ie=jpim1
674      is = mjg(1) - jpjzoom + 2         ! if monotasking and no zoom, is=2
675      in = mjg(1) + nlcj - jpjzoom - 1  ! if monotasking and no zoom, in=jpjm1
676
677      ALLOCATE( nbondi_bdy(nb_bdy))
678      ALLOCATE( nbondj_bdy(nb_bdy))
679      nbondi_bdy(:)=2
680      nbondj_bdy(:)=2
681      ALLOCATE( nbondi_bdy_b(nb_bdy))
682      ALLOCATE( nbondj_bdy_b(nb_bdy))
683      nbondi_bdy_b(:)=2
684      nbondj_bdy_b(:)=2
685
686      ! Work out dimensions of boundary data on each neighbour process
687      IF(nbondi .eq. 0) THEN
688         iw_b(1) = jpizoom + nimppt(nowe+1)
689         ie_b(1) = jpizoom + nimppt(nowe+1)+nlcit(nowe+1)-3
690         is_b(1) = jpjzoom + njmppt(nowe+1)
691         in_b(1) = jpjzoom + njmppt(nowe+1)+nlcjt(nowe+1)-3
692
693         iw_b(2) = jpizoom + nimppt(noea+1)
694         ie_b(2) = jpizoom + nimppt(noea+1)+nlcit(noea+1)-3
695         is_b(2) = jpjzoom + njmppt(noea+1)
696         in_b(2) = jpjzoom + njmppt(noea+1)+nlcjt(noea+1)-3
697      ELSEIF(nbondi .eq. 1) THEN
698         iw_b(1) = jpizoom + nimppt(nowe+1)
699         ie_b(1) = jpizoom + nimppt(nowe+1)+nlcit(nowe+1)-3
700         is_b(1) = jpjzoom + njmppt(nowe+1)
701         in_b(1) = jpjzoom + njmppt(nowe+1)+nlcjt(nowe+1)-3
702      ELSEIF(nbondi .eq. -1) THEN
703         iw_b(2) = jpizoom + nimppt(noea+1)
704         ie_b(2) = jpizoom + nimppt(noea+1)+nlcit(noea+1)-3
705         is_b(2) = jpjzoom + njmppt(noea+1)
706         in_b(2) = jpjzoom + njmppt(noea+1)+nlcjt(noea+1)-3
707      ENDIF
708
709      IF(nbondj .eq. 0) THEN
710         iw_b(3) = jpizoom + nimppt(noso+1)
711         ie_b(3) = jpizoom + nimppt(noso+1)+nlcit(noso+1)-3
712         is_b(3) = jpjzoom + njmppt(noso+1)
713         in_b(3) = jpjzoom + njmppt(noso+1)+nlcjt(noso+1)-3
714
715         iw_b(4) = jpizoom + nimppt(nono+1)
716         ie_b(4) = jpizoom + nimppt(nono+1)+nlcit(nono+1)-3
717         is_b(4) = jpjzoom + njmppt(nono+1)
718         in_b(4) = jpjzoom + njmppt(nono+1)+nlcjt(nono+1)-3
719      ELSEIF(nbondj .eq. 1) THEN
720         iw_b(3) = jpizoom + nimppt(noso+1)
721         ie_b(3) = jpizoom + nimppt(noso+1)+nlcit(noso+1)-3
722         is_b(3) = jpjzoom + njmppt(noso+1)
723         in_b(3) = jpjzoom + njmppt(noso+1)+nlcjt(noso+1)-3
724      ELSEIF(nbondj .eq. -1) THEN
725         iw_b(4) = jpizoom + nimppt(nono+1)
726         ie_b(4) = jpizoom + nimppt(nono+1)+nlcit(nono+1)-3
727         is_b(4) = jpjzoom + njmppt(nono+1)
728         in_b(4) = jpjzoom + njmppt(nono+1)+nlcjt(nono+1)-3
729      ENDIF
730
731      DO ib_bdy = 1, nb_bdy
732         DO igrd = 1, jpbgrd
733            icount  = 0
734            icountr = 0
735            idx_bdy(ib_bdy)%nblen(igrd)    = 0
736            idx_bdy(ib_bdy)%nblenrim(igrd) = 0
737            DO ib = 1, nblendta(igrd,ib_bdy)
738               ! check that data is in correct order in file
739               ibm1 = MAX(1,ib-1)
740               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc...
741                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
742                     CALL ctl_stop('bdy_init : ERROR : boundary data in file  &
743                                    must be defined in order of distance from edge nbr.', &
744                                   'A utility for re-ordering boundary coordinates and data &
745                                    files exists in the TOOLS/OBC directory')
746                  ENDIF   
747               ENDIF
748               ! check if point is in local domain
749               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
750                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN
751                  !
752                  icount = icount  + 1
753                  !
754                  IF( nbrdta(ib,igrd,ib_bdy) == 1 )   icountr = icountr+1
755               ENDIF
756            ENDDO
757            idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
758            idx_bdy(ib_bdy)%nblen   (igrd) = icount  !: length of boundary data on each proc       
759         ENDDO  ! igrd
760
761         ! Allocate index arrays for this boundary set
762         !--------------------------------------------
763         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
764         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
765         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
766         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
767         ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) )
768         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
769         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
770         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) )
771         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )
772
773         ! Dispatch mapping indices and discrete distances on each processor
774         ! -----------------------------------------------------------------
775
776         com_east = 0
777         com_west = 0
778         com_south = 0
779         com_north = 0
780
781         com_east_b = 0
782         com_west_b = 0
783         com_south_b = 0
784         com_north_b = 0
785         DO igrd = 1, jpbgrd
786            icount  = 0
787            ! Loop on rimwidth to ensure outermost points come first in the local arrays.
788            DO ir=1, nn_rimwidth(ib_bdy)
789               DO ib = 1, nblendta(igrd,ib_bdy)
790                  ! check if point is in local domain and equals ir
791                  IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
792                     & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   &
793                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
794                     !
795                     icount = icount  + 1
796
797                     ! Rather assume that boundary data indices are given on global domain
798                     ! TO BE DISCUSSED ?
799!                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1
800!                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
801                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+jpizoom
802                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+jpjzoom
803                     ! check if point has to be sent
804                     ii = idx_bdy(ib_bdy)%nbi(icount,igrd)
805                     ij = idx_bdy(ib_bdy)%nbj(icount,igrd)
806                     if((com_east .ne. 1) .and. (ii .eq. (nlci-1)) .and. (nbondi .le. 0)) then
807                        com_east = 1
808                     elseif((com_west .ne. 1) .and. (ii .eq. 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then
809                        com_west = 1
810                     endif
811                     if((com_south .ne. 1) .and. (ij .eq. 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then
812                        com_south = 1
813                     elseif((com_north .ne. 1) .and. (ij .eq. (nlcj-1)) .and. (nbondj .le. 0)) then
814                        com_north = 1
815                     endif
816                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
817                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
818                  ENDIF
819                  ! check if point has to be received from a neighbour
820                  IF(nbondi .eq. 0) THEN
821                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   &
822                       & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   &
823                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
824                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
825                       if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then
826                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
827                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
828                            com_south = 1
829                          elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
830                            com_north = 1
831                          endif
832                          com_west_b = 1
833                       endif
834                     ENDIF
835                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND.   &
836                       & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND.   &
837                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
838                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
839                       if((com_east_b .ne. 1) .and. (ii .eq. 2)) then
840                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
841                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
842                            com_south = 1
843                          elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
844                            com_north = 1
845                          endif
846                          com_east_b = 1
847                       endif
848                     ENDIF
849                  ELSEIF(nbondi .eq. 1) THEN
850                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   &
851                       & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   &
852                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
853                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
854                       if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then
855                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
856                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
857                            com_south = 1
858                          elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
859                            com_north = 1
860                          endif
861                          com_west_b = 1
862                       endif
863                     ENDIF
864                  ELSEIF(nbondi .eq. -1) THEN
865                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND.   &
866                       & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND.   &
867                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
868                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
869                       if((com_east_b .ne. 1) .and. (ii .eq. 2)) then
870                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
871                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
872                            com_south = 1
873                          elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
874                            com_north = 1
875                          endif
876                          com_east_b = 1
877                       endif
878                     ENDIF
879                  ENDIF
880                  IF(nbondj .eq. 0) THEN
881                     IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1  &
882                       & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
883                       & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
884                       com_north_b = 1 
885                     ENDIF
886                     IF(com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1  &
887                       &.OR. nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
888                       & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
889                       com_south_b = 1 
890                     ENDIF
891                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND.   &
892                       & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND.   &
893                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
894                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
895                       if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then
896                          com_south_b = 1
897                       endif
898                     ENDIF
899                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND.   &
900                       & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND.   &
901                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
902                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
903                       if((com_north_b .ne. 1) .and. (ij .eq. 2)) then
904                          com_north_b = 1
905                       endif
906                     ENDIF
907                  ELSEIF(nbondj .eq. 1) THEN
908                     IF( com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. &
909                       & nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
910                       & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
911                       com_south_b = 1 
912                     ENDIF
913                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND.   &
914                       & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND.   &
915                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
916                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
917                       if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then
918                          com_south_b = 1
919                       endif
920                     ENDIF
921                  ELSEIF(nbondj .eq. -1) THEN
922                     IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1  &
923                       & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
924                       & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
925                       com_north_b = 1 
926                     ENDIF
927                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND.   &
928                       & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND.   &
929                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
930                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
931                       if((com_north_b .ne. 1) .and. (ij .eq. 2)) then
932                          com_north_b = 1
933                       endif
934                     ENDIF
935                  ENDIF
936               ENDDO
937            ENDDO
938         ENDDO 
939         ! definition of the i- and j- direction local boundaries arrays
940         ! used for sending the boudaries
941         IF((com_east .eq. 1) .and. (com_west .eq. 1)) THEN
942            nbondi_bdy(ib_bdy) = 0
943         ELSEIF ((com_east .eq. 1) .and. (com_west .eq. 0)) THEN
944            nbondi_bdy(ib_bdy) = -1
945         ELSEIF ((com_east .eq. 0) .and. (com_west .eq. 1)) THEN
946            nbondi_bdy(ib_bdy) = 1
947         ENDIF
948
949         IF((com_north .eq. 1) .and. (com_south .eq. 1)) THEN
950            nbondj_bdy(ib_bdy) = 0
951         ELSEIF ((com_north .eq. 1) .and. (com_south .eq. 0)) THEN
952            nbondj_bdy(ib_bdy) = -1
953         ELSEIF ((com_north .eq. 0) .and. (com_south .eq. 1)) THEN
954            nbondj_bdy(ib_bdy) = 1
955         ENDIF
956
957         ! definition of the i- and j- direction local boundaries arrays
958         ! used for receiving the boudaries
959         IF((com_east_b .eq. 1) .and. (com_west_b .eq. 1)) THEN
960            nbondi_bdy_b(ib_bdy) = 0
961         ELSEIF ((com_east_b .eq. 1) .and. (com_west_b .eq. 0)) THEN
962            nbondi_bdy_b(ib_bdy) = -1
963         ELSEIF ((com_east_b .eq. 0) .and. (com_west_b .eq. 1)) THEN
964            nbondi_bdy_b(ib_bdy) = 1
965         ENDIF
966
967         IF((com_north_b .eq. 1) .and. (com_south_b .eq. 1)) THEN
968            nbondj_bdy_b(ib_bdy) = 0
969         ELSEIF ((com_north_b .eq. 1) .and. (com_south_b .eq. 0)) THEN
970            nbondj_bdy_b(ib_bdy) = -1
971         ELSEIF ((com_north_b .eq. 0) .and. (com_south_b .eq. 1)) THEN
972            nbondj_bdy_b(ib_bdy) = 1
973         ENDIF
974
975         ! Compute rim weights for FRS scheme
976         ! ----------------------------------
977         DO igrd = 1, jpbgrd
978            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
979               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
980               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
981!               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.  ! quadratic
982!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy))       ! linear
983            END DO
984         END DO 
985
986         ! Compute damping coefficients
987         ! ----------------------------
988         DO igrd = 1, jpbgrd
989            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
990               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
991               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 
992               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic
993            END DO
994         END DO
995
996      ENDDO
997
998      ! ------------------------------------------------------
999      ! Initialise masks and find normal/tangential directions
1000      ! ------------------------------------------------------
1001
1002      ! Read global 2D mask at T-points: bdytmask
1003      ! -----------------------------------------
1004      ! bdytmask = 1  on the computational domain AND on open boundaries
1005      !          = 0  elsewhere   
1006 
1007      IF( ln_mask_file ) THEN
1008         CALL iom_open( cn_mask_file, inum )
1009         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
1010         CALL iom_close( inum )
1011
1012         ! Derive mask on U and V grid from mask on T grid
1013         bdyumask(:,:) = 0.e0
1014         bdyvmask(:,:) = 0.e0
1015         DO ij=1, jpjm1
1016            DO ii=1, jpim1
1017               bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
1018               bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1) 
1019            END DO
1020         END DO
1021         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
1022
1023
1024         ! Mask corrections
1025         ! ----------------
1026         DO ik = 1, jpkm1
1027            DO ij = 1, jpj
1028               DO ii = 1, jpi
1029                  tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
1030                  umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
1031                  vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
1032                  bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij)
1033               END DO     
1034            END DO
1035         END DO
1036
1037         DO ik = 1, jpkm1
1038            DO ij = 2, jpjm1
1039               DO ii = 2, jpim1
1040                  fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   &
1041                     &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
1042               END DO     
1043            END DO
1044         END DO
1045
1046         tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)
1047
1048      ENDIF ! ln_mask_file=.TRUE.
1049     
1050      bdytmask(:,:) = tmask(:,:,1)
1051
1052      ! bdy masks and bmask are now set to zero on boundary points:
1053      igrd = 1       ! In the free surface case, bmask is at T-points
1054      DO ib_bdy = 1, nb_bdy
1055        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
1056          bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1057        ENDDO
1058      ENDDO
1059      !
1060      igrd = 1
1061      DO ib_bdy = 1, nb_bdy
1062        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
1063          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1064        ENDDO
1065      ENDDO
1066      !
1067      igrd = 2
1068      DO ib_bdy = 1, nb_bdy
1069        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1070          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1071        ENDDO
1072      ENDDO
1073      !
1074      igrd = 3
1075      DO ib_bdy = 1, nb_bdy
1076        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1077          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1078        ENDDO
1079      ENDDO
1080
1081      ! Lateral boundary conditions
1082      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
1083      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
1084
1085      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
1086
1087         idx_bdy(ib_bdy)%flagu(:) = 0.e0
1088         idx_bdy(ib_bdy)%flagv(:) = 0.e0
1089         icount = 0 
1090
1091         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
1092         !flagu =  0 : u is tangential
1093         !flagu =  1 : u is normal to the boundary and is direction is inward
1094 
1095         igrd = 2      ! u-component
1096         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
1097            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1098            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1099            zefl = bdytmask(nbi  ,nbj)
1100            zwfl = bdytmask(nbi+1,nbj)
1101            IF( zefl + zwfl == 2 ) THEN
1102               icount = icount + 1
1103            ELSE
1104               idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl
1105            ENDIF
1106         END DO
1107
1108         !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
1109         !flagv =  0 : u is tangential
1110         !flagv =  1 : u is normal to the boundary and is direction is inward
1111
1112         igrd = 3      ! v-component
1113         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
1114            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1115            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1116            znfl = bdytmask(nbi,nbj  )
1117            zsfl = bdytmask(nbi,nbj+1)
1118            IF( znfl + zsfl == 2 ) THEN
1119               icount = icount + 1
1120            ELSE
1121               idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl
1122            END IF
1123         END DO
1124
1125         IF( icount /= 0 ) THEN
1126            IF(lwp) WRITE(numout,*)
1127            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   &
1128               ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy
1129            IF(lwp) WRITE(numout,*) ' ========== '
1130            IF(lwp) WRITE(numout,*)
1131            nstop = nstop + 1
1132         ENDIF
1133   
1134      ENDDO
1135
1136      ! Compute total lateral surface for volume correction:
1137      ! ----------------------------------------------------
1138      ! JC: this must be done at each time step with key_vvl
1139      bdysurftot = 0.e0 
1140      IF( ln_vol ) THEN 
1141         igrd = 2      ! Lateral surface at U-points
1142         DO ib_bdy = 1, nb_bdy
1143            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1144               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1145               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1146               flagu => idx_bdy(ib_bdy)%flagu(ib)
1147               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           &
1148                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
1149                  &                    * tmask_i(nbi  , nbj)                           &
1150                  &                    * tmask_i(nbi+1, nbj)                   
1151            ENDDO
1152         ENDDO
1153
1154         igrd=3 ! Add lateral surface at V-points
1155         DO ib_bdy = 1, nb_bdy
1156            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1157               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1158               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1159               flagv => idx_bdy(ib_bdy)%flagv(ib)
1160               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           &
1161                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
1162                  &                    * tmask_i(nbi, nbj  )                           &
1163                  &                    * tmask_i(nbi, nbj+1)
1164            ENDDO
1165         ENDDO
1166         !
1167         IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain
1168      END IF   
1169      !
1170      ! Tidy up
1171      !--------
1172      IF (nb_bdy>0) THEN
1173         DEALLOCATE(nbidta, nbjdta, nbrdta)
1174      ENDIF
1175
1176      IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
1177
1178   END SUBROUTINE bdy_init
1179
1180   SUBROUTINE bdy_ctl_seg
1181      !!----------------------------------------------------------------------
1182      !!                 ***  ROUTINE bdy_ctl_seg  ***
1183      !!
1184      !! ** Purpose :   Check straight open boundary segments location
1185      !!
1186      !! ** Method  :   - Look for open boundary corners
1187      !!                - Check that segments start or end on land
1188      !!----------------------------------------------------------------------
1189      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest 
1190      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns 
1191      REAL(wp), DIMENSION(2) ::   ztestmask
1192      !!----------------------------------------------------------------------
1193      !
1194      IF (lwp) WRITE(numout,*) ' '
1195      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1196      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1197      !
1198      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege
1199      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw
1200      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn
1201      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs
1202      ! 1. Check bounds
1203      !----------------
1204      DO ib = 1, nbdysegn
1205         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
1206         IF ((jpjnob(ib).ge.jpjglo-1).or.& 
1207            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1208         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1209         IF (jpindt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1210         IF (jpinft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
1211      END DO
1212      !
1213      DO ib = 1, nbdysegs
1214         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
1215         IF ((jpjsob(ib).ge.jpjglo-1).or.& 
1216            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1217         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1218         IF (jpisdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1219         IF (jpisft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
1220      END DO
1221      !
1222      DO ib = 1, nbdysege
1223         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib)
1224         IF ((jpieob(ib).ge.jpiglo-1).or.& 
1225            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1226         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1227         IF (jpjedt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1228         IF (jpjeft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
1229      END DO
1230      !
1231      DO ib = 1, nbdysegw
1232         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib)
1233         IF ((jpiwob(ib).ge.jpiglo-1).or.& 
1234            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1235         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1236         IF (jpjwdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1237         IF (jpjwft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
1238      ENDDO
1239      !
1240      !     
1241      ! 2. Look for segment crossings
1242      !------------------------------
1243      IF (lwp) WRITE(numout,*) '**Look for segments corners  :'
1244      !
1245      itest = 0 ! corner number
1246      !
1247      ! flag to detect if start or end of open boundary belongs to a corner
1248      ! if not (=0), it must be on land.
1249      ! if a corner is detected, save bdy package number for further tests
1250      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1251      ! South/West crossings
1252      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1253         DO ib1 = 1, nbdysegw       
1254            DO ib2 = 1, nbdysegs
1255               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1256                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1257                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1258                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1259                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN 
1260                     ! We have a possible South-West corner                     
1261!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1262!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1263                     icornw(ib1,1) = npckgs(ib2)
1264                     icorns(ib2,1) = npckgw(ib1)
1265                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1266                     IF(lwp) WRITE(numout,*)
1267                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1268                     &                                     jpisft(ib2), jpjwft(ib1)
1269                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1270                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), & 
1271                     &                                                    ' and South segment: ',npckgs(ib2)
1272                     IF(lwp) WRITE(numout,*)
1273                     nstop = nstop + 1
1274                  ELSE
1275                     IF(lwp) WRITE(numout,*)
1276                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices'
1277                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1) , &
1278                     &                                                    ' and South segment: ',npckgs(ib2)
1279                     IF(lwp) WRITE(numout,*)
1280                     nstop = nstop+1
1281                  END IF
1282               END IF
1283            END DO
1284         END DO
1285      END IF
1286      !
1287      ! South/East crossings
1288      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1289         DO ib1 = 1, nbdysege
1290            DO ib2 = 1, nbdysegs
1291               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1292                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1293                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. &
1294                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN
1295                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1296                     ! We have a possible South-East corner
1297!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1298!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1299                     icorne(ib1,1) = npckgs(ib2)
1300                     icorns(ib2,2) = npckge(ib1)
1301                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1302                     IF(lwp) WRITE(numout,*)
1303                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1304                     &                                     jpisdt(ib2), jpjeft(ib1)
1305                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1306                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), &
1307                     &                                                    ' and South segment: ',npckgs(ib2)
1308                     IF(lwp) WRITE(numout,*)
1309                     nstop = nstop + 1
1310                  ELSE
1311                     IF(lwp) WRITE(numout,*)
1312                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices'
1313                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
1314                     &                                                    ' and South segment: ',npckgs(ib2)
1315                     IF(lwp) WRITE(numout,*)
1316                     nstop = nstop + 1
1317                  END IF
1318               END IF
1319            END DO
1320         END DO
1321      END IF
1322      !
1323      ! North/West crossings
1324      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1325         DO ib1 = 1, nbdysegw       
1326            DO ib2 = 1, nbdysegn
1327               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. &
1328                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. &
1329                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1330                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1331                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1332                     ! We have a possible North-West corner
1333!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1334!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1335                     icornw(ib1,2) = npckgn(ib2)
1336                     icornn(ib2,1) = npckgw(ib1)
1337                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1338                     IF(lwp) WRITE(numout,*)
1339                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1340                     &                                     jpinft(ib2), jpjwdt(ib1)
1341                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1342                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), &
1343                     &                                                    ' and North segment: ',npckgn(ib2)
1344                     IF(lwp) WRITE(numout,*)
1345                     nstop = nstop + 1
1346                  ELSE
1347                     IF(lwp) WRITE(numout,*)
1348                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices'
1349                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1), &
1350                     &                                                    ' and North segment: ',npckgn(ib2)
1351                     IF(lwp) WRITE(numout,*)
1352                     nstop = nstop + 1
1353                  END IF
1354               END IF
1355            END DO
1356         END DO
1357      END IF
1358      !
1359      ! North/East crossings
1360      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1361         DO ib1 = 1, nbdysege       
1362            DO ib2 = 1, nbdysegn
1363               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1364                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1365                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1366                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1367                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1368                     ! We have a possible North-East corner
1369!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1370!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1371                     icorne(ib1,2) = npckgn(ib2)
1372                     icornn(ib2,2) = npckge(ib1)
1373                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1374                     IF(lwp) WRITE(numout,*)
1375                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1376                     &                                     jpindt(ib2), jpjedt(ib1)
1377                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1378                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), &
1379                     &                                                    ' and North segment: ',npckgn(ib2)
1380                     IF(lwp) WRITE(numout,*)
1381                     nstop = nstop + 1
1382                  ELSE
1383                     IF(lwp) WRITE(numout,*)
1384                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices'
1385                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
1386                     &                                                    ' and North segment: ',npckgn(ib2)
1387                     IF(lwp) WRITE(numout,*)
1388                     nstop = nstop + 1
1389                  END IF
1390               END IF
1391            END DO
1392         END DO
1393      END IF
1394      !
1395      ! 3. Check if segment extremities are on land
1396      !--------------------------------------------
1397      !
1398      ! West segments
1399      DO ib = 1, nbdysegw
1400         ! get mask at boundary extremities:
1401         ztestmask(1:2)=0.
1402         DO ji = 1, jpi
1403            DO jj = 1, jpj             
1404              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1405               &  ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1406              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1407               &  ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1408            END DO
1409         END DO
1410         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1411
1412         IF (ztestmask(1)==1) THEN
1413            IF (icornw(ib,1)==0) THEN
1414               IF(lwp) WRITE(numout,*)
1415               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
1416               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1417               IF(lwp) WRITE(numout,*)
1418               nstop = nstop + 1
1419            ELSE
1420               ! This is a corner
1421               WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1422               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1423               itest=itest+1
1424            ENDIF
1425         ENDIF
1426         IF (ztestmask(2)==1) THEN
1427            IF (icornw(ib,2)==0) THEN
1428               IF(lwp) WRITE(numout,*)
1429               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
1430               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1431               IF(lwp) WRITE(numout,*)
1432               nstop = nstop + 1
1433            ELSE
1434               ! This is a corner
1435               WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1436               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1437               itest=itest+1
1438            ENDIF
1439         ENDIF
1440      END DO
1441      !
1442      ! East segments
1443      DO ib = 1, nbdysege
1444         ! get mask at boundary extremities:
1445         ztestmask(1:2)=0.
1446         DO ji = 1, jpi
1447            DO jj = 1, jpj             
1448              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1449               &  ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
1450              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1451               &  ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1452            END DO
1453         END DO
1454         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1455
1456         IF (ztestmask(1)==1) THEN
1457            IF (icorne(ib,1)==0) THEN
1458               IF(lwp) WRITE(numout,*)
1459               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
1460               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1461               IF(lwp) WRITE(numout,*)
1462               nstop = nstop + 1 
1463            ELSE
1464               ! This is a corner
1465               WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1466               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1467               itest=itest+1
1468            ENDIF
1469         ENDIF
1470         IF (ztestmask(2)==1) THEN
1471            IF (icorne(ib,2)==0) THEN
1472               IF(lwp) WRITE(numout,*)
1473               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
1474               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1475               IF(lwp) WRITE(numout,*)
1476               nstop = nstop + 1
1477            ELSE
1478               ! This is a corner
1479               WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1480               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1481               itest=itest+1
1482            ENDIF
1483         ENDIF
1484      END DO
1485      !
1486      ! South segments
1487      DO ib = 1, nbdysegs
1488         ! get mask at boundary extremities:
1489         ztestmask(1:2)=0.
1490         DO ji = 1, jpi
1491            DO jj = 1, jpj             
1492              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1493               &  ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1494              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1495               &  ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1496            END DO
1497         END DO
1498         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1499
1500         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1501            IF(lwp) WRITE(numout,*)
1502            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
1503            IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1504            IF(lwp) WRITE(numout,*)
1505            nstop = nstop + 1
1506         ENDIF
1507         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1508            IF(lwp) WRITE(numout,*)
1509            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
1510            IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1511            IF(lwp) WRITE(numout,*)
1512            nstop = nstop + 1
1513         ENDIF
1514      END DO
1515      !
1516      ! North segments
1517      DO ib = 1, nbdysegn
1518         ! get mask at boundary extremities:
1519         ztestmask(1:2)=0.
1520         DO ji = 1, jpi
1521            DO jj = 1, jpj             
1522              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1523               &  ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
1524              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1525               &  ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1526            END DO
1527         END DO
1528         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1529
1530         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1531            IF(lwp) WRITE(numout,*)
1532            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
1533            IF(lwp) WRITE(numout,*) ' ==========  does not start on land'                                                 
1534            IF(lwp) WRITE(numout,*)
1535            nstop = nstop + 1
1536         ENDIF
1537         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1538            IF(lwp) WRITE(numout,*)
1539            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
1540            IF(lwp) WRITE(numout,*) ' ==========  does not end on land'                                                 
1541            IF(lwp) WRITE(numout,*)
1542            nstop = nstop + 1
1543         ENDIF
1544      END DO
1545      !
1546      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1547      !
1548      ! Other tests TBD:
1549      ! segments completly on land
1550      ! optimized open boundary array length according to landmask
1551      ! Nudging layers that overlap with interior domain
1552      !
1553   END SUBROUTINE bdy_ctl_seg
1554
1555   SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1556      !!----------------------------------------------------------------------
1557      !!                 ***  ROUTINE bdy_ctl_corn  ***
1558      !!
1559      !! ** Purpose :   Check numerical schemes consistency between
1560      !!                segments having a common corner
1561      !!
1562      !! ** Method  :   
1563      !!----------------------------------------------------------------------
1564      INTEGER, INTENT(in)  ::   ib1, ib2
1565      INTEGER :: itest
1566      !!----------------------------------------------------------------------
1567      itest = 0
1568
1569      IF (nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 1
1570      IF (nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 1
1571      IF (nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 1
1572      !
1573      IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1
1574      IF (nn_dyn3d_dta(ib1)/=nn_dyn3d_dta(ib2)) itest = itest + 1
1575      IF (nn_tra_dta(ib1)/=nn_tra_dta(ib2)) itest = itest + 1
1576      !
1577      IF (nn_rimwidth(ib1)/=nn_rimwidth(ib2)) itest = itest + 1   
1578      !
1579      IF ( itest>0 ) THEN
1580         IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2
1581         IF(lwp) WRITE(numout,*) ' ==========  have different open bdy schemes'                                                 
1582         IF(lwp) WRITE(numout,*)
1583         nstop = nstop + 1
1584      ENDIF
1585      !
1586   END SUBROUTINE bdy_ctl_corn
1587
1588#else
1589   !!---------------------------------------------------------------------------------
1590   !!   Dummy module                                   NO open boundaries
1591   !!---------------------------------------------------------------------------------
1592CONTAINS
1593   SUBROUTINE bdy_init      ! Dummy routine
1594   END SUBROUTINE bdy_init
1595#endif
1596
1597   !!=================================================================================
1598END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.