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

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 83.9 KB

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

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

Download in other formats: