Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

bdyini.F90

Last change on this file was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

Download in other formats: