Context Navigation

bdyini.F90 @ 7567

Last change on this file since 7567 was 7567, checked in by hadjt, 7 years ago
CO6 version adapted for shelf seas climate projections, including added diagnostics
File size: 86.2 KB

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

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

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

Context Navigation

source: branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 7567

Download in other formats: