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

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

source: branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 5878

Last change on this file since 5878 was 5878, checked in by deazer, 8 years ago

Adds simple ramp of diffusion coefficient in BDY rim using TANH weighting

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