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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

bdyini.F90 in branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 8882

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

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

Download in other formats: