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

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

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

Download in other formats: