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 NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90 @ 11049

Last change on this file since 11049 was 11049, checked in by girrmann, 5 years ago
dev_r10984_HPC-13 : CYCLE instruction is not systematic anymore, computation is done on the halo whenever possible and overwritten by lbc_bdy instruction, see #2285
Property svn:keywords set to `Id`
File size: 82.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: