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_merge_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/bdyini.F90 @ 9125

Last change on this file since 9125 was 9125, checked in by timgraham, 6 years ago
Removed wrk_arrays from whole code. No change in SETTE results from this.
File size: 82.5 KB

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

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

Download in other formats: