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

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

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

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 4332

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

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

Download in other formats: