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/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 9169

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

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

Download in other formats: