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_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 4291

Last change on this file since 4291 was 4148, checked in by cetlod, 10 years ago

merge in trunk changes between r3853 and r3940 and commit the changes, see ticket #1169

  • Property svn:keywords set to Id
File size: 75.1 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      IF( .not. ln_mask_file ) THEN
1052         ! If .not. ln_mask_file then we need to derive mask on U and V grid
1053         ! from mask on T grid here.
1054         bdyumask(:,:) = 0.e0
1055         bdyvmask(:,:) = 0.e0
1056         DO ij=1, jpjm1
1057            DO ii=1, jpim1
1058               bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
1059               bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1) 
1060            END DO
1061         END DO
1062         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
1063      ENDIF
1064
1065      ! bdy masks and bmask are now set to zero on boundary points:
1066      igrd = 1       ! In the free surface case, bmask is at T-points
1067      DO ib_bdy = 1, nb_bdy
1068        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
1069          bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1070        ENDDO
1071      ENDDO
1072      !
1073      igrd = 1
1074      DO ib_bdy = 1, nb_bdy
1075        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
1076          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1077        ENDDO
1078      ENDDO
1079      !
1080      igrd = 2
1081      DO ib_bdy = 1, nb_bdy
1082        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1083          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1084        ENDDO
1085      ENDDO
1086      !
1087      igrd = 3
1088      DO ib_bdy = 1, nb_bdy
1089        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1090          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1091        ENDDO
1092      ENDDO
1093
1094      ! Lateral boundary conditions
1095      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
1096      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
1097
1098      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
1099
1100         idx_bdy(ib_bdy)%flagu(:) = 0.e0
1101         idx_bdy(ib_bdy)%flagv(:) = 0.e0
1102         icount = 0 
1103
1104         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
1105         !flagu =  0 : u is tangential
1106         !flagu =  1 : u is normal to the boundary and is direction is inward
1107 
1108         igrd = 2      ! u-component
1109         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
1110            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1111            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1112            zefl = bdytmask(nbi  ,nbj)
1113            zwfl = bdytmask(nbi+1,nbj)
1114            IF( zefl + zwfl == 2 ) THEN
1115               icount = icount + 1
1116            ELSE
1117               idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl
1118            ENDIF
1119         END DO
1120
1121         !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
1122         !flagv =  0 : u is tangential
1123         !flagv =  1 : u is normal to the boundary and is direction is inward
1124
1125         igrd = 3      ! v-component
1126         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
1127            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1128            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1129            znfl = bdytmask(nbi,nbj  )
1130            zsfl = bdytmask(nbi,nbj+1)
1131            IF( znfl + zsfl == 2 ) THEN
1132               icount = icount + 1
1133            ELSE
1134               idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl
1135            END IF
1136         END DO
1137
1138         IF( icount /= 0 ) THEN
1139            IF(lwp) WRITE(numout,*)
1140            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   &
1141               ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy
1142            IF(lwp) WRITE(numout,*) ' ========== '
1143            IF(lwp) WRITE(numout,*)
1144            nstop = nstop + 1
1145         ENDIF
1146   
1147      ENDDO
1148
1149      ! Compute total lateral surface for volume correction:
1150      ! ----------------------------------------------------
1151      ! JC: this must be done at each time step with key_vvl
1152      bdysurftot = 0.e0 
1153      IF( ln_vol ) THEN 
1154         igrd = 2      ! Lateral surface at U-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               flagu => idx_bdy(ib_bdy)%flagu(ib)
1160               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           &
1161                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
1162                  &                    * tmask_i(nbi  , nbj)                           &
1163                  &                    * tmask_i(nbi+1, nbj)                   
1164            ENDDO
1165         ENDDO
1166
1167         igrd=3 ! Add lateral surface at V-points
1168         DO ib_bdy = 1, nb_bdy
1169            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1170               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1171               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1172               flagv => idx_bdy(ib_bdy)%flagv(ib)
1173               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           &
1174                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
1175                  &                    * tmask_i(nbi, nbj  )                           &
1176                  &                    * tmask_i(nbi, nbj+1)
1177            ENDDO
1178         ENDDO
1179         !
1180         IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain
1181      END IF   
1182      !
1183      ! Tidy up
1184      !--------
1185      IF (nb_bdy>0) THEN
1186         DEALLOCATE(nbidta, nbjdta, nbrdta)
1187      ENDIF
1188
1189      IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
1190
1191   END SUBROUTINE bdy_init
1192
1193   SUBROUTINE bdy_ctl_seg
1194      !!----------------------------------------------------------------------
1195      !!                 ***  ROUTINE bdy_ctl_seg  ***
1196      !!
1197      !! ** Purpose :   Check straight open boundary segments location
1198      !!
1199      !! ** Method  :   - Look for open boundary corners
1200      !!                - Check that segments start or end on land
1201      !!----------------------------------------------------------------------
1202      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest 
1203      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns 
1204      REAL(wp), DIMENSION(2) ::   ztestmask
1205      !!----------------------------------------------------------------------
1206      !
1207      IF (lwp) WRITE(numout,*) ' '
1208      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1209      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1210      !
1211      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege
1212      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw
1213      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn
1214      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs
1215      ! 1. Check bounds
1216      !----------------
1217      DO ib = 1, nbdysegn
1218         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
1219         IF ((jpjnob(ib).ge.jpjglo-1).or.& 
1220            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1221         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1222         IF (jpindt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1223         IF (jpinft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
1224      END DO
1225      !
1226      DO ib = 1, nbdysegs
1227         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
1228         IF ((jpjsob(ib).ge.jpjglo-1).or.& 
1229            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1230         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1231         IF (jpisdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1232         IF (jpisft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
1233      END DO
1234      !
1235      DO ib = 1, nbdysege
1236         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib)
1237         IF ((jpieob(ib).ge.jpiglo-1).or.& 
1238            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1239         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1240         IF (jpjedt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1241         IF (jpjeft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
1242      END DO
1243      !
1244      DO ib = 1, nbdysegw
1245         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib)
1246         IF ((jpiwob(ib).ge.jpiglo-1).or.& 
1247            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1248         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1249         IF (jpjwdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1250         IF (jpjwft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
1251      ENDDO
1252      !
1253      !     
1254      ! 2. Look for segment crossings
1255      !------------------------------
1256      IF (lwp) WRITE(numout,*) '**Look for segments corners  :'
1257      !
1258      itest = 0 ! corner number
1259      !
1260      ! flag to detect if start or end of open boundary belongs to a corner
1261      ! if not (=0), it must be on land.
1262      ! if a corner is detected, save bdy package number for further tests
1263      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1264      ! South/West crossings
1265      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1266         DO ib1 = 1, nbdysegw       
1267            DO ib2 = 1, nbdysegs
1268               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1269                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1270                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1271                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1272                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN 
1273                     ! We have a possible South-West corner                     
1274!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1275!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1276                     icornw(ib1,1) = npckgs(ib2)
1277                     icorns(ib2,1) = npckgw(ib1)
1278                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1279                     IF(lwp) WRITE(numout,*)
1280                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1281                     &                                     jpisft(ib2), jpjwft(ib1)
1282                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1283                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), & 
1284                     &                                                    ' and South segment: ',npckgs(ib2)
1285                     IF(lwp) WRITE(numout,*)
1286                     nstop = nstop + 1
1287                  ELSE
1288                     IF(lwp) WRITE(numout,*)
1289                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices'
1290                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1) , &
1291                     &                                                    ' and South segment: ',npckgs(ib2)
1292                     IF(lwp) WRITE(numout,*)
1293                     nstop = nstop+1
1294                  END IF
1295               END IF
1296            END DO
1297         END DO
1298      END IF
1299      !
1300      ! South/East crossings
1301      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1302         DO ib1 = 1, nbdysege
1303            DO ib2 = 1, nbdysegs
1304               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1305                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1306                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. &
1307                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN
1308                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1309                     ! We have a possible South-East corner
1310!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1311!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1312                     icorne(ib1,1) = npckgs(ib2)
1313                     icorns(ib2,2) = npckge(ib1)
1314                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1315                     IF(lwp) WRITE(numout,*)
1316                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1317                     &                                     jpisdt(ib2), jpjeft(ib1)
1318                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1319                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), &
1320                     &                                                    ' and South segment: ',npckgs(ib2)
1321                     IF(lwp) WRITE(numout,*)
1322                     nstop = nstop + 1
1323                  ELSE
1324                     IF(lwp) WRITE(numout,*)
1325                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices'
1326                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
1327                     &                                                    ' and South segment: ',npckgs(ib2)
1328                     IF(lwp) WRITE(numout,*)
1329                     nstop = nstop + 1
1330                  END IF
1331               END IF
1332            END DO
1333         END DO
1334      END IF
1335      !
1336      ! North/West crossings
1337      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1338         DO ib1 = 1, nbdysegw       
1339            DO ib2 = 1, nbdysegn
1340               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. &
1341                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. &
1342                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1343                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1344                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1345                     ! We have a possible North-West corner
1346!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1347!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1348                     icornw(ib1,2) = npckgn(ib2)
1349                     icornn(ib2,1) = npckgw(ib1)
1350                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1351                     IF(lwp) WRITE(numout,*)
1352                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1353                     &                                     jpinft(ib2), jpjwdt(ib1)
1354                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1355                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), &
1356                     &                                                    ' and North segment: ',npckgn(ib2)
1357                     IF(lwp) WRITE(numout,*)
1358                     nstop = nstop + 1
1359                  ELSE
1360                     IF(lwp) WRITE(numout,*)
1361                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices'
1362                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1), &
1363                     &                                                    ' and North segment: ',npckgn(ib2)
1364                     IF(lwp) WRITE(numout,*)
1365                     nstop = nstop + 1
1366                  END IF
1367               END IF
1368            END DO
1369         END DO
1370      END IF
1371      !
1372      ! North/East crossings
1373      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1374         DO ib1 = 1, nbdysege       
1375            DO ib2 = 1, nbdysegn
1376               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1377                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1378                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1379                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1380                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1381                     ! We have a possible North-East corner
1382!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1383!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1384                     icorne(ib1,2) = npckgn(ib2)
1385                     icornn(ib2,2) = npckge(ib1)
1386                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1387                     IF(lwp) WRITE(numout,*)
1388                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1389                     &                                     jpindt(ib2), jpjedt(ib1)
1390                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1391                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), &
1392                     &                                                    ' and North segment: ',npckgn(ib2)
1393                     IF(lwp) WRITE(numout,*)
1394                     nstop = nstop + 1
1395                  ELSE
1396                     IF(lwp) WRITE(numout,*)
1397                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices'
1398                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
1399                     &                                                    ' and North segment: ',npckgn(ib2)
1400                     IF(lwp) WRITE(numout,*)
1401                     nstop = nstop + 1
1402                  END IF
1403               END IF
1404            END DO
1405         END DO
1406      END IF
1407      !
1408      ! 3. Check if segment extremities are on land
1409      !--------------------------------------------
1410      !
1411      ! West segments
1412      DO ib = 1, nbdysegw
1413         ! get mask at boundary extremities:
1414         ztestmask(1:2)=0.
1415         DO ji = 1, jpi
1416            DO jj = 1, jpj             
1417              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1418               &  ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1419              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1420               &  ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1421            END DO
1422         END DO
1423         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1424
1425         IF (ztestmask(1)==1) THEN
1426            IF (icornw(ib,1)==0) THEN
1427               IF(lwp) WRITE(numout,*)
1428               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
1429               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1430               IF(lwp) WRITE(numout,*)
1431               nstop = nstop + 1
1432            ELSE
1433               ! This is a corner
1434               WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1435               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1436               itest=itest+1
1437            ENDIF
1438         ENDIF
1439         IF (ztestmask(2)==1) THEN
1440            IF (icornw(ib,2)==0) THEN
1441               IF(lwp) WRITE(numout,*)
1442               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
1443               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1444               IF(lwp) WRITE(numout,*)
1445               nstop = nstop + 1
1446            ELSE
1447               ! This is a corner
1448               WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1449               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1450               itest=itest+1
1451            ENDIF
1452         ENDIF
1453      END DO
1454      !
1455      ! East segments
1456      DO ib = 1, nbdysege
1457         ! get mask at boundary extremities:
1458         ztestmask(1:2)=0.
1459         DO ji = 1, jpi
1460            DO jj = 1, jpj             
1461              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1462               &  ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
1463              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1464               &  ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1465            END DO
1466         END DO
1467         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1468
1469         IF (ztestmask(1)==1) THEN
1470            IF (icorne(ib,1)==0) THEN
1471               IF(lwp) WRITE(numout,*)
1472               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
1473               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1474               IF(lwp) WRITE(numout,*)
1475               nstop = nstop + 1 
1476            ELSE
1477               ! This is a corner
1478               WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1479               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1480               itest=itest+1
1481            ENDIF
1482         ENDIF
1483         IF (ztestmask(2)==1) THEN
1484            IF (icorne(ib,2)==0) THEN
1485               IF(lwp) WRITE(numout,*)
1486               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
1487               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1488               IF(lwp) WRITE(numout,*)
1489               nstop = nstop + 1
1490            ELSE
1491               ! This is a corner
1492               WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1493               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1494               itest=itest+1
1495            ENDIF
1496         ENDIF
1497      END DO
1498      !
1499      ! South segments
1500      DO ib = 1, nbdysegs
1501         ! get mask at boundary extremities:
1502         ztestmask(1:2)=0.
1503         DO ji = 1, jpi
1504            DO jj = 1, jpj             
1505              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1506               &  ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1507              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1508               &  ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1509            END DO
1510         END DO
1511         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1512
1513         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1514            IF(lwp) WRITE(numout,*)
1515            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
1516            IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1517            IF(lwp) WRITE(numout,*)
1518            nstop = nstop + 1
1519         ENDIF
1520         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1521            IF(lwp) WRITE(numout,*)
1522            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
1523            IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1524            IF(lwp) WRITE(numout,*)
1525            nstop = nstop + 1
1526         ENDIF
1527      END DO
1528      !
1529      ! North segments
1530      DO ib = 1, nbdysegn
1531         ! get mask at boundary extremities:
1532         ztestmask(1:2)=0.
1533         DO ji = 1, jpi
1534            DO jj = 1, jpj             
1535              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1536               &  ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
1537              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1538               &  ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1539            END DO
1540         END DO
1541         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1542
1543         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1544            IF(lwp) WRITE(numout,*)
1545            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
1546            IF(lwp) WRITE(numout,*) ' ==========  does not start on land'                                                 
1547            IF(lwp) WRITE(numout,*)
1548            nstop = nstop + 1
1549         ENDIF
1550         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1551            IF(lwp) WRITE(numout,*)
1552            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
1553            IF(lwp) WRITE(numout,*) ' ==========  does not end on land'                                                 
1554            IF(lwp) WRITE(numout,*)
1555            nstop = nstop + 1
1556         ENDIF
1557      END DO
1558      !
1559      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1560      !
1561      ! Other tests TBD:
1562      ! segments completly on land
1563      ! optimized open boundary array length according to landmask
1564      ! Nudging layers that overlap with interior domain
1565      !
1566   END SUBROUTINE bdy_ctl_seg
1567
1568   SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1569      !!----------------------------------------------------------------------
1570      !!                 ***  ROUTINE bdy_ctl_corn  ***
1571      !!
1572      !! ** Purpose :   Check numerical schemes consistency between
1573      !!                segments having a common corner
1574      !!
1575      !! ** Method  :   
1576      !!----------------------------------------------------------------------
1577      INTEGER, INTENT(in)  ::   ib1, ib2
1578      INTEGER :: itest
1579      !!----------------------------------------------------------------------
1580      itest = 0
1581
1582      IF (nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 1
1583      IF (nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 1
1584      IF (nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 1
1585      !
1586      IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1
1587      IF (nn_dyn3d_dta(ib1)/=nn_dyn3d_dta(ib2)) itest = itest + 1
1588      IF (nn_tra_dta(ib1)/=nn_tra_dta(ib2)) itest = itest + 1
1589      !
1590      IF (nn_rimwidth(ib1)/=nn_rimwidth(ib2)) itest = itest + 1   
1591      !
1592      IF ( itest>0 ) THEN
1593         IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2
1594         IF(lwp) WRITE(numout,*) ' ==========  have different open bdy schemes'                                                 
1595         IF(lwp) WRITE(numout,*)
1596         nstop = nstop + 1
1597      ENDIF
1598      !
1599   END SUBROUTINE bdy_ctl_corn
1600
1601#else
1602   !!---------------------------------------------------------------------------------
1603   !!   Dummy module                                   NO open boundaries
1604   !!---------------------------------------------------------------------------------
1605CONTAINS
1606   SUBROUTINE bdy_init      ! Dummy routine
1607   END SUBROUTINE bdy_init
1608#endif
1609
1610   !!=================================================================================
1611END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.