source: NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/BDY/bdyini.F90 @ 10455

Last change on this file since 10455 was 10455, checked in by mathiot, 22 months ago

suggestion for ticket #2200 (bdyvol)

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