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, 7 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.