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

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 8790

Last change on this file since 8790 was 8790, checked in by deazer, 7 years ago
Add optional ssh switch to del with Baltic in CO7 Also added lbclnk in dynkeg for reproducible runs.
File size: 85.3 KB

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

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

Download in other formats: