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

Context Navigation

source: branches/2016/dev_r6522_SIMPLIF_3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 6864

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

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

Download in other formats: