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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago
#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
Property svn:keywords set to `Id`
File size: 82.6 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	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 ', &
848	& ' in order of distance from edge nbr A utility for re-ordering ', &
849	& ' boundary coordinates and data files exists in the TOOLS/OBC directory')
850	ENDIF
851	ENDIF
852	! check if point is in local domain
853	IF( nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND. &
854	& nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in ) THEN
855	!
856	icount = icount + 1
857	!
858	IF( nbrdta(ib,igrd,ib_bdy) == 1 ) icountr = icountr+1
859	ENDIF
860	ENDDO
861	idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
862	idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc
863	ENDDO ! igrd
864
865	! Allocate index arrays for this boundary set
866	!--------------------------------------------
867	ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
868	ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
869	ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
870	ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
871	ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) )
872	ALLOCATE( idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) )
873	ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
874	ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
875	ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1,jpbgrd) )
876	ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1,jpbgrd) )
877
878	! Dispatch mapping indices and discrete distances on each processor
879	! -----------------------------------------------------------------
880
881	com_east = 0
882	com_west = 0
883	com_south = 0
884	com_north = 0
885
886	com_east_b = 0
887	com_west_b = 0
888	com_south_b = 0
889	com_north_b = 0
890	DO igrd = 1, jpbgrd
891	icount = 0
892	! Loop on rimwidth to ensure outermost points come first in the local arrays.
893	DO ir=1, nn_rimwidth(ib_bdy)
894	DO ib = 1, nblendta(igrd,ib_bdy)
895	! check if point is in local domain and equals ir
896	IF( nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND. &
897	& nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND. &
898	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
899	!
900	icount = icount + 1
901
902	! Rather assume that boundary data indices are given on global domain
903	! TO BE DISCUSSED ?
904	! idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1
905	! idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
906	idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+jpizoom
907	idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+jpjzoom
908	! check if point has to be sent
909	ii = idx_bdy(ib_bdy)%nbi(icount,igrd)
910	ij = idx_bdy(ib_bdy)%nbj(icount,igrd)
911	if((com_east .ne. 1) .and. (ii .eq. (nlci-1)) .and. (nbondi .le. 0)) then
912	com_east = 1
913	elseif((com_west .ne. 1) .and. (ii .eq. 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then
914	com_west = 1
915	endif
916	if((com_south .ne. 1) .and. (ij .eq. 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then
917	com_south = 1
918	elseif((com_north .ne. 1) .and. (ij .eq. (nlcj-1)) .and. (nbondj .le. 0)) then
919	com_north = 1
920	endif
921	idx_bdy(ib_bdy)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_bdy)
922	idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
923	ENDIF
924	! check if point has to be received from a neighbour
925	IF(nbondi .eq. 0) THEN
926	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND. &
927	& nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND. &
928	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
929	ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
930	if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then
931	ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
932	if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
933	com_south = 1
934	elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
935	com_north = 1
936	endif
937	com_west_b = 1
938	endif
939	ENDIF
940	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND. &
941	& nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND. &
942	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
943	ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
944	if((com_east_b .ne. 1) .and. (ii .eq. 2)) then
945	ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
946	if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
947	com_south = 1
948	elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
949	com_north = 1
950	endif
951	com_east_b = 1
952	endif
953	ENDIF
954	ELSEIF(nbondi .eq. 1) THEN
955	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND. &
956	& nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND. &
957	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
958	ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
959	if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then
960	ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
961	if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
962	com_south = 1
963	elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
964	com_north = 1
965	endif
966	com_west_b = 1
967	endif
968	ENDIF
969	ELSEIF(nbondi .eq. -1) THEN
970	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND. &
971	& nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND. &
972	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
973	ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
974	if((com_east_b .ne. 1) .and. (ii .eq. 2)) then
975	ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
976	if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
977	com_south = 1
978	elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
979	com_north = 1
980	endif
981	com_east_b = 1
982	endif
983	ENDIF
984	ENDIF
985	IF(nbondj .eq. 0) THEN
986	IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 &
987	& .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
988	& nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
989	com_north_b = 1
990	ENDIF
991	IF(com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 &
992	&.OR. nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
993	& nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
994	com_south_b = 1
995	ENDIF
996	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND. &
997	& nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND. &
998	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
999	ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
1000	if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then
1001	com_south_b = 1
1002	endif
1003	ENDIF
1004	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND. &
1005	& nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND. &
1006	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
1007	ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
1008	if((com_north_b .ne. 1) .and. (ij .eq. 2)) then
1009	com_north_b = 1
1010	endif
1011	ENDIF
1012	ELSEIF(nbondj .eq. 1) THEN
1013	IF( com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. &
1014	& nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
1015	& nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
1016	com_south_b = 1
1017	ENDIF
1018	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND. &
1019	& nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND. &
1020	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
1021	ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
1022	if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then
1023	com_south_b = 1
1024	endif
1025	ENDIF
1026	ELSEIF(nbondj .eq. -1) THEN
1027	IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 &
1028	& .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
1029	& nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
1030	com_north_b = 1
1031	ENDIF
1032	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND. &
1033	& nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND. &
1034	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
1035	ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
1036	if((com_north_b .ne. 1) .and. (ij .eq. 2)) then
1037	com_north_b = 1
1038	endif
1039	ENDIF
1040	ENDIF
1041	ENDDO
1042	ENDDO
1043	ENDDO
1044
1045	! definition of the i- and j- direction local boundaries arrays
1046	! used for sending the boudaries
1047	IF((com_east .eq. 1) .and. (com_west .eq. 1)) THEN
1048	nbondi_bdy(ib_bdy) = 0
1049	ELSEIF ((com_east .eq. 1) .and. (com_west .eq. 0)) THEN
1050	nbondi_bdy(ib_bdy) = -1
1051	ELSEIF ((com_east .eq. 0) .and. (com_west .eq. 1)) THEN
1052	nbondi_bdy(ib_bdy) = 1
1053	ENDIF
1054
1055	IF((com_north .eq. 1) .and. (com_south .eq. 1)) THEN
1056	nbondj_bdy(ib_bdy) = 0
1057	ELSEIF ((com_north .eq. 1) .and. (com_south .eq. 0)) THEN
1058	nbondj_bdy(ib_bdy) = -1
1059	ELSEIF ((com_north .eq. 0) .and. (com_south .eq. 1)) THEN
1060	nbondj_bdy(ib_bdy) = 1
1061	ENDIF
1062
1063	! definition of the i- and j- direction local boundaries arrays
1064	! used for receiving the boudaries
1065	IF((com_east_b .eq. 1) .and. (com_west_b .eq. 1)) THEN
1066	nbondi_bdy_b(ib_bdy) = 0
1067	ELSEIF ((com_east_b .eq. 1) .and. (com_west_b .eq. 0)) THEN
1068	nbondi_bdy_b(ib_bdy) = -1
1069	ELSEIF ((com_east_b .eq. 0) .and. (com_west_b .eq. 1)) THEN
1070	nbondi_bdy_b(ib_bdy) = 1
1071	ENDIF
1072
1073	IF((com_north_b .eq. 1) .and. (com_south_b .eq. 1)) THEN
1074	nbondj_bdy_b(ib_bdy) = 0
1075	ELSEIF ((com_north_b .eq. 1) .and. (com_south_b .eq. 0)) THEN
1076	nbondj_bdy_b(ib_bdy) = -1
1077	ELSEIF ((com_north_b .eq. 0) .and. (com_south_b .eq. 1)) THEN
1078	nbondj_bdy_b(ib_bdy) = 1
1079	ENDIF
1080
1081	! Compute rim weights for FRS scheme
1082	! ----------------------------------
1083	DO igrd = 1, jpbgrd
1084	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
1085	nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
1086	idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation
1087	! idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2. ! quadratic
1088	! idx_bdy(ib_bdy)%nbw(ib,igrd) = FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)) ! linear
1089	END DO
1090	END DO
1091
1092	! Compute damping coefficients
1093	! ----------------------------
1094	DO igrd = 1, jpbgrd
1095	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
1096	nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
1097	idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &
1098	& (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))*2. ! quadratic
1099	idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) &
1100	& (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))*2. ! quadratic
1101	END DO
1102	END DO
1103
1104	ENDDO
1105
1106	! ------------------------------------------------------
1107	! Initialise masks and find normal/tangential directions
1108	! ------------------------------------------------------
1109
1110	! Read global 2D mask at T-points: bdytmask
1111	! -----------------------------------------
1112	! bdytmask = 1 on the computational domain AND on open boundaries
1113	! = 0 elsewhere
1114
1115	IF( ln_mask_file ) THEN
1116	CALL iom_open( cn_mask_file, inum )
1117	CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
1118	CALL iom_close( inum )
1119
1120	! Derive mask on U and V grid from mask on T grid
1121	bdyumask(:,:) = 0.e0
1122	bdyvmask(:,:) = 0.e0
1123	DO ij=1, jpjm1
1124	DO ii=1, jpim1
1125	bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
1126	bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii ,ij+1)
1127	END DO
1128	END DO
1129	CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) ! Lateral boundary cond.
1130
1131
1132	! Mask corrections
1133	! ----------------
1134	DO ik = 1, jpkm1
1135	DO ij = 1, jpj
1136	DO ii = 1, jpi
1137	tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
1138	umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
1139	vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
1140	bmask(ii,ij) = bmask(ii,ij) * bdytmask(ii,ij)
1141	END DO
1142	END DO
1143	END DO
1144
1145	DO ik = 1, jpkm1
1146	DO ij = 2, jpjm1
1147	DO ii = 2, jpim1
1148	fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij ) * bdytmask(ii+1,ij ) &
1149	& * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
1150	END DO
1151	END DO
1152	END DO
1153
1154	tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)
1155
1156	ENDIF ! ln_mask_file=.TRUE.
1157
1158	bdytmask(:,:) = tmask(:,:,1)
1159	IF( .not. ln_mask_file ) THEN
1160	! If .not. ln_mask_file then we need to derive mask on U and V grid
1161	! from mask on T grid here.
1162	bdyumask(:,:) = 0.e0
1163	bdyvmask(:,:) = 0.e0
1164	DO ij=1, jpjm1
1165	DO ii=1, jpim1
1166	bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
1167	bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii ,ij+1)
1168	END DO
1169	END DO
1170	CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) ! Lateral boundary cond.
1171	ENDIF
1172
1173	! bdy masks and bmask are now set to zero on boundary points:
1174	igrd = 1 ! In the free surface case, bmask is at T-points
1175	DO ib_bdy = 1, nb_bdy
1176	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1177	bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1178	ENDDO
1179	ENDDO
1180	!
1181	igrd = 1
1182	DO ib_bdy = 1, nb_bdy
1183	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1184	bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1185	ENDDO
1186	ENDDO
1187	!
1188	igrd = 2
1189	DO ib_bdy = 1, nb_bdy
1190	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1191	bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1192	ENDDO
1193	ENDDO
1194	!
1195	igrd = 3
1196	DO ib_bdy = 1, nb_bdy
1197	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1198	bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
1199	ENDDO
1200	ENDDO
1201
1202	! For the flagu/flagv calculation below we require a version of fmask without
1203	! the land boundary condition (shlat) included:
1204	CALL wrk_alloc(jpi,jpj,zfmask)
1205	DO ij = 2, jpjm1
1206	DO ii = 2, jpim1
1207	zfmask(ii,ij) = tmask(ii,ij ,1) * tmask(ii+1,ij ,1) &
1208	& * tmask(ii,ij+1,1) * tmask(ii+1,ij+1,1)
1209	END DO
1210	END DO
1211
1212	! Lateral boundary conditions
1213	CALL lbc_lnk( zfmask , 'F', 1. )
1214	CALL lbc_lnk( fmask , 'F', 1. ) ; CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
1215	CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
1216
1217	DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components
1218
1219	idx_bdy(ib_bdy)%flagu(:,:) = 0.e0
1220	idx_bdy(ib_bdy)%flagv(:,:) = 0.e0
1221	icount = 0
1222
1223	! Calculate relationship of U direction to the local orientation of the boundary
1224	! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
1225	! flagu = 0 : u is tangential
1226	! flagu = 1 : u is normal to the boundary and is direction is inward
1227
1228	DO igrd = 1,jpbgrd
1229	SELECT CASE( igrd )
1230	CASE( 1 )
1231	pmask => umask(:,:,1)
1232	i_offset = 0
1233	CASE( 2 )
1234	pmask => bdytmask
1235	i_offset = 1
1236	CASE( 3 )
1237	pmask => zfmask(:,:)
1238	i_offset = 0
1239	END SELECT
1240	icount = 0
1241	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1242	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1243	nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1244	zefl = pmask(nbi+i_offset-1,nbj)
1245	zwfl = pmask(nbi+i_offset,nbj)
1246	! This error check only works if you are using the bdyXmask arrays
1247	IF( i_offset == 1 .and. zefl + zwfl == 2 ) THEN
1248	icount = icount + 1
1249	IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj)
1250	ELSE
1251	idx_bdy(ib_bdy)%flagu(ib,igrd) = -zefl + zwfl
1252	ENDIF
1253	END DO
1254	IF( icount /= 0 ) THEN
1255	IF(lwp) WRITE(numout,*)
1256	IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,', &
1257	' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1258	IF(lwp) WRITE(numout,*) ' ========== '
1259	IF(lwp) WRITE(numout,*)
1260	nstop = nstop + 1
1261	ENDIF
1262	END DO
1263
1264	! Calculate relationship of V direction to the local orientation of the boundary
1265	! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
1266	! flagv = 0 : v is tangential
1267	! flagv = 1 : v is normal to the boundary and is direction is inward
1268
1269	DO igrd = 1,jpbgrd
1270	SELECT CASE( igrd )
1271	CASE( 1 )
1272	pmask => vmask(:,:,1)
1273	j_offset = 0
1274	CASE( 2 )
1275	pmask => zfmask(:,:)
1276	j_offset = 0
1277	CASE( 3 )
1278	pmask => bdytmask
1279	j_offset = 1
1280	END SELECT
1281	icount = 0
1282	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1283	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1284	nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1285	znfl = pmask(nbi,nbj+j_offset-1 )
1286	zsfl = pmask(nbi,nbj+j_offset)
1287	! This error check only works if you are using the bdyXmask arrays
1288	IF( j_offset == 1 .and. znfl + zsfl == 2 ) THEN
1289	IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj)
1290	icount = icount + 1
1291	ELSE
1292	idx_bdy(ib_bdy)%flagv(ib,igrd) = -znfl + zsfl
1293	END IF
1294	END DO
1295	IF( icount /= 0 ) THEN
1296	IF(lwp) WRITE(numout,*)
1297	IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,', &
1298	' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1299	IF(lwp) WRITE(numout,*) ' ========== '
1300	IF(lwp) WRITE(numout,*)
1301	nstop = nstop + 1
1302	ENDIF
1303	END DO
1304
1305	END DO
1306
1307	! Compute total lateral surface for volume correction:
1308	! ----------------------------------------------------
1309	! JC: this must be done at each time step with key_vvl
1310	bdysurftot = 0.e0
1311	IF( ln_vol ) THEN
1312	igrd = 2 ! Lateral surface at U-points
1313	DO ib_bdy = 1, nb_bdy
1314	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1315	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1316	nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1317	flagu => idx_bdy(ib_bdy)%flagu(ib,igrd)
1318	bdysurftot = bdysurftot + hu (nbi , nbj) &
1319	& * e2u (nbi , nbj) * ABS( flagu ) &
1320	& * tmask_i(nbi , nbj) &
1321	& * tmask_i(nbi+1, nbj)
1322	ENDDO
1323	ENDDO
1324
1325	igrd=3 ! Add lateral surface at V-points
1326	DO ib_bdy = 1, nb_bdy
1327	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1328	nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1329	nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1330	flagv => idx_bdy(ib_bdy)%flagv(ib,igrd)
1331	bdysurftot = bdysurftot + hv (nbi, nbj ) &
1332	& * e1v (nbi, nbj ) * ABS( flagv ) &
1333	& * tmask_i(nbi, nbj ) &
1334	& * tmask_i(nbi, nbj+1)
1335	ENDDO
1336	ENDDO
1337	!
1338	IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain
1339	END IF
1340	!
1341	! Tidy up
1342	!--------
1343	IF (nb_bdy>0) THEN
1344	DEALLOCATE(nbidta, nbjdta, nbrdta)
1345	ENDIF
1346
1347	CALL wrk_dealloc(jpi,jpj,zfmask)
1348
1349	IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
1350
1351	END SUBROUTINE bdy_init
1352
1353	SUBROUTINE bdy_ctl_seg
1354	!!----------------------------------------------------------------------
1355	!! * ROUTINE bdy_ctl_seg *
1356	!!
1357	!! ** Purpose : Check straight open boundary segments location
1358	!!
1359	!! ** Method : - Look for open boundary corners
1360	!! - Check that segments start or end on land
1361	!!----------------------------------------------------------------------
1362	INTEGER :: ib, ib1, ib2, ji ,jj, itest
1363	INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns
1364	REAL(wp), DIMENSION(2) :: ztestmask
1365	!!----------------------------------------------------------------------
1366	!
1367	IF (lwp) WRITE(numout,*) ' '
1368	IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1369	IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1370	!
1371	IF(lwp) WRITE(numout,*) 'Number of east segments : ', nbdysege
1372	IF(lwp) WRITE(numout,*) 'Number of west segments : ', nbdysegw
1373	IF(lwp) WRITE(numout,*) 'Number of north segments : ', nbdysegn
1374	IF(lwp) WRITE(numout,*) 'Number of south segments : ', nbdysegs
1375	! 1. Check bounds
1376	!----------------
1377	DO ib = 1, nbdysegn
1378	IF (lwp) WRITE(numout,) '*check north seg bounds pckg: ', npckgn(ib)
1379	IF ((jpjnob(ib).ge.jpjglo-1).or.&
1380	&(jpjnob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1381	IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1382	IF (jpindt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
1383	IF (jpinft(ib).ge.jpiglo) CALL ctl_stop( 'End index out of domain' )
1384	END DO
1385	!
1386	DO ib = 1, nbdysegs
1387	IF (lwp) WRITE(numout,) '*check south seg bounds pckg: ', npckgs(ib)
1388	IF ((jpjsob(ib).ge.jpjglo-1).or.&
1389	&(jpjsob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1390	IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1391	IF (jpisdt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
1392	IF (jpisft(ib).ge.jpiglo) CALL ctl_stop( 'End index out of domain' )
1393	END DO
1394	!
1395	DO ib = 1, nbdysege
1396	IF (lwp) WRITE(numout,) '*check east seg bounds pckg: ', npckge(ib)
1397	IF ((jpieob(ib).ge.jpiglo-1).or.&
1398	&(jpieob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1399	IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1400	IF (jpjedt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
1401	IF (jpjeft(ib).ge.jpjglo) CALL ctl_stop( 'End index out of domain' )
1402	END DO
1403	!
1404	DO ib = 1, nbdysegw
1405	IF (lwp) WRITE(numout,) '*check west seg bounds pckg: ', npckgw(ib)
1406	IF ((jpiwob(ib).ge.jpiglo-1).or.&
1407	&(jpiwob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1408	IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1409	IF (jpjwdt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
1410	IF (jpjwft(ib).ge.jpjglo) CALL ctl_stop( 'End index out of domain' )
1411	ENDDO
1412	!
1413	!
1414	! 2. Look for segment crossings
1415	!------------------------------
1416	IF (lwp) WRITE(numout,) '*Look for segments corners :'
1417	!
1418	itest = 0 ! corner number
1419	!
1420	! flag to detect if start or end of open boundary belongs to a corner
1421	! if not (=0), it must be on land.
1422	! if a corner is detected, save bdy package number for further tests
1423	icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1424	! South/West crossings
1425	IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1426	DO ib1 = 1, nbdysegw
1427	DO ib2 = 1, nbdysegs
1428	IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1429	& ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1430	& ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1431	& ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1432	IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN
1433	! We have a possible South-West corner
1434	! WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1435	! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1436	icornw(ib1,1) = npckgs(ib2)
1437	icorns(ib2,1) = npckgw(ib1)
1438	ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1439	IF(lwp) WRITE(numout,*)
1440	IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1441	& jpisft(ib2), jpjwft(ib1)
1442	IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
1443	IF(lwp) WRITE(numout,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1444	& ' and South segment: ',npckgs(ib2)
1445	IF(lwp) WRITE(numout,*)
1446	nstop = nstop + 1
1447	ELSE
1448	IF(lwp) WRITE(numout,*)
1449	IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices'
1450	IF(lwp) WRITE(numout,*) ' ========== Crossing problem with West segment: ',npckgw(ib1) , &
1451	& ' and South segment: ',npckgs(ib2)
1452	IF(lwp) WRITE(numout,*)
1453	nstop = nstop+1
1454	END IF
1455	END IF
1456	END DO
1457	END DO
1458	END IF
1459	!
1460	! South/East crossings
1461	IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1462	DO ib1 = 1, nbdysege
1463	DO ib2 = 1, nbdysegs
1464	IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1465	& ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1466	& ( jpjedt(ib1)<=jpjsob(ib2) ).AND. &
1467	& ( jpjeft(ib1)>=jpjsob(ib2) )) THEN
1468	IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1469	! We have a possible South-East corner
1470	! WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1471	! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1472	icorne(ib1,1) = npckgs(ib2)
1473	icorns(ib2,2) = npckge(ib1)
1474	ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1475	IF(lwp) WRITE(numout,*)
1476	IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1477	& jpisdt(ib2), jpjeft(ib1)
1478	IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
1479	IF(lwp) WRITE(numout,*) ' Crossing problem with East segment: ',npckge(ib1), &
1480	& ' and South segment: ',npckgs(ib2)
1481	IF(lwp) WRITE(numout,*)
1482	nstop = nstop + 1
1483	ELSE
1484	IF(lwp) WRITE(numout,*)
1485	IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices'
1486	IF(lwp) WRITE(numout,*) ' ========== Crossing problem with East segment: ',npckge(ib1), &
1487	& ' and South segment: ',npckgs(ib2)
1488	IF(lwp) WRITE(numout,*)
1489	nstop = nstop + 1
1490	END IF
1491	END IF
1492	END DO
1493	END DO
1494	END IF
1495	!
1496	! North/West crossings
1497	IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1498	DO ib1 = 1, nbdysegw
1499	DO ib2 = 1, nbdysegn
1500	IF (( jpindt(ib2)<=jpiwob(ib1) ).AND. &
1501	& ( jpinft(ib2)>=jpiwob(ib1) ).AND. &
1502	& ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1503	& ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1504	IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1505	! We have a possible North-West corner
1506	! WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1507	! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1508	icornw(ib1,2) = npckgn(ib2)
1509	icornn(ib2,1) = npckgw(ib1)
1510	ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1511	IF(lwp) WRITE(numout,*)
1512	IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1513	& jpinft(ib2), jpjwdt(ib1)
1514	IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
1515	IF(lwp) WRITE(numout,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1516	& ' and North segment: ',npckgn(ib2)
1517	IF(lwp) WRITE(numout,*)
1518	nstop = nstop + 1
1519	ELSE
1520	IF(lwp) WRITE(numout,*)
1521	IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices'
1522	IF(lwp) WRITE(numout,*) ' ========== Crossing problem with West segment: ',npckgw(ib1), &
1523	& ' and North segment: ',npckgn(ib2)
1524	IF(lwp) WRITE(numout,*)
1525	nstop = nstop + 1
1526	END IF
1527	END IF
1528	END DO
1529	END DO
1530	END IF
1531	!
1532	! North/East crossings
1533	IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1534	DO ib1 = 1, nbdysege
1535	DO ib2 = 1, nbdysegn
1536	IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1537	& ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1538	& ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1539	& ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1540	IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1541	! We have a possible North-East corner
1542	! WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1543	! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1544	icorne(ib1,2) = npckgn(ib2)
1545	icornn(ib2,2) = npckge(ib1)
1546	ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1547	IF(lwp) WRITE(numout,*)
1548	IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1549	& jpindt(ib2), jpjedt(ib1)
1550	IF(lwp) WRITE(numout,*) ' ========== Not allowed yet'
1551	IF(lwp) WRITE(numout,*) ' Crossing problem with East segment: ',npckge(ib1), &
1552	& ' and North segment: ',npckgn(ib2)
1553	IF(lwp) WRITE(numout,*)
1554	nstop = nstop + 1
1555	ELSE
1556	IF(lwp) WRITE(numout,*)
1557	IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices'
1558	IF(lwp) WRITE(numout,*) ' ========== Crossing problem with East segment: ',npckge(ib1), &
1559	& ' and North segment: ',npckgn(ib2)
1560	IF(lwp) WRITE(numout,*)
1561	nstop = nstop + 1
1562	END IF
1563	END IF
1564	END DO
1565	END DO
1566	END IF
1567	!
1568	! 3. Check if segment extremities are on land
1569	!--------------------------------------------
1570	!
1571	! West segments
1572	DO ib = 1, nbdysegw
1573	! get mask at boundary extremities:
1574	ztestmask(1:2)=0.
1575	DO ji = 1, jpi
1576	DO jj = 1, jpj
1577	IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &
1578	& ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1579	IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &
1580	& ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1)
1581	END DO
1582	END DO
1583	IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
1584
1585	IF (ztestmask(1)==1) THEN
1586	IF (icornw(ib,1)==0) THEN
1587	IF(lwp) WRITE(numout,*)
1588	IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
1589	IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner'
1590	IF(lwp) WRITE(numout,*)
1591	nstop = nstop + 1
1592	ELSE
1593	! This is a corner
1594	WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1595	CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1596	itest=itest+1
1597	ENDIF
1598	ENDIF
1599	IF (ztestmask(2)==1) THEN
1600	IF (icornw(ib,2)==0) THEN
1601	IF(lwp) WRITE(numout,*)
1602	IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
1603	IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner'
1604	IF(lwp) WRITE(numout,*)
1605	nstop = nstop + 1
1606	ELSE
1607	! This is a corner
1608	WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1609	CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1610	itest=itest+1
1611	ENDIF
1612	ENDIF
1613	END DO
1614	!
1615	! East segments
1616	DO ib = 1, nbdysege
1617	! get mask at boundary extremities:
1618	ztestmask(1:2)=0.
1619	DO ji = 1, jpi
1620	DO jj = 1, jpj
1621	IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &
1622	& ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
1623	IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &
1624	& ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1)
1625	END DO
1626	END DO
1627	IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
1628
1629	IF (ztestmask(1)==1) THEN
1630	IF (icorne(ib,1)==0) THEN
1631	IF(lwp) WRITE(numout,*)
1632	IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
1633	IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner'
1634	IF(lwp) WRITE(numout,*)
1635	nstop = nstop + 1
1636	ELSE
1637	! This is a corner
1638	WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1639	CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1640	itest=itest+1
1641	ENDIF
1642	ENDIF
1643	IF (ztestmask(2)==1) THEN
1644	IF (icorne(ib,2)==0) THEN
1645	IF(lwp) WRITE(numout,*)
1646	IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
1647	IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner'
1648	IF(lwp) WRITE(numout,*)
1649	nstop = nstop + 1
1650	ELSE
1651	! This is a corner
1652	WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1653	CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1654	itest=itest+1
1655	ENDIF
1656	ENDIF
1657	END DO
1658	!
1659	! South segments
1660	DO ib = 1, nbdysegs
1661	! get mask at boundary extremities:
1662	ztestmask(1:2)=0.
1663	DO ji = 1, jpi
1664	DO jj = 1, jpj
1665	IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &
1666	& ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1667	IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &
1668	& ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1)
1669	END DO
1670	END DO
1671	IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
1672
1673	IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1674	IF(lwp) WRITE(numout,*)
1675	IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
1676	IF(lwp) WRITE(numout,*) ' ========== does not start on land or on a corner'
1677	IF(lwp) WRITE(numout,*)
1678	nstop = nstop + 1
1679	ENDIF
1680	IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1681	IF(lwp) WRITE(numout,*)
1682	IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
1683	IF(lwp) WRITE(numout,*) ' ========== does not end on land or on a corner'
1684	IF(lwp) WRITE(numout,*)
1685	nstop = nstop + 1
1686	ENDIF
1687	END DO
1688	!
1689	! North segments
1690	DO ib = 1, nbdysegn
1691	! get mask at boundary extremities:
1692	ztestmask(1:2)=0.
1693	DO ji = 1, jpi
1694	DO jj = 1, jpj
1695	IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &
1696	& ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
1697	IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &
1698	& ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1)
1699	END DO
1700	END DO
1701	IF( lk_mpp ) CALL mpp_sum( ztestmask, 2 ) ! sum over the global domain
1702
1703	IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1704	IF(lwp) WRITE(numout,*)
1705	IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
1706	IF(lwp) WRITE(numout,*) ' ========== does not start on land'
1707	IF(lwp) WRITE(numout,*)
1708	nstop = nstop + 1
1709	ENDIF
1710	IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1711	IF(lwp) WRITE(numout,*)
1712	IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
1713	IF(lwp) WRITE(numout,*) ' ========== does not end on land'
1714	IF(lwp) WRITE(numout,*)
1715	nstop = nstop + 1
1716	ENDIF
1717	END DO
1718	!
1719	IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1720	!
1721	! Other tests TBD:
1722	! segments completly on land
1723	! optimized open boundary array length according to landmask
1724	! Nudging layers that overlap with interior domain
1725	!
1726	END SUBROUTINE bdy_ctl_seg
1727
1728	SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1729	!!----------------------------------------------------------------------
1730	!! * ROUTINE bdy_ctl_corn *
1731	!!
1732	!! ** Purpose : Check numerical schemes consistency between
1733	!! segments having a common corner
1734	!!
1735	!! ** Method :
1736	!!----------------------------------------------------------------------
1737	INTEGER, INTENT(in) :: ib1, ib2
1738	INTEGER :: itest
1739	!!----------------------------------------------------------------------
1740	itest = 0
1741
1742	IF (cn_dyn2d(ib1)/=cn_dyn2d(ib2)) itest = itest + 1
1743	IF (cn_dyn3d(ib1)/=cn_dyn3d(ib2)) itest = itest + 1
1744	IF (cn_tra(ib1)/=cn_tra(ib2)) itest = itest + 1
1745	!
1746	IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1
1747	IF (nn_dyn3d_dta(ib1)/=nn_dyn3d_dta(ib2)) itest = itest + 1
1748	IF (nn_tra_dta(ib1)/=nn_tra_dta(ib2)) itest = itest + 1
1749	!
1750	IF (nn_rimwidth(ib1)/=nn_rimwidth(ib2)) itest = itest + 1
1751	!
1752	IF ( itest>0 ) THEN
1753	IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2
1754	IF(lwp) WRITE(numout,*) ' ========== have different open bdy schemes'
1755	IF(lwp) WRITE(numout,*)
1756	nstop = nstop + 1
1757	ENDIF
1758	!
1759	END SUBROUTINE bdy_ctl_corn
1760
1761	#else
1762	!!---------------------------------------------------------------------------------
1763	!! Dummy module NO open boundaries
1764	!!---------------------------------------------------------------------------------
1765	CONTAINS
1766	SUBROUTINE bdy_init ! Dummy routine
1767	END SUBROUTINE bdy_init
1768	#endif
1769
1770	!!=================================================================================
1771	END MODULE bdyini

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

Download in other formats: