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 @ 11183

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

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

Download in other formats: