Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

bdyini.F90 @ 4743

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

Download in other formats: