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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

bdyini.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 4292

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

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

Download in other formats: