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

Context Navigation

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 5626

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

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

Download in other formats: