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 NEMO/trunk/src/OCE/BDY – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/BDY/bdyini.F90 @ 15363

Last change on this file since 15363 was 15363, checked in by smasson, 2 years ago
trunk: continuing with ticket #2731
Property svn:keywords set to `Id`
File size: 107.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) optimization of BDY communications
15	!! 3.7 ! 2016 (T. Lovato) Remove bdy macro, call here init for dta and tides
16	!!----------------------------------------------------------------------
17	!! bdy_init : Initialization of unstructured open boundaries
18	!!----------------------------------------------------------------------
19	USE oce ! ocean dynamics and tracers variables
20	USE dom_oce ! ocean space and time domain
21	USE sbc_oce , ONLY: nn_ice
22	USE bdy_oce ! unstructured open boundary conditions
23	USE bdydta ! open boundary cond. setting (bdy_dta_init routine)
24	USE bdytides ! open boundary cond. setting (bdytide_init routine)
25	USE tide_mod, ONLY: ln_tide ! tidal forcing
26	USE phycst , ONLY: rday
27	!
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
33	IMPLICIT NONE
34	PRIVATE
35
36	PUBLIC bdy_init ! routine called in nemo_init
37	PUBLIC find_neib ! routine called in bdy_nmn
38
39	INTEGER, PARAMETER :: jp_nseg = 100 !
40	! Straight open boundary segment parameters:
41	INTEGER :: nbdysege, nbdysegw, nbdysegn, nbdysegs
42	INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge !
43	INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw !
44	INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn !
45	INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs !
46
47	!! * Substitutions
48	# include "do_loop_substitute.h90"
49	!!----------------------------------------------------------------------
50	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
51	!! $Id$
52	!! Software governed by the CeCILL license (see ./LICENSE)
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/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file, &
69	& ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, &
70	& cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, &
71	& ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
72	& cn_ice, nn_ice_dta, &
73	& ln_vol, nn_volctl, nn_rimwidth
74	!
75	INTEGER :: ios ! Local integer output status for namelist read
76	!!----------------------------------------------------------------------
77
78	! ------------------------
79	! Read namelist parameters
80	! ------------------------
81	READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901)
82	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist' )
83	! make sur that all elements of the namelist variables have a default definition from namelist_ref
84	ln_coords_file (2:jp_bdy) = ln_coords_file (1)
85	cn_coords_file (2:jp_bdy) = cn_coords_file (1)
86	cn_dyn2d (2:jp_bdy) = cn_dyn2d (1)
87	nn_dyn2d_dta (2:jp_bdy) = nn_dyn2d_dta (1)
88	cn_dyn3d (2:jp_bdy) = cn_dyn3d (1)
89	nn_dyn3d_dta (2:jp_bdy) = nn_dyn3d_dta (1)
90	cn_tra (2:jp_bdy) = cn_tra (1)
91	nn_tra_dta (2:jp_bdy) = nn_tra_dta (1)
92	ln_tra_dmp (2:jp_bdy) = ln_tra_dmp (1)
93	ln_dyn3d_dmp (2:jp_bdy) = ln_dyn3d_dmp (1)
94	rn_time_dmp (2:jp_bdy) = rn_time_dmp (1)
95	rn_time_dmp_out(2:jp_bdy) = rn_time_dmp_out(1)
96	cn_ice (2:jp_bdy) = cn_ice (1)
97	nn_ice_dta (2:jp_bdy) = nn_ice_dta (1)
98	READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
99	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist' )
100	IF(lwm) WRITE ( numond, nambdy )
101
102	IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE. ! forced for Agrif children
103
104	IF( nb_bdy == 0 ) ln_bdy = .FALSE.
105
106	! -----------------------------------------
107	! unstructured open boundaries use control
108	! -----------------------------------------
109	IF ( ln_bdy ) THEN
110	IF(lwp) WRITE(numout,*)
111	IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
112	IF(lwp) WRITE(numout,*) '~~~~~~~~'
113	!
114	! Open boundaries definition (arrays and masks)
115	CALL bdy_def
116	IF( ln_meshmask ) CALL bdy_meshwri()
117	!
118	! Open boundaries initialisation of external data arrays
119	CALL bdy_dta_init
120	!
121	! Open boundaries initialisation of tidal harmonic forcing
122	IF( ln_tide ) CALL bdytide_init
123	!
124	ELSE
125	IF(lwp) WRITE(numout,*)
126	IF(lwp) WRITE(numout,*) 'bdy_init : open boundaries not used (ln_bdy = F)'
127	IF(lwp) WRITE(numout,*) '~~~~~~~~'
128	!
129	ENDIF
130	!
131	END SUBROUTINE bdy_init
132
133
134	SUBROUTINE bdy_def
135	!!----------------------------------------------------------------------
136	!! * ROUTINE bdy_init *
137	!!
138	!! ** Purpose : Definition of unstructured open boundaries.
139	!!
140	!! ** Method : Read initialization arrays (mask, indices) to identify
141	!! an unstructured open boundary
142	!!
143	!! ** Input : bdy_init.nc, input file for unstructured open boundaries
144	!!----------------------------------------------------------------------
145	INTEGER :: ji, jj ! dummy loop indices
146	INTEGER :: ib_bdy, ii, ij, igrd, ib, ir, iseg ! dummy loop indices
147	INTEGER :: icount, icountr, icountr0, ibr_max ! local integers
148	INTEGER :: ilen1 ! - -
149	INTEGER :: iiRst, iiRnd, iiSst, iiSnd, iiRcorn, iiSdiag, iiSsono
150	INTEGER :: ijRst, ijRnd, ijSst, ijSnd, ijRcorn, ijSdiag, ijSsono
151	INTEGER :: iiout, ijout, iioutdir, ijoutdir, icnt
152	INTEGER :: iRnei, iRdiag, iRsono
153	INTEGER :: iSnei, iSdiag, iSsono ! - -
154	INTEGER :: iwe, ies, iso, ino, inum, id_dummy ! - -
155	INTEGER :: jpbdta ! - -
156	INTEGER :: ib_bdy1, ib_bdy2, ib1, ib2 ! - -
157	INTEGER :: ii1, ii2, ii3, ij1, ij2, ij3 ! - -
158	INTEGER :: iibe, ijbe, iibi, ijbi ! - -
159	INTEGER :: flagu, flagv ! short cuts
160	INTEGER :: nbdyind, nbdybeg, nbdyend
161	INTEGER , DIMENSION(4) :: kdimsz
162	INTEGER , DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays
163	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta
164	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points
165	CHARACTER(LEN=1) , DIMENSION(jpbgrd) :: cgrid
166	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zz_read ! work space for 2D global boundary data
167	REAL(wp), POINTER , DIMENSION(:,:) :: zmask ! pointer to 2D mask fields
168	REAL(wp) , DIMENSION(jpi,jpj) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat)
169	REAL(wp) , DIMENSION(jpi,jpj) :: ztmask, zumask, zvmask ! temporary u/v mask array
170	REAL(wp) , DIMENSION(jpi,jpj) :: zzbdy
171	!!----------------------------------------------------------------------
172	!
173	cgrid = (/'t','u','v'/)
174
175	! -----------------------------------------
176	! Check and write out namelist parameters
177	! -----------------------------------------
178
179	IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy
180
181	DO ib_bdy = 1,nb_bdy
182
183	IF(lwp) THEN
184	WRITE(numout,*) ' '
185	WRITE(numout,*) '------ Open boundary data set ',ib_bdy,' ------'
186	IF( ln_coords_file(ib_bdy) ) THEN
187	WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
188	ELSE
189	WRITE(numout,*) 'Boundary defined in namelist.'
190	ENDIF
191	WRITE(numout,*)
192	ENDIF
193
194	! barotropic bdy
195	!----------------
196	IF(lwp) THEN
197	WRITE(numout,*) 'Boundary conditions for barotropic solution: '
198	SELECT CASE( cn_dyn2d(ib_bdy) )
199	CASE( 'none' ) ; WRITE(numout,*) ' no open boundary condition'
200	CASE( 'frs' ) ; WRITE(numout,*) ' Flow Relaxation Scheme'
201	CASE( 'flather' ) ; WRITE(numout,*) ' Flather radiation condition'
202	CASE( 'orlanski' ) ; WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'
203	CASE( 'orlanski_npo' ) ; WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'
204	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn2d' )
205	END SELECT
206	ENDIF
207
208	dta_bdy(ib_bdy)%lneed_ssh = cn_dyn2d(ib_bdy) == 'flather'
209	dta_bdy(ib_bdy)%lneed_dyn2d = cn_dyn2d(ib_bdy) /= 'none'
210
211	IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
212	SELECT CASE( nn_dyn2d_dta(ib_bdy) ) !
213	CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data'
214	CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file'
215	CASE( 2 ) ; WRITE(numout,*) ' tidal harmonic forcing taken from file'
216	CASE( 3 ) ; WRITE(numout,*) ' boundary data AND tidal harmonic forcing taken from files'
217	CASE DEFAULT ; CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
218	END SELECT
219	ENDIF
220	IF ( dta_bdy(ib_bdy)%lneed_dyn2d .AND. nn_dyn2d_dta(ib_bdy) .GE. 2 .AND. .NOT.ln_tide ) THEN
221	CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' )
222	ENDIF
223	IF(lwp) WRITE(numout,*)
224
225	! baroclinic bdy
226	!----------------
227	IF(lwp) THEN
228	WRITE(numout,*) 'Boundary conditions for baroclinic velocities: '
229	SELECT CASE( cn_dyn3d(ib_bdy) )
230	CASE('none') ; WRITE(numout,*) ' no open boundary condition'
231	CASE('frs') ; WRITE(numout,*) ' Flow Relaxation Scheme'
232	CASE('specified') ; WRITE(numout,*) ' Specified value'
233	CASE('neumann') ; WRITE(numout,*) ' Neumann conditions'
234	CASE('zerograd') ; WRITE(numout,*) ' Zero gradient for baroclinic velocities'
235	CASE('zero') ; WRITE(numout,*) ' Zero baroclinic velocities (runoff case)'
236	CASE('orlanski') ; WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'
237	CASE('orlanski_npo') ; WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'
238	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn3d' )
239	END SELECT
240	ENDIF
241
242	dta_bdy(ib_bdy)%lneed_dyn3d = cn_dyn3d(ib_bdy) == 'frs' .OR. cn_dyn3d(ib_bdy) == 'specified' &
243	& .OR. cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo'
244
245	IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn3d ) THEN
246	SELECT CASE( nn_dyn3d_dta(ib_bdy) ) !
247	CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data'
248	CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file'
249	CASE DEFAULT ; CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
250	END SELECT
251	END IF
252
253	IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
254	IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN
255	IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
256	ln_dyn3d_dmp(ib_bdy) = .false.
257	ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN
258	CALL ctl_stop( 'Use FRS OR relaxation' )
259	ELSE
260	IF(lwp) WRITE(numout,*) ' + baroclinic velocities relaxation zone'
261	IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days'
262	IF(rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
263	dta_bdy(ib_bdy)%lneed_dyn3d = .TRUE.
264	ENDIF
265	ELSE
266	IF(lwp) WRITE(numout,*) ' NO relaxation on baroclinic velocities'
267	ENDIF
268	IF(lwp) WRITE(numout,*)
269
270	! tra bdy
271	!----------------
272	IF(lwp) THEN
273	WRITE(numout,*) 'Boundary conditions for temperature and salinity: '
274	SELECT CASE( cn_tra(ib_bdy) )
275	CASE('none') ; WRITE(numout,*) ' no open boundary condition'
276	CASE('frs') ; WRITE(numout,*) ' Flow Relaxation Scheme'
277	CASE('specified') ; WRITE(numout,*) ' Specified value'
278	CASE('neumann') ; WRITE(numout,*) ' Neumann conditions'
279	CASE('runoff') ; WRITE(numout,*) ' Runoff conditions : Neumann for T and specified to 0.1 for salinity'
280	CASE('orlanski') ; WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging'
281	CASE('orlanski_npo') ; WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging'
282	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_tra' )
283	END SELECT
284	ENDIF
285
286	dta_bdy(ib_bdy)%lneed_tra = cn_tra(ib_bdy) == 'frs' .OR. cn_tra(ib_bdy) == 'specified' &
287	& .OR. cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo'
288
289	IF( lwp .AND. dta_bdy(ib_bdy)%lneed_tra ) THEN
290	SELECT CASE( nn_tra_dta(ib_bdy) ) !
291	CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data'
292	CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file'
293	CASE DEFAULT ; CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
294	END SELECT
295	ENDIF
296
297	IF ( ln_tra_dmp(ib_bdy) ) THEN
298	IF ( cn_tra(ib_bdy) == 'none' ) THEN
299	IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
300	ln_tra_dmp(ib_bdy) = .false.
301	ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN
302	CALL ctl_stop( 'Use FRS OR relaxation' )
303	ELSE
304	IF(lwp) WRITE(numout,*) ' + T/S relaxation zone'
305	IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days'
306	IF(lwp) WRITE(numout,*) ' Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'
307	IF(lwp.AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
308	dta_bdy(ib_bdy)%lneed_tra = .TRUE.
309	ENDIF
310	ELSE
311	IF(lwp) WRITE(numout,*) ' NO T/S relaxation'
312	ENDIF
313	IF(lwp) WRITE(numout,*)
314
315	#if defined key_si3
316	IF(lwp) THEN
317	WRITE(numout,*) 'Boundary conditions for sea ice: '
318	SELECT CASE( cn_ice(ib_bdy) )
319	CASE('none') ; WRITE(numout,*) ' no open boundary condition'
320	CASE('frs') ; WRITE(numout,*) ' Flow Relaxation Scheme'
321	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_ice' )
322	END SELECT
323	ENDIF
324
325	dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none'
326
327	IF( dta_bdy(ib_bdy)%lneed_ice .AND. nn_ice /= 2 ) THEN
328	WRITE(ctmp1,*) 'bdy number ', ib_bdy,', needs ice model but nn_ice = ', nn_ice
329	CALL ctl_stop( ctmp1 )
330	ENDIF
331
332	IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN
333	SELECT CASE( nn_ice_dta(ib_bdy) ) !
334	CASE( 0 ) ; WRITE(numout,*) ' initial state used for bdy data'
335	CASE( 1 ) ; WRITE(numout,*) ' boundary data taken from file'
336	CASE DEFAULT ; CALL ctl_stop( 'nn_ice_dta must be 0 or 1' )
337	END SELECT
338	ENDIF
339	#else
340	dta_bdy(ib_bdy)%lneed_ice = .FALSE.
341	#endif
342	!
343	IF(lwp) WRITE(numout,*)
344	IF(lwp) WRITE(numout,*) ' Width of relaxation zone = ', nn_rimwidth(ib_bdy)
345	IF(lwp) WRITE(numout,*)
346	!
347	END DO ! nb_bdy
348
349	IF( lwp ) THEN
350	IF( ln_vol ) THEN ! check volume conservation (nn_volctl value)
351	WRITE(numout,*) 'Volume correction applied at open boundaries'
352	WRITE(numout,*)
353	SELECT CASE ( nn_volctl )
354	CASE( 1 ) ; WRITE(numout,*) ' The total volume will be constant'
355	CASE( 0 ) ; WRITE(numout,*) ' The total volume will vary according to the surface E-P flux'
356	CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' )
357	END SELECT
358	WRITE(numout,*)
359	!
360	! sanity check if used with tides
361	IF( ln_tide ) THEN
362	WRITE(numout,*) ' The total volume correction is not working with tides. '
363	WRITE(numout,*) ' Set ln_vol to .FALSE. '
364	WRITE(numout,*) ' or '
365	WRITE(numout,*) ' equilibriate your bdy input files '
366	CALL ctl_stop( 'The total volume correction is not working with tides.' )
367	END IF
368	ELSE
369	WRITE(numout,*) 'No volume correction applied at open boundaries'
370	WRITE(numout,*)
371	ENDIF
372	ENDIF
373
374	! -------------------------------------------------
375	! Initialise indices arrays for open boundaries
376	! -------------------------------------------------
377
378	nblendta(:,:) = 0
379	nbdysege = 0
380	nbdysegw = 0
381	nbdysegn = 0
382	nbdysegs = 0
383
384	! Define all boundaries
385	! ---------------------
386	DO ib_bdy = 1, nb_bdy
387	!
388	IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! build bdy coordinates with segments defined in namelist
389
390	CALL bdy_read_seg( ib_bdy, nblendta(:,ib_bdy) )
391
392	ELSE ! Read size of arrays in boundary coordinates file.
393
394	CALL iom_open( cn_coords_file(ib_bdy), inum )
395	DO igrd = 1, jpbgrd
396	id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )
397	nblendta(igrd,ib_bdy) = MAXVAL(kdimsz)
398	END DO
399	CALL iom_close( inum )
400	ENDIF
401	!
402	END DO ! ib_bdy
403
404	! Now look for crossings in user (namelist) defined open boundary segments:
405	IF( nbdysege > 0 .OR. nbdysegw > 0 .OR. nbdysegn > 0 .OR. nbdysegs > 0) CALL bdy_ctl_seg
406
407	! Allocate arrays
408	!---------------
409	jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
410	ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), nbrdta(jpbdta, jpbgrd, nb_bdy) )
411	nbrdta(:,:,:) = 0 ! initialize nbrdta as it may not be completely defined for each bdy
412
413	! Calculate global boundary index arrays or read in from file
414	!------------------------------------------------------------
415	! 1. Read global index arrays from boundary coordinates file.
416	DO ib_bdy = 1, nb_bdy
417	!
418	IF( ln_coords_file(ib_bdy) ) THEN
419	!
420	ALLOCATE( zz_read( MAXVAL(nblendta), 1 ) )
421	CALL iom_open( cn_coords_file(ib_bdy), inum )
422	!
423	DO igrd = 1, jpbgrd
424	CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
425	DO ii = 1,nblendta(igrd,ib_bdy)
426	nbidta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) + nn_hls
427	END DO
428	CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
429	DO ii = 1,nblendta(igrd,ib_bdy)
430	nbjdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) + nn_hls
431	END DO
432	CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
433	DO ii = 1,nblendta(igrd,ib_bdy)
434	nbrdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
435	END DO
436	!
437	ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
438	IF(lwp) WRITE(numout,*)
439	IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
440	IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
441	IF (ibr_max < nn_rimwidth(ib_bdy)) &
442	CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
443	END DO
444	!
445	CALL iom_close( inum )
446	DEALLOCATE( zz_read )
447	!
448	ENDIF
449	!
450	END DO
451
452	! 2. Now fill indices corresponding to straight open boundary arrays:
453	CALL bdy_coords_seg( nbidta, nbjdta, nbrdta )
454
455	! Deal with duplicated points
456	!-----------------------------
457	! We assign negative indices to duplicated points (to remove them from bdy points to be updated)
458	! if their distance to the bdy is greater than the other
459	! If their distance are the same, just keep only one to avoid updating a point twice
460	DO igrd = 1, jpbgrd
461	DO ib_bdy1 = 1, nb_bdy
462	DO ib_bdy2 = 1, nb_bdy
463	IF (ib_bdy1/=ib_bdy2) THEN
464	DO ib1 = 1, nblendta(igrd,ib_bdy1)
465	DO ib2 = 1, nblendta(igrd,ib_bdy2)
466	IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
467	& (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
468	! IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &
469	! & nbidta(ib1, igrd, ib_bdy1), &
470	! & nbjdta(ib2, igrd, ib_bdy2)
471	! keep only points with the lowest distance to boundary:
472	IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN
473	nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
474	nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
475	ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN
476	nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
477	nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
478	! Arbitrary choice if distances are the same:
479	ELSE
480	nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
481	nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
482	ENDIF
483	END IF
484	END DO
485	END DO
486	ENDIF
487	END DO
488	END DO
489	END DO
490	!
491	! Find lenght of boundaries and rim on local mpi domain
492	!------------------------------------------------------
493	!
494	iwe = mig(1)
495	ies = mig(jpi)
496	iso = mjg(1)
497	ino = mjg(jpj)
498	!
499	DO ib_bdy = 1, nb_bdy
500	DO igrd = 1, jpbgrd
501	icount = 0 ! initialization of local bdy length
502	icountr = 0 ! initialization of local rim 0 and rim 1 bdy length
503	icountr0 = 0 ! initialization of local rim 0 bdy length
504	idx_bdy(ib_bdy)%nblen(igrd) = 0
505	idx_bdy(ib_bdy)%nblenrim(igrd) = 0
506	idx_bdy(ib_bdy)%nblenrim0(igrd) = 0
507	DO ib = 1, nblendta(igrd,ib_bdy)
508	! check that data is in correct order in file
509	IF( ib > 1 ) THEN
510	IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ib-1,igrd,ib_bdy) ) THEN
511	CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', &
512	& ' in order of distance from edge nbr A utility for re-ordering ', &
513	& ' boundary coordinates and data files exists in the TOOLS/OBC directory')
514	ENDIF
515	ENDIF
516	! check if point is in local domain
517	IF( nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND. &
518	& nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino ) THEN
519	!
520	icount = icount + 1
521	IF( nbrdta(ib,igrd,ib_bdy) == 1 .OR. nbrdta(ib,igrd,ib_bdy) == 0 ) icountr = icountr + 1
522	IF( nbrdta(ib,igrd,ib_bdy) == 0 ) icountr0 = icountr0 + 1
523	ENDIF
524	END DO
525	idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc
526	idx_bdy(ib_bdy)%nblenrim (igrd) = icountr !: length of rim 0 and rim 1 boundary data on each proc
527	idx_bdy(ib_bdy)%nblenrim0(igrd) = icountr0 !: length of rim 0 boundary data on each proc
528	END DO ! igrd
529
530	! Allocate index arrays for this boundary set
531	!--------------------------------------------
532	ilen1 = MAXVAL( idx_bdy(ib_bdy)%nblen(:) )
533	ALLOCATE( idx_bdy(ib_bdy)%nbi (ilen1,jpbgrd) , &
534	& idx_bdy(ib_bdy)%nbj (ilen1,jpbgrd) , &
535	& idx_bdy(ib_bdy)%nbr (ilen1,jpbgrd) , &
536	& idx_bdy(ib_bdy)%nbd (ilen1,jpbgrd) , &
537	& idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) , &
538	& idx_bdy(ib_bdy)%ntreat(ilen1,jpbgrd) , &
539	& idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) , &
540	& idx_bdy(ib_bdy)%nbw (ilen1,jpbgrd) , &
541	& idx_bdy(ib_bdy)%flagu (ilen1,jpbgrd) , &
542	& idx_bdy(ib_bdy)%flagv (ilen1,jpbgrd) )
543
544	! Dispatch mapping indices and discrete distances on each processor
545	! -----------------------------------------------------------------
546	DO igrd = 1, jpbgrd
547	icount = 0
548	! Outer loop on rimwidth to ensure outermost points come first in the local arrays.
549	DO ir = 0, nn_rimwidth(ib_bdy)
550	DO ib = 1, nblendta(igrd,ib_bdy)
551	! check if point is in local domain and equals ir
552	IF( nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND. &
553	& nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino .AND. &
554	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
555	!
556	icount = icount + 1
557	idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy) - mig(1) + 1 ! global to local indexes
558	idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy) - mjg(1) + 1 ! global to local indexes
559	idx_bdy(ib_bdy)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_bdy)
560	idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
561	ENDIF
562	END DO
563	END DO
564	END DO ! igrd
565
566	END DO ! ib_bdy
567
568	! Initialize array indicating communications in bdy
569	! -------------------------------------------------
570	ALLOCATE( lsend_bdyolr(nb_bdy,jpbgrd,8,0:1), lrecv_bdyolr(nb_bdy,jpbgrd,8,0:1) )
571	lsend_bdyolr(:,:,:,:) = .false.
572	lrecv_bdyolr(:,:,:,:) = .false.
573
574	DO ib_bdy = 1, nb_bdy
575	DO igrd = 1, jpbgrd
576	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! only the rim triggers communications, see bdy routines
577	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
578	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
579	IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN ; ir = 0
580	ELSE ; ir = 1
581	END IF
582	!
583	! check if point has to be sent to a neighbour
584	IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. ij >= Njs0 .AND. ij <= Nje0 ) THEN ! we inner side
585	IF( mpiSnei(nn_hls,jpwe) > -1 ) lsend_bdyolr(ib_bdy,igrd,jpwe,ir) = .TRUE.
586	ENDIF
587	IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. ij >= Njs0 .AND. ij <= Nje0 ) THEN ! ea inner side
588	IF( mpiSnei(nn_hls,jpea) > -1 ) lsend_bdyolr(ib_bdy,igrd,jpea,ir) = .TRUE.
589	ENDIF
590	IF( ii >= Nis0 .AND. ii <= Nie0 .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN ! so inner side
591	IF( mpiSnei(nn_hls,jpso) > -1 ) lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
592	ENDIF
593	IF( ii < Nis0 .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN ! so side we-halo
594	IF( mpiSnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 ) lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
595	ENDIF
596	IF( ii > Nie0 .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN ! so side ea-halo
597	IF( mpiSnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 ) lsend_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
598	ENDIF
599	IF( ii >= Nis0 .AND. ii <= Nie0 .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN ! no inner side
600	IF( mpiSnei(nn_hls,jpno) > -1 ) lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
601	ENDIF
602	IF( ii < Nis0 .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN ! no side we-halo
603	IF( mpiSnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 ) lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
604	ENDIF
605	IF( ii > Nie0 .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN ! no side ea-halo
606	IF( mpiSnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 ) lsend_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
607	ENDIF
608	IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN ! sw inner corner
609	IF( mpiSnei(nn_hls,jpsw) > -1 ) lsend_bdyolr(ib_bdy,igrd,jpsw,ir) = .TRUE.
610	ENDIF
611	IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. ij >= Njs0 .AND. ij < Njs0 + nn_hls ) THEN ! se inner corner
612	IF( mpiSnei(nn_hls,jpse) > -1 ) lsend_bdyolr(ib_bdy,igrd,jpse,ir) = .TRUE.
613	ENDIF
614	IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN ! nw inner corner
615	IF( mpiSnei(nn_hls,jpnw) > -1 ) lsend_bdyolr(ib_bdy,igrd,jpnw,ir) = .TRUE.
616	ENDIF
617	IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. ij <= Nje0 .AND. ij > Nje0 - nn_hls ) THEN ! ne inner corner
618	IF( mpiSnei(nn_hls,jpne) > -1 ) lsend_bdyolr(ib_bdy,igrd,jpne,ir) = .TRUE.
619	ENDIF
620	!
621	! check if point has to be received from a neighbour
622	IF( ii < Nis0 .AND. ij >= Njs0 .AND. ij <= Nje0 ) THEN ! we side
623	IF( mpiRnei(nn_hls,jpwe) > -1 ) lrecv_bdyolr(ib_bdy,igrd,jpwe,ir) = .TRUE.
624	ENDIF
625	IF( ii > Nie0 .AND. ij >= Njs0 .AND. ij <= Nje0 ) THEN ! ea side
626	IF( mpiRnei(nn_hls,jpea) > -1 ) lrecv_bdyolr(ib_bdy,igrd,jpea,ir) = .TRUE.
627	ENDIF
628	IF( ii >= Nis0 .AND. ii <= Nie0 .AND. ij < Njs0 ) THEN ! so side
629	IF( mpiRnei(nn_hls,jpso) > -1 ) lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
630	ENDIF
631	IF( ii >= Nis0 .AND. ii <= Nie0 .AND. ij > Nje0 ) THEN ! no side
632	IF( mpiRnei(nn_hls,jpno) > -1 ) lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
633	ENDIF
634	IF( ii < Nis0 .AND. ij < Njs0 ) THEN ! sw corner
635	IF( mpiRnei(nn_hls,jpsw) > -1 ) lrecv_bdyolr(ib_bdy,igrd,jpsw,ir) = .TRUE.
636	IF( mpiRnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 ) lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
637	ENDIF
638	IF( ii > Nie0 .AND. ij < Njs0 ) THEN ! se corner
639	IF( mpiRnei(nn_hls,jpse) > -1 ) lrecv_bdyolr(ib_bdy,igrd,jpse,ir) = .TRUE.
640	IF( mpiRnei(nn_hls,jpso) > -1 .AND. nn_comm == 1 ) lrecv_bdyolr(ib_bdy,igrd,jpso,ir) = .TRUE.
641	ENDIF
642	IF( ii < Nis0 .AND. ij > Nje0 ) THEN ! nw corner
643	IF( mpiRnei(nn_hls,jpnw) > -1 ) lrecv_bdyolr(ib_bdy,igrd,jpnw,ir) = .TRUE.
644	IF( mpiRnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 ) lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
645	ENDIF
646	IF( ii > Nie0 .AND. ij > Nje0 ) THEN ! ne corner
647	IF( mpiRnei(nn_hls,jpne) > -1 ) lrecv_bdyolr(ib_bdy,igrd,jpne,ir) = .TRUE.
648	IF( mpiRnei(nn_hls,jpno) > -1 .AND. nn_comm == 1 ) lrecv_bdyolr(ib_bdy,igrd,jpno,ir) = .TRUE.
649	ENDIF
650	!
651	END DO
652	END DO ! igrd
653
654	! Comment out for debug
655	!!$ DO ir = 0,1
656	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing, &
657	!!$ & lsend = lsend_bdyolr(ib_bdy,1,:,ir), lrecv = lrecv_bdyolr(ib_bdy,1,:,ir) )
658	!!$ IF(lwp) WRITE(numout,*) ' seb bdy debug olr T', ir ; CALL FLUSH(numout)
659	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing, &
660	!!$ & lsend = lsend_bdyolr(ib_bdy,2,:,ir), lrecv = lrecv_bdyolr(ib_bdy,2,:,ir) )
661	!!$ IF(lwp) WRITE(numout,*) ' seb bdy debug olr U', ir ; CALL FLUSH(numout)
662	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing, &
663	!!$ & lsend = lsend_bdyolr(ib_bdy,3,:,ir), lrecv = lrecv_bdyolr(ib_bdy,3,:,ir) )
664	!!$ IF(lwp) WRITE(numout,*) ' seb bdy debug olr V', ir ; CALL FLUSH(numout)
665	!!$ END DO
666
667	! Compute rim weights for FRS scheme
668	! ----------------------------------
669	DO igrd = 1, jpbgrd
670	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
671	ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) ) ! both rim 0 and rim 1 have the same weights
672	idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 ) ! tanh formulation
673	! idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic
674	! idx_bdy(ib_bdy)%nbw(ib,igrd) = REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)) ! linear
675	END DO
676	END DO
677
678	! Compute damping coefficients
679	! ----------------------------
680	DO igrd = 1, jpbgrd
681	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
682	ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) ) ! both rim 0 and rim 1 have the same damping coefficients
683	idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &
684	& (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))*2. ! quadratic
685	idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) &
686	& (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))*2. ! quadratic
687	END DO
688	END DO
689
690	END DO ! ib_bdy
691
692	! ------------------------------------------------------
693	! Initialise masks and find normal/tangential directions
694	! ------------------------------------------------------
695
696	! ------------------------------------------
697	! handle rim0, do as if rim 1 was free ocean
698	! ------------------------------------------
699
700	ztmask(:,:) = tmask(:,:,1) ; zumask(:,:) = umask(:,:,1) ; zvmask(:,:) = vmask(:,:,1)
701	! For the flagu/flagv calculation below we require a version of fmask without
702	! the land boundary condition (shlat) included:
703	DO_2D( 0, 0, 0, 0 )
704	zfmask(ji,jj) = ztmask(ji,jj ) * ztmask(ji+1,jj ) &
705	& * ztmask(ji,jj+1) * ztmask(ji+1,jj+1)
706	END_2D
707	CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp )
708
709	! Read global 2D mask at T-points: bdytmask
710	! -----------------------------------------
711	! bdytmask = 1 on the computational domain but not on open boundaries
712	! = 0 elsewhere
713
714	bdytmask(:,:) = ssmask(:,:)
715
716	! Derive mask on U and V grid from mask on T grid
717	DO_2D( 0, 0, 0, 0 )
718	bdyumask(ji,jj) = bdytmask(ji,jj) * bdytmask(ji+1,jj )
719	bdyvmask(ji,jj) = bdytmask(ji,jj) * bdytmask(ji ,jj+1)
720	END_2D
721	CALL lbc_lnk( 'bdyini', bdyumask, 'U', 1.0_wp , bdyvmask, 'V', 1.0_wp ) ! Lateral boundary cond.
722
723	! bdy masks are now set to zero on rim 0 points:
724	DO ib_bdy = 1, nb_bdy
725	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1) ! extent of rim 0
726	bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
727	END DO
728	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2) ! extent of rim 0
729	bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
730	END DO
731	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3) ! extent of rim 0
732	bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
733	END DO
734	END DO
735
736	CALL bdy_rim_treat( zumask, zvmask, zfmask, .true. ) ! compute flagu, flagv, ntreat on rim 0
737
738	! ------------------------------------
739	! handle rim1, do as if rim 0 was land
740	! ------------------------------------
741
742	! z[tuv]mask are now set to zero on rim 0 points:
743	DO ib_bdy = 1, nb_bdy
744	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1) ! extent of rim 0
745	ztmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
746	END DO
747	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2) ! extent of rim 0
748	zumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
749	END DO
750	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3) ! extent of rim 0
751	zvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
752	END DO
753	END DO
754
755	! Recompute zfmask
756	DO_2D( 0, 0, 0, 0 )
757	zfmask(ji,jj) = ztmask(ji,jj ) * ztmask(ji+1,jj ) &
758	& * ztmask(ji,jj+1) * ztmask(ji+1,jj+1)
759	END_2D
760	CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp )
761
762	! bdy masks are now set to zero on rim1 points:
763	DO ib_bdy = 1, nb_bdy
764	DO ib = idx_bdy(ib_bdy)%nblenrim0(1) + 1, idx_bdy(ib_bdy)%nblenrim(1) ! extent of rim 1
765	bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
766	END DO
767	DO ib = idx_bdy(ib_bdy)%nblenrim0(2) + 1, idx_bdy(ib_bdy)%nblenrim(2) ! extent of rim 1
768	bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
769	END DO
770	DO ib = idx_bdy(ib_bdy)%nblenrim0(3) + 1, idx_bdy(ib_bdy)%nblenrim(3) ! extent of rim 1
771	bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
772	END DO
773	END DO
774
775	CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. ) ! compute flagu, flagv, ntreat on rim 1
776	!
777	! Check which boundaries might need communication
778	ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,8,0:1), lrecv_bdyint(nb_bdy,jpbgrd,8,0:1) )
779	lsend_bdyint(:,:,:,:) = .false.
780	lrecv_bdyint(:,:,:,:) = .false.
781	ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,8,0:1), lrecv_bdyext(nb_bdy,jpbgrd,8,0:1) )
782	lsend_bdyext(:,:,:,:) = .false.
783	lrecv_bdyext(:,:,:,:) = .false.
784	!
785	DO ib_bdy = 1, nb_bdy
786	DO igrd = 1, jpbgrd
787	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
788	IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE
789	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
790	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
791	ir = idx_bdy(ib_bdy)%nbr(ib,igrd)
792	flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd))
793	flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd))
794	iibe = ii - flagu ! neighbouring point towards the exterior of the computational domain
795	ijbe = ij - flagv
796	iibi = ii + flagu ! neighbouring point towards the interior of the computational domain
797	ijbi = ij + flagv
798	CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 ) ! free ocean neighbours
799	!
800	! take care of the 4 sides
801	!
802	DO icnt = 1, 4
803	SELECT CASE( icnt )
804	! ... _____
805	CASE( 1 ) ! x: rim on rcvwe/sndea-side o\| :
806	! o: potential neighbour(s) o\|x :
807	! outside of the MPI domain ..o\|__:__
808	iRnei = jpwe ; iSnei = jpea
809	iiRst = 1 ; ijRst = 2 ! Rcv we-side starting point, excluding sw-corner
810	iiRnd = 1 ; ijRnd = jpj-1 ! Rcv we-side ending point, excluding nw-corner
811	iiSst = jpi-2*nn_hls+1 ; ijSst = 2 ! Snd ea-side starting point, excluding se-corner
812	iiSnd = jpi-2*nn_hls+1 ; ijSnd = jpj-1 ! Snd ea-side ending point, excluding ne-corner
813	iioutdir = -1 ; ijoutdir = -999 ! outside MPI domain: westward
814	! ______....
815	CASE( 2 ) ! x: rim on rcvea/sndwe-side : \|o
816	! o: potential neighbour(s) : x\|o
817	! outside of the MPI domain ___:__\|o..
818	iRnei = jpea ; iSnei = jpwe
819	iiRst = jpi ; ijRst = 2 ! Rcv ea-side starting point, excluding se-corner
820	iiRnd = jpi ; ijRnd = jpj-1 ! Rcv ea-side ending point, excluding ne-corner
821	iiSst = 2*nn_hls ; ijSst = 2 ! Snd we-side starting point, excluding sw-corner
822	iiSnd = 2*nn_hls ; ijSnd = jpj-1 ! Snd we-side ending point, excluding nw-corner
823	iioutdir = 1 ; ijoutdir = -999 ! outside MPI domain: eastward
824	!
825	CASE( 3 ) ! x: rim on rcvso/sndno-side \| \|
826	! o: potential neighbour(s) \|¨¨¨¨¨¨¨\|
827	! outside of the MPI domain \|___x___\|
828	! : o o o :
829	! : :
830	iRnei = jpso ; iSnei = jpno
831	iiRst = 2 ; ijRst = 1 ! Rcv so-side starting point, excluding sw-corner
832	iiRnd = jpi-1 ; ijRnd = 1 ! Rcv so-side ending point, excluding se-corner
833	iiSst = 2 ; ijSst = jpj-2*nn_hls+1 ! Snd no-side starting point, excluding nw-corner
834	iiSnd = jpi-1 ; ijSnd = jpj-2*nn_hls+1 ! Snd no-side ending point, excluding ne-corner
835	iioutdir = -999 ; ijoutdir = -1 ! outside MPI domain: southward
836	! : :
837	CASE( 4 ) ! x: rim on rcvno/sndso-side :_o_o_o_:
838	! o: potential neighbour(s) \| x \|
839	! outside of the MPI domain \| \|
840	! \|¨¨¨¨¨¨¨\|
841	iRnei = jpno ; iSnei = jpso
842	iiRst = 2 ; ijRst = jpj ! Rcv no-side starting point, excluding nw-corner
843	iiRnd = jpi-1 ; ijRnd = jpj ! Rcv no-side ending point, excluding ne-corner
844	iiSst = 2 ; ijSst = 2*nn_hls ! Snd so-side starting point, excluding sw-corner
845	iiSnd = jpi-1 ; ijSnd = 2*nn_hls ! Snd so-side ending point, excluding se-corner
846	iioutdir = -999 ; ijoutdir = 1 ! outside MPI domain: northward
847	END SELECT
848	!
849	IF( ii >= iiRst .AND. ii <= iiRnd .AND. ij >= ijRst .AND. ij <= ijRnd ) THEN ! rim point in recv side
850	iiout = ii+iioutdir ; ijout = ij+ijoutdir ! in which direction do we go outside of the MPI domain?
851	! take care of neighbourg(s) in the interior of the computational domain
852	IF( iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR. & ! Neib outside of the MPI domain
853	& ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN ! -> I cannot compute it -> recv it
854	IF( mpiRnei(nn_hls,iRnei) > -1 ) lrecv_bdyint(ib_bdy,igrd,iRnei,ir) = .TRUE.
855	ENDIF
856	! take care of neighbourg in the exterior of the computational domain
857	IF( iibe==iiout .OR. ijbe==ijout ) THEN ! Neib outside of the MPI domain -> I cannot compute it -> recv it
858	IF( mpiRnei(nn_hls,iRnei) > -1 ) lrecv_bdyext(ib_bdy,igrd,iRnei,ir) = .TRUE.
859	ENDIF
860	ENDIF
861
862	IF( ii >= iiSst .AND. ii <= iiSnd .AND. ij >= ijSst .AND. ij <= ijSnd ) THEN ! rim point in send side
863	iiout = ii+iioutdir ; ijout = ij+ijoutdir ! in which direction do we go outside of the nei MPI domain?
864	! take care of neighbourg(s) in the interior of the computational domain
865	IF( iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR. & ! Neib outside of nei MPI domain
866	& ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN ! -> nei cannot compute it
867	IF( mpiSnei(nn_hls,iSnei) > -1 ) lsend_bdyint(ib_bdy,igrd,iSnei,ir) = .TRUE. ! -> send to nei
868	ENDIF
869	! take care of neighbourg in the exterior of the computational domain
870	IF( iibe == iiout .OR. ijbe == ijout ) THEN ! Neib outside of the nei MPI domain -> nei cannot compute it
871	IF( mpiSnei(nn_hls,iSnei) > -1 ) lsend_bdyext(ib_bdy,igrd,iSnei,ir) = .TRUE. ! -> send to nei
872	ENDIF
873	END IF
874
875	END DO ! 4 sides
876	!
877	! specific treatment for the corners
878	!
879	DO icnt = 1, 4
880	SELECT CASE( icnt )
881	! ...\|....
882	CASE( 1 ) ! x: rim on sw-corner o\| :
883	! o: potential neighbour(s) o\|x__:__
884	! outside of the MPI domain o o o:
885	! :
886	iRdiag = jpsw ; iRsono = jpso ! Recv: for sw or so
887	iSdiag = jpne ; iSsono = jpno ! Send: to ne or no
888	iiRcorn = 1 ; ijRcorn = 1 ! receiving sw-corner
889	iiSdiag = jpi-2nn_hls+1 ; ijSdiag = jpj-2nn_hls+1 ! send to sw-corner of ne neighbourg
890	iiSsono = 1 ; ijSsono = jpj-2*nn_hls+1 ! send to sw-corner of no neighbourg
891	iioutdir = -1 ; ijoutdir = -1 ! outside MPI domain: westward or southward
892	! ....\|...
893	CASE( 2 ) ! x: rim on se-corner : \|o
894	! o: potential neighbour(s) __:__x\|o
895	! outside of the MPI domain :o o o
896	! :
897	iRdiag = jpse ; iRsono = jpso ! Recv: for se or so
898	iSdiag = jpnw ; iSsono = jpno ! Send: to nw or no
899	iiRcorn = jpi ; ijRcorn = 1 ! receiving se-corner
900	iiSdiag = 2nn_hls ; ijSdiag = jpj-2nn_hls+1 ! send to se-corner of nw neighbourg
901	iiSsono = jpi ; ijSsono = jpj-2*nn_hls+1 ! send to se-corner of no neighbourg
902	iioutdir = 1 ; ijoutdir = -1 ! outside MPI domain: eastward or southward
903	! :
904	! o o_o:___
905	CASE( 3 ) ! x: rim on nw-corner o\|x :
906	! o: potential neighbour(s) ..o\|...:
907	! outside of the MPI domain \|
908	iRdiag = jpnw ; iRsono = jpno ! Recv: for nw or no
909	iSdiag = jpse ; iSsono = jpso ! Send: to se or so
910	iiRcorn = 1 ; ijRcorn = jpj ! receiving nw-corner
911	iiSdiag = jpi-2nn_hls+1 ; ijSdiag = 2nn_hls ! send to nw-corner of se neighbourg
912	iiSsono = 1 ; ijSsono = 2*nn_hls ! send to nw-corner of so neighbourg
913	iioutdir = -1 ; ijoutdir = 1 ! outside MPI domain: westward or northward
914	! :
915	! ___:o_o o
916	CASE( 4 ) ! x: rim on ne-corner : x\|o
917	! o: potential neighbour(s) :...\|o...
918	! outside of the MPI domain \|
919	iRdiag = jpne ; iRsono = jpno ! Recv: for ne or no
920	iSdiag = jpsw ; iSsono = jpso ! Send: to sw or so
921	iiRcorn = jpi ; ijRcorn = jpj ! receiving ne-corner
922	iiSdiag = 2nn_hls ; ijSdiag = 2nn_hls ! send to ne-corner of sw neighbourg
923	iiSsono = jpi ; ijSsono = 2*nn_hls ! send to ne-corner of so neighbourg
924	iioutdir = 1 ; ijoutdir = 1 ! outside MPI domain: eastward or southward
925	END SELECT
926	!
927	! Check if we need to receive data for this rim point
928	IF( ii == iiRcorn .AND. ij == ijRcorn ) THEN ! the rim point is located on the corner for the MPI domain
929	iiout = ii+iioutdir ; ijout = ij+ijoutdir ! in which direction do we go outside of the MPI domain?
930	! take care of neighbourg(s) in the interior of the computational domain
931	IF( iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR. & ! Neib outside of the MPI domain
932	& ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN ! -> I cannot compute it -> recv it
933	IF( mpiRnei(nn_hls,iRdiag) > -1 ) lrecv_bdyint(ib_bdy,igrd,iRdiag,ir) = .TRUE. ! Receive directly from diagonal neighbourg
934	IF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 ) lrecv_bdyint(ib_bdy,igrd,iRsono,ir) = .TRUE. ! Receive through the South/North neighbourg
935	ENDIF
936	! take care of neighbourg in the exterior of the computational domain
937	IF( iibe==iiout .OR. ijbe==ijout ) THEN ! Neib outside of the MPI domain -> I cannot compute it -> recv it
938	IF( mpiRnei(nn_hls,iRdiag) > -1 ) lrecv_bdyext(ib_bdy,igrd,iRdiag,ir) = .TRUE. ! Receive directly from diagonal neighbourg
939	IF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 ) lrecv_bdyext(ib_bdy,igrd,iRsono,ir) = .TRUE. ! Receive through the South/North neighbourg
940	ENDIF
941	ENDIF
942	!
943	! Check if this rim point corresponds to the corner of one neighbourg. if yes, do we need to send data?
944	! Direct send to diag: Is this rim point the corner point of a diag neighbour with which we communicate?
945	IF( ii == iiSdiag .AND. ij == ijSdiag .AND. mpiSnei(nn_hls,iSdiag) > -1 ) THEN
946	iiout = ii+iioutdir ; ijout = ij+ijoutdir ! in which direction do we go outside of the nei MPI domain?
947	! take care of neighbourg(s) in the interior of the computational domain
948	IF( iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR. & ! Neib outside of diag nei MPI
949	& ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) & ! domain -> nei cannot compute it
950	& lsend_bdyint(ib_bdy,igrd,iSdiag,ir) = .TRUE. ! send rim point data to diag nei
951	! take care of neighbourg in the exterior of the computational domain
952	IF( iibe==iiout .OR. ijbe==ijout ) &
953	& lsend_bdyext(ib_bdy,igrd,iSdiag,ir) = .TRUE.
954	ENDIF
955	! Indirect send to diag (through so/no): rim point is the corner point of a so/no nei with which we communicate
956	IF( ii == iiSsono .AND. ij == ijSsono .AND. mpiSnei(nn_hls,iSsono) > -1 .AND. nn_comm == 1 ) THEN
957	iiout = ii+iioutdir ; ijout = ij+ijoutdir ! in which direction do we go outside of the nei MPI domain?
958	! take care of neighbourg(s) in the interior of the computational domain
959	IF( iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR. & ! Neib outside of so/no nei MPI
960	& ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) & ! domain -> nei cannot compute it
961	& lsend_bdyint(ib_bdy,igrd,iSsono,ir) = .TRUE. ! send rim point data to so/no nei
962	! take care of neighbourg in the exterior of the computational domain
963	IF( iibe==iiout .OR. ijbe==ijout ) &
964	& lsend_bdyext(ib_bdy,igrd,iSsono,ir) = .TRUE.
965	ENDIF
966	!
967	END DO ! 4 corners
968	END DO ! ib
969	END DO ! igrd
970
971	! Comment out for debug
972	!!$ DO ir = 0,1
973	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing, &
974	!!$ & lsend = lsend_bdyint(ib_bdy,1,:,ir), lrecv = lrecv_bdyint(ib_bdy,1,:,ir) )
975	!!$ IF(lwp) WRITE(numout,*) ' bdy debug int T', ir ; CALL FLUSH(numout)
976	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing, &
977	!!$ & lsend = lsend_bdyint(ib_bdy,2,:,ir), lrecv = lrecv_bdyint(ib_bdy,2,:,ir) )
978	!!$ IF(lwp) WRITE(numout,*) ' bdy debug int U', ir ; CALL FLUSH(numout)
979	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing, &
980	!!$ & lsend = lsend_bdyint(ib_bdy,3,:,ir), lrecv = lrecv_bdyint(ib_bdy,3,:,ir) )
981	!!$ IF(lwp) WRITE(numout,*) ' bdy debug int V', ir ; CALL FLUSH(numout)
982	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing, &
983	!!$ & lsend = lsend_bdyext(ib_bdy,1,:,ir), lrecv = lrecv_bdyext(ib_bdy,1,:,ir) )
984	!!$ IF(lwp) WRITE(numout,*) ' bdy debug ext T', ir ; CALL FLUSH(numout)
985	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing, &
986	!!$ & lsend = lsend_bdyext(ib_bdy,2,:,ir), lrecv = lrecv_bdyext(ib_bdy,2,:,ir) )
987	!!$ IF(lwp) WRITE(numout,*) ' bdy debug ext U', ir ; CALL FLUSH(numout)
988	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing, &
989	!!$ & lsend = lsend_bdyext(ib_bdy,3,:,ir), lrecv = lrecv_bdyext(ib_bdy,3,:,ir) )
990	!!$ IF(lwp) WRITE(numout,*) ' bdy debug ext V', ir ; CALL FLUSH(numout)
991	!!$ END DO
992
993	END DO ! ib_bdy
994
995	DO ib_bdy = 1,nb_bdy
996	IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &
997	& cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &
998	& cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo' ) THEN
999	DO igrd = 1, jpbgrd
1000	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1001	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1002	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1003	IF( mig0(ii) > 2 .AND. mig0(ii) < Ni0glo-2 .AND. mjg0(ij) > 2 .AND. mjg0(ij) < Nj0glo-2 ) THEN
1004	WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain'
1005	CALL ctl_stop( ctmp1 )
1006	END IF
1007	END DO
1008	END DO
1009	END IF
1010	END DO
1011	!
1012	DEALLOCATE( nbidta, nbjdta, nbrdta )
1013	!
1014	END SUBROUTINE bdy_def
1015
1016
1017	SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 )
1018	!!----------------------------------------------------------------------
1019	!! * ROUTINE bdy_rim_treat *
1020	!!
1021	!! ** Purpose : Initialize structures ( flagu, flagv, ntreat ) indicating how rim points
1022	!! are to be handled in the boundary condition treatment
1023	!!
1024	!! ** Method : - to handle rim 0 zmasks must indicate ocean points (set at one on rim 0 and rim 1 and interior)
1025	!! and bdymasks must be set at 0 on rim 0 (set at one on rim 1 and interior)
1026	!! (as if rim 1 was free ocean)
1027	!! - to handle rim 1 zmasks must be set at 0 on rim 0 (set at one on rim 1 and interior)
1028	!! and bdymasks must indicate free ocean points (set at one on interior)
1029	!! (as if rim 0 was land)
1030	!! - we can then check in which direction the interior of the computational domain is with the difference
1031	!! mask array values on both sides to compute flagu and flagv
1032	!! - and look at the ocean neighbours to compute ntreat
1033	!!----------------------------------------------------------------------
1034	REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: pumask, pvmask ! temporary u/v mask array
1035	REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: pfmask ! temporary fmask excluding coastal boundary condition (shlat)
1036	LOGICAL , INTENT (in ) :: lrim0 ! .true. -> rim 0 .false. -> rim 1
1037	INTEGER :: ib_bdy, ii, ij, igrd, ib, icount ! dummy loop indices
1038	INTEGER :: i_offset, j_offset, inn ! local integer
1039	INTEGER :: ibeg, iend ! local integer
1040	LOGICAL :: llnon, llson, llean, llwen ! local logicals indicating the presence of a ocean neighbour
1041	REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! pointer to 2D mask fields
1042	REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars
1043	CHARACTER(LEN=1), DIMENSION(jpbgrd) :: cgrid
1044	REAL(wp) , DIMENSION(jpi,jpj) :: ztmp
1045	!!----------------------------------------------------------------------
1046
1047	cgrid = (/'t','u','v'/)
1048
1049	DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components
1050
1051	DO igrd = 1, jpbgrd
1052
1053	IF( lrim0 ) THEN ! extent of rim 0
1054	ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
1055	ELSE ! extent of rim 1
1056	ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1 ; iend = idx_bdy(ib_bdy)%nblenrim(igrd)
1057	END IF
1058
1059	! Calculate relationship of U direction to the local orientation of the boundary
1060	! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
1061	! flagu = 0 : u is tangential
1062	! flagu = 1 : u is normal to the boundary and is direction is inward
1063	SELECT CASE( igrd )
1064	CASE( 1 ) ; zmask => pumask ; i_offset = 0 ! U(i-1) T(i) U(i )
1065	CASE( 2 ) ; zmask => bdytmask ; i_offset = 1 ! T(i ) U(i) T(i+1)
1066	CASE( 3 ) ; zmask => pfmask ; i_offset = 0 ! F(i-1) V(i) F(i )
1067	END SELECT
1068	icount = 0
1069	ztmp(:,:) = -999._wp
1070	DO ib = ibeg, iend
1071	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1072	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1073	IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 ) CYCLE ! call lbc_lnk -> no need to compute these pts
1074	zwfl = zmask(ii+i_offset-1,ij)
1075	zefl = zmask(ii+i_offset ,ij)
1076	! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
1077	IF( i_offset == 1 .and. zefl + zwfl == 2._wp ) THEN
1078	icount = icount + 1
1079	IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
1080	ELSE
1081	ztmp(ii,ij) = -zwfl + zefl
1082	ENDIF
1083	END DO
1084	IF( icount /= 0 ) THEN
1085	WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,', &
1086	' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1087	CALL ctl_stop( ctmp1 )
1088	ENDIF
1089	SELECT CASE( igrd )
1090	CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1091	CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1092	CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1093	END SELECT
1094	DO ib = ibeg, iend
1095	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1096	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1097	idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij)
1098	END DO
1099
1100	! Calculate relationship of V direction to the local orientation of the boundary
1101	! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
1102	! flagv = 0 : v is tangential
1103	! flagv = 1 : v is normal to the boundary and is direction is inward
1104	SELECT CASE( igrd )
1105	CASE( 1 ) ; zmask => pvmask ; j_offset = 0
1106	CASE( 2 ) ; zmask => pfmask ; j_offset = 0
1107	CASE( 3 ) ; zmask => bdytmask ; j_offset = 1
1108	END SELECT
1109	icount = 0
1110	ztmp(:,:) = -999._wp
1111	DO ib = ibeg, iend
1112	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1113	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1114	IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 ) CYCLE ! call lbc_lnk -> no need to compute these pts
1115	zsfl = zmask(ii,ij+j_offset-1)
1116	znfl = zmask(ii,ij+j_offset )
1117	! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
1118	IF( j_offset == 1 .and. znfl + zsfl == 2._wp ) THEN
1119	IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
1120	icount = icount + 1
1121	ELSE
1122	ztmp(ii,ij) = -zsfl + znfl
1123	END IF
1124	END DO
1125	IF( icount /= 0 ) THEN
1126	WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,', &
1127	' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1128	CALL ctl_stop( ctmp1 )
1129	ENDIF
1130	SELECT CASE( igrd )
1131	CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1132	CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1133	CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1134	END SELECT
1135	DO ib = ibeg, iend
1136	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1137	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1138	idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij)
1139	END DO
1140
1141	! Calculate ntreat
1142	SELECT CASE( igrd )
1143	CASE( 1 ) ; zmask => bdytmask
1144	CASE( 2 ) ; zmask => bdyumask
1145	CASE( 3 ) ; zmask => bdyvmask
1146	END SELECT
1147	ztmp(:,:) = -999._wp
1148	DO ib = ibeg, iend
1149	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1150	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1151	IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 ) CYCLE ! call lbc_lnk -> no need to compute these pts
1152	llnon = zmask(ii ,ij+1) == 1._wp
1153	llson = zmask(ii ,ij-1) == 1._wp
1154	llean = zmask(ii+1,ij ) == 1._wp
1155	llwen = zmask(ii-1,ij ) == 1._wp
1156	inn = COUNT( (/ llnon, llson, llean, llwen /) )
1157	IF( inn == 0 ) THEN ! no neighbours -> interior of a corner or cluster of rim points
1158	! ! ! _____ ! _____ ! __ __
1159	! 1 \| o ! 2 o \| ! 3 \| x ! 4 x \| ! \| \| -> error
1160	! \|_x_ _ ! _ _x_\| ! \| o ! o \| ! \|x_x\|
1161	IF( zmask(ii+1,ij+1) == 1._wp ) THEN ; ztmp(ii,ij) = 1._wp
1162	ELSEIF( zmask(ii-1,ij+1) == 1._wp ) THEN ; ztmp(ii,ij) = 2._wp
1163	ELSEIF( zmask(ii+1,ij-1) == 1._wp ) THEN ; ztmp(ii,ij) = 3._wp
1164	ELSEIF( zmask(ii-1,ij-1) == 1._wp ) THEN ; ztmp(ii,ij) = 4._wp
1165	ELSE ; ztmp(ii,ij) = -1._wp
1166	WRITE(ctmp1,*) 'Problem with ',cgrid(igrd) ,' grid point', ii, ij, &
1167	' on boundary set ', ib_bdy, ' has no free ocean neighbour'
1168	IF( lrim0 ) THEN
1169	WRITE(ctmp2,*) ' There seems to be a cluster of rim 0 points.'
1170	ELSE
1171	WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.'
1172	END IF
1173	CALL ctl_warn( ctmp1, ctmp2 )
1174	END IF
1175	END IF
1176	IF( inn == 1 ) THEN ! middle of linear bdy or incomplete corner ! ___ o
1177	! \| ! \| ! o ! ______ ! \|x___
1178	! 5 \| x o ! 6 o x \| ! 7 __x__ ! 8 x
1179	! \| ! \| ! ! o
1180	IF( llean ) ztmp(ii,ij) = 5._wp
1181	IF( llwen ) ztmp(ii,ij) = 6._wp
1182	IF( llnon ) ztmp(ii,ij) = 7._wp
1183	IF( llson ) ztmp(ii,ij) = 8._wp
1184	END IF
1185	IF( inn == 2 ) THEN ! exterior of a corner
1186	! o ! o ! _____\| ! \|_____
1187	! 9 ____x o ! 10 o x___ ! 11 x o ! 12 o x
1188	! \| ! \| ! o ! o
1189	IF( llnon .AND. llean ) ztmp(ii,ij) = 9._wp
1190	IF( llnon .AND. llwen ) ztmp(ii,ij) = 10._wp
1191	IF( llson .AND. llean ) ztmp(ii,ij) = 11._wp
1192	IF( llson .AND. llwen ) ztmp(ii,ij) = 12._wp
1193	END IF
1194	IF( inn == 3 ) THEN ! 3 neighbours __ __
1195	! \|_ o ! o _\| ! \|_\| ! o
1196	! 13 _\| x o ! 14 o x \|_ ! 15 o x o ! 16 o x o
1197	! \| o ! o \| ! o ! __\|¨\|__
1198	IF( llnon .AND. llean .AND. llson ) ztmp(ii,ij) = 13._wp
1199	IF( llnon .AND. llwen .AND. llson ) ztmp(ii,ij) = 14._wp
1200	IF( llwen .AND. llson .AND. llean ) ztmp(ii,ij) = 15._wp
1201	IF( llwen .AND. llnon .AND. llean ) ztmp(ii,ij) = 16._wp
1202	END IF
1203	IF( inn == 4 ) THEN
1204	WRITE(ctmp1,*) 'Problem with ',cgrid(igrd) ,' grid point', ii, ij, &
1205	' on boundary set ', ib_bdy, ' have 4 neighbours'
1206	CALL ctl_stop( ctmp1 )
1207	END IF
1208	END DO
1209	SELECT CASE( igrd )
1210	CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1211	CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1212	CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1213	END SELECT
1214	DO ib = ibeg, iend
1215	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1216	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1217	idx_bdy(ib_bdy)%ntreat(ib,igrd) = NINT(ztmp(ii,ij))
1218	END DO
1219	!
1220	END DO ! jpbgrd
1221	!
1222	END DO ! ib_bdy
1223
1224	END SUBROUTINE bdy_rim_treat
1225
1226
1227	SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )
1228	!!----------------------------------------------------------------------
1229	!! * ROUTINE find_neib *
1230	!!
1231	!! ** Purpose : get ii1, ij1, ii2, ij2, ii3, ij3, the indices of
1232	!! the free ocean neighbours of (ii,ij) for bdy treatment
1233	!!
1234	!! ** Method : use itreat input to select a case
1235	!! N.B. ntreat is defined for all bdy points in routine bdy_rim_treat
1236	!!
1237	!!----------------------------------------------------------------------
1238	INTEGER, INTENT(in ) :: ii, ij, itreat
1239	INTEGER, INTENT( out) :: ii1, ij1, ii2, ij2, ii3, ij3
1240	!!----------------------------------------------------------------------
1241	SELECT CASE( itreat ) ! points that will be used by bdy routines, -1 will be discarded
1242	! ! ! _____ ! _____
1243	! 1 \| o ! 2 o \| ! 3 \| x ! 4 x \|
1244	! \|_x_ _ ! _ _x_\| ! \| o ! o \|
1245	CASE( 1 ) ; ii1 = ii+1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1246	CASE( 2 ) ; ii1 = ii-1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1247	CASE( 3 ) ; ii1 = ii+1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1248	CASE( 4 ) ; ii1 = ii-1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1249	! \| ! \| ! o ! ______ ! or incomplete corner
1250	! 5 \| x o ! 6 o x \| ! 7 __x__ ! 8 x ! 7 ____ o
1251	! \| ! \| ! ! o ! \|x___
1252	CASE( 5 ) ; ii1 = ii+1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1253	CASE( 6 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1254	CASE( 7 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1255	CASE( 8 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1256	! o ! o ! _____\| ! \|_____
1257	! 9 ____x o ! 10 o x___ ! 11 x o ! 12 o x
1258	! \| ! \| ! o ! o
1259	CASE( 9 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1260	CASE( 10 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1261	CASE( 11 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1262	CASE( 12 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1263	! \|_ o ! o _\| ! ¨¨\|_\|¨¨ ! o
1264	! 13 _\| x o ! 14 o x \|_ ! 15 o x o ! 16 o x o
1265	! \| o ! o \| ! o ! __\|¨\|__
1266	CASE( 13 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1
1267	CASE( 14 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1
1268	CASE( 15 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij-1 ; ii3 = ii+1 ; ij3 = ij
1269	CASE( 16 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij+1 ; ii3 = ii+1 ; ij3 = ij
1270	END SELECT
1271	END SUBROUTINE find_neib
1272
1273
1274	SUBROUTINE bdy_read_seg( kb_bdy, knblendta )
1275	!!----------------------------------------------------------------------
1276	!! * ROUTINE bdy_read_seg *
1277	!!
1278	!! ** Purpose : build bdy coordinates with segments defined in namelist
1279	!!
1280	!! ** Method : read namelist nambdy_index blocks
1281	!!
1282	!!----------------------------------------------------------------------
1283	INTEGER , INTENT (in ) :: kb_bdy ! bdy number
1284	INTEGER, DIMENSION(jpbgrd), INTENT ( out) :: knblendta ! length of index arrays
1285	!!
1286	INTEGER :: ios ! Local integer output status for namelist read
1287	INTEGER :: nbdyind, nbdybeg, nbdyend
1288	INTEGER :: nbdy_count, nbdy_rdstart, nbdy_loc
1289	CHARACTER(LEN=1) :: ctypebdy ! - -
1290	CHARACTER(LEN=50):: cerrmsg ! - -
1291	NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
1292	!!----------------------------------------------------------------------
1293	! Need to support possibility of reading more than one nambdy_index from
1294	! the namelist_cfg internal file.
1295	! Do this by finding the kb_bdy'th occurence of nambdy_index in the
1296	! character buffer as the starting point.
1297	nbdy_rdstart = 1
1298	DO nbdy_count = 1, kb_bdy
1299	nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_index' )
1300	IF( nbdy_loc .GT. 0 ) THEN
1301	nbdy_rdstart = nbdy_rdstart + nbdy_loc
1302	ELSE
1303	WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',kb_bdy,' of nambdy_index not found'
1304	ios = -1
1305	CALL ctl_nam ( ios , cerrmsg )
1306	ENDIF
1307	END DO
1308	nbdy_rdstart = MAX( 1, nbdy_rdstart - 2 )
1309	READ ( numnam_cfg( nbdy_rdstart: ), nambdy_index, IOSTAT = ios, ERR = 904)
1310	904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' )
1311	IF(lwm) WRITE ( numond, nambdy_index )
1312
1313	SELECT CASE ( TRIM(ctypebdy) )
1314	CASE( 'N' )
1315	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1316	nbdyind = Nj0glo - 2 ! set boundary to whole side of model domain.
1317	nbdybeg = 2
1318	nbdyend = Ni0glo - 1
1319	ENDIF
1320	nbdysegn = nbdysegn + 1
1321	npckgn(nbdysegn) = kb_bdy ! Save bdy package number
1322	jpjnob(nbdysegn) = nbdyind
1323	jpindt(nbdysegn) = nbdybeg
1324	jpinft(nbdysegn) = nbdyend
1325	!
1326	CASE( 'S' )
1327	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1328	nbdyind = 2 ! set boundary to whole side of model domain.
1329	nbdybeg = 2
1330	nbdyend = Ni0glo - 1
1331	ENDIF
1332	nbdysegs = nbdysegs + 1
1333	npckgs(nbdysegs) = kb_bdy ! Save bdy package number
1334	jpjsob(nbdysegs) = nbdyind
1335	jpisdt(nbdysegs) = nbdybeg
1336	jpisft(nbdysegs) = nbdyend
1337	!
1338	CASE( 'E' )
1339	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1340	nbdyind = Ni0glo - 2 ! set boundary to whole side of model domain.
1341	nbdybeg = 2
1342	nbdyend = Nj0glo - 1
1343	ENDIF
1344	nbdysege = nbdysege + 1
1345	npckge(nbdysege) = kb_bdy ! Save bdy package number
1346	jpieob(nbdysege) = nbdyind
1347	jpjedt(nbdysege) = nbdybeg
1348	jpjeft(nbdysege) = nbdyend
1349	!
1350	CASE( 'W' )
1351	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1352	nbdyind = 2 ! set boundary to whole side of model domain.
1353	nbdybeg = 2
1354	nbdyend = Nj0glo - 1
1355	ENDIF
1356	nbdysegw = nbdysegw + 1
1357	npckgw(nbdysegw) = kb_bdy ! Save bdy package number
1358	jpiwob(nbdysegw) = nbdyind
1359	jpjwdt(nbdysegw) = nbdybeg
1360	jpjwft(nbdysegw) = nbdyend
1361	!
1362	CASE DEFAULT ; CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
1363	END SELECT
1364
1365	! For simplicity we assume that in case of straight bdy, arrays have the same length
1366	! (even if it is true that last tangential velocity points
1367	! are useless). This simplifies a little bit boundary data format (and agrees with format
1368	! used so far in obc package)
1369
1370	knblendta(1:jpbgrd) = (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy)
1371
1372	END SUBROUTINE bdy_read_seg
1373
1374
1375	SUBROUTINE bdy_ctl_seg
1376	!!----------------------------------------------------------------------
1377	!! * ROUTINE bdy_ctl_seg *
1378	!!
1379	!! ** Purpose : Check straight open boundary segments location
1380	!!
1381	!! ** Method : - Look for open boundary corners
1382	!! - Check that segments start or end on land
1383	!!----------------------------------------------------------------------
1384	INTEGER :: ib, ib1, ib2, ji ,jj, itest
1385	INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns
1386	REAL(wp), DIMENSION(2) :: ztestmask
1387	!!----------------------------------------------------------------------
1388	!
1389	IF (lwp) WRITE(numout,*) ' '
1390	IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1391	IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1392	!
1393	IF(lwp) WRITE(numout,*) 'Number of east segments : ', nbdysege
1394	IF(lwp) WRITE(numout,*) 'Number of west segments : ', nbdysegw
1395	IF(lwp) WRITE(numout,*) 'Number of north segments : ', nbdysegn
1396	IF(lwp) WRITE(numout,*) 'Number of south segments : ', nbdysegs
1397	!
1398	! 1. Check bounds
1399	!----------------
1400	DO ib = 1, nbdysegn
1401	IF (lwp) WRITE(numout,) '*check north seg bounds pckg: ', npckgn(ib)
1402	IF ((jpjnob(ib).ge.Nj0glo-1).or.&
1403	&(jpjnob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1404	IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1405	IF (jpindt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1406	IF (jpinft(ib).gt.Ni0glo) CALL ctl_stop( 'End index out of domain' )
1407	END DO
1408	!
1409	DO ib = 1, nbdysegs
1410	IF (lwp) WRITE(numout,) '*check south seg bounds pckg: ', npckgs(ib)
1411	IF ((jpjsob(ib).ge.Nj0glo-1).or.&
1412	&(jpjsob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1413	IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1414	IF (jpisdt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1415	IF (jpisft(ib).gt.Ni0glo) CALL ctl_stop( 'End index out of domain' )
1416	END DO
1417	!
1418	DO ib = 1, nbdysege
1419	IF (lwp) WRITE(numout,) '*check east seg bounds pckg: ', npckge(ib)
1420	IF ((jpieob(ib).ge.Ni0glo-1).or.&
1421	&(jpieob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1422	IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1423	IF (jpjedt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1424	IF (jpjeft(ib).gt.Nj0glo) CALL ctl_stop( 'End index out of domain' )
1425	END DO
1426	!
1427	DO ib = 1, nbdysegw
1428	IF (lwp) WRITE(numout,) '*check west seg bounds pckg: ', npckgw(ib)
1429	IF ((jpiwob(ib).ge.Ni0glo-1).or.&
1430	&(jpiwob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1431	IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1432	IF (jpjwdt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1433	IF (jpjwft(ib).gt.Nj0glo) CALL ctl_stop( 'End index out of domain' )
1434	ENDDO
1435	!
1436	! 2. Look for segment crossings
1437	!------------------------------
1438	IF (lwp) WRITE(numout,) '*Look for segments corners :'
1439	!
1440	itest = 0 ! corner number
1441	!
1442	! flag to detect if start or end of open boundary belongs to a corner
1443	! if not (=0), it must be on land.
1444	! if a corner is detected, save bdy package number for further tests
1445	icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1446	! South/West crossings
1447	IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1448	DO ib1 = 1, nbdysegw
1449	DO ib2 = 1, nbdysegs
1450	IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1451	& ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1452	& ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1453	& ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1454	IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN
1455	! We have a possible South-West corner
1456	! WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1457	! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1458	icornw(ib1,1) = npckgs(ib2)
1459	icorns(ib2,1) = npckgw(ib1)
1460	ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1461	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1462	& jpisft(ib2), jpjwft(ib1)
1463	WRITE(ctmp2,*) ' Not allowed yet'
1464	WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1465	& ' and South segment: ',npckgs(ib2)
1466	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1467	ELSE
1468	WRITE(ctmp1,*) ' Check South and West Open boundary indices'
1469	WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , &
1470	& ' and South segment: ',npckgs(ib2)
1471	CALL ctl_stop( ctmp1, ctmp2 )
1472	END IF
1473	END IF
1474	END DO
1475	END DO
1476	END IF
1477	!
1478	! South/East crossings
1479	IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1480	DO ib1 = 1, nbdysege
1481	DO ib2 = 1, nbdysegs
1482	IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1483	& ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1484	& ( jpjedt(ib1)<=jpjsob(ib2) ).AND. &
1485	& ( jpjeft(ib1)>=jpjsob(ib2) )) THEN
1486	IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1487	! We have a possible South-East corner
1488	! WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1489	! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1490	icorne(ib1,1) = npckgs(ib2)
1491	icorns(ib2,2) = npckge(ib1)
1492	ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1493	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1494	& jpisdt(ib2), jpjeft(ib1)
1495	WRITE(ctmp2,*) ' Not allowed yet'
1496	WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1497	& ' and South segment: ',npckgs(ib2)
1498	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1499	ELSE
1500	WRITE(ctmp1,*) ' Check South and East Open boundary indices'
1501	WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1502	& ' and South segment: ',npckgs(ib2)
1503	CALL ctl_stop( ctmp1, ctmp2 )
1504	END IF
1505	END IF
1506	END DO
1507	END DO
1508	END IF
1509	!
1510	! North/West crossings
1511	IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1512	DO ib1 = 1, nbdysegw
1513	DO ib2 = 1, nbdysegn
1514	IF (( jpindt(ib2)<=jpiwob(ib1) ).AND. &
1515	& ( jpinft(ib2)>=jpiwob(ib1) ).AND. &
1516	& ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1517	& ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1518	IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1519	! We have a possible North-West corner
1520	! WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1521	! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1522	icornw(ib1,2) = npckgn(ib2)
1523	icornn(ib2,1) = npckgw(ib1)
1524	ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1525	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1526	& jpinft(ib2), jpjwdt(ib1)
1527	WRITE(ctmp2,*) ' Not allowed yet'
1528	WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1529	& ' and North segment: ',npckgn(ib2)
1530	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1531	ELSE
1532	WRITE(ctmp1,*) ' Check North and West Open boundary indices'
1533	WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1534	& ' and North segment: ',npckgn(ib2)
1535	CALL ctl_stop( ctmp1, ctmp2 )
1536	END IF
1537	END IF
1538	END DO
1539	END DO
1540	END IF
1541	!
1542	! North/East crossings
1543	IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1544	DO ib1 = 1, nbdysege
1545	DO ib2 = 1, nbdysegn
1546	IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1547	& ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1548	& ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1549	& ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1550	IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1551	! We have a possible North-East corner
1552	! WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1553	! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1554	icorne(ib1,2) = npckgn(ib2)
1555	icornn(ib2,2) = npckge(ib1)
1556	ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1557	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1558	& jpindt(ib2), jpjedt(ib1)
1559	WRITE(ctmp2,*) ' Not allowed yet'
1560	WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1561	& ' and North segment: ',npckgn(ib2)
1562	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1563	ELSE
1564	WRITE(ctmp1,*) ' Check North and East Open boundary indices'
1565	WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1566	& ' and North segment: ',npckgn(ib2)
1567	CALL ctl_stop( ctmp1, ctmp2 )
1568	END IF
1569	END IF
1570	END DO
1571	END DO
1572	END IF
1573	!
1574	! 3. Check if segment extremities are on land
1575	!--------------------------------------------
1576	!
1577	! West segments
1578	DO ib = 1, nbdysegw
1579	! get mask at boundary extremities:
1580	ztestmask(1:2)=0.
1581	DO ji = 1, jpi
1582	DO jj = 1, jpj
1583	IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwdt(ib) ) ztestmask(1) = tmask(ji,jj,1)
1584	IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwft(ib) ) ztestmask(2) = tmask(ji,jj,1)
1585	END DO
1586	END DO
1587	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1588
1589	IF (ztestmask(1)==1) THEN
1590	IF (icornw(ib,1)==0) THEN
1591	WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1592	CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1593	ELSE
1594	! This is a corner
1595	IF(lwp) WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1596	CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1597	itest=itest+1
1598	ENDIF
1599	ENDIF
1600	IF (ztestmask(2)==1) THEN
1601	IF (icornw(ib,2)==0) THEN
1602	WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1603	CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' )
1604	ELSE
1605	! This is a corner
1606	IF(lwp) WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1607	CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1608	itest=itest+1
1609	ENDIF
1610	ENDIF
1611	END DO
1612	!
1613	! East segments
1614	DO ib = 1, nbdysege
1615	! get mask at boundary extremities:
1616	ztestmask(1:2)=0.
1617	DO ji = 1, jpi
1618	DO jj = 1, jpj
1619	IF( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjedt(ib) ) ztestmask(1) = tmask(ji,jj,1)
1620	IF( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjeft(ib) ) ztestmask(2) = tmask(ji,jj,1)
1621	END DO
1622	END DO
1623	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1624
1625	IF (ztestmask(1)==1) THEN
1626	IF (icorne(ib,1)==0) THEN
1627	WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1628	CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1629	ELSE
1630	! This is a corner
1631	IF(lwp) WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1632	CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1633	itest=itest+1
1634	ENDIF
1635	ENDIF
1636	IF (ztestmask(2)==1) THEN
1637	IF (icorne(ib,2)==0) THEN
1638	WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1639	CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1640	ELSE
1641	! This is a corner
1642	IF(lwp) WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1643	CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1644	itest=itest+1
1645	ENDIF
1646	ENDIF
1647	END DO
1648	!
1649	! South segments
1650	DO ib = 1, nbdysegs
1651	! get mask at boundary extremities:
1652	ztestmask(1:2)=0.
1653	DO ji = 1, jpi
1654	DO jj = 1, jpj
1655	IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisdt(ib) ) ztestmask(1) = tmask(ji,jj,1)
1656	IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisft(ib) ) ztestmask(2) = tmask(ji,jj,1)
1657	END DO
1658	END DO
1659	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1660
1661	IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1662	WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1663	CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1664	ENDIF
1665	IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1666	WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1667	CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1668	ENDIF
1669	END DO
1670	!
1671	! North segments
1672	DO ib = 1, nbdysegn
1673	! get mask at boundary extremities:
1674	ztestmask(1:2)=0.
1675	DO ji = 1, jpi
1676	DO jj = 1, jpj
1677	IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpindt(ib) ) ztestmask(1) = tmask(ji,jj,1)
1678	IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpinft(ib) ) ztestmask(2) = tmask(ji,jj,1)
1679	END DO
1680	END DO
1681	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1682
1683	IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1684	WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1685	CALL ctl_stop( ctmp1, ' does not start on land' )
1686	ENDIF
1687	IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1688	WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1689	CALL ctl_stop( ctmp1, ' does not end on land' )
1690	ENDIF
1691	END DO
1692	!
1693	IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1694	!
1695	! Other tests TBD:
1696	! segments completly on land
1697	! optimized open boundary array length according to landmask
1698	! Nudging layers that overlap with interior domain
1699	!
1700	END SUBROUTINE bdy_ctl_seg
1701
1702
1703	SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta )
1704	!!----------------------------------------------------------------------
1705	!! * ROUTINE bdy_coords_seg *
1706	!!
1707	!! ** Purpose : build nbidta, nbidta, nbrdta for bdy built with segments
1708	!!
1709	!! ** Method :
1710	!!
1711	!!----------------------------------------------------------------------
1712	INTEGER, DIMENSION(:,:,:), intent( out) :: nbidta, nbjdta, nbrdta ! Index arrays: i and j indices of bdy dta
1713	!!
1714	INTEGER :: ii, ij, ir, iseg
1715	INTEGER :: igrd ! grid type (t=1, u=2, v=3)
1716	INTEGER :: icount !
1717	INTEGER :: ib_bdy ! bdy number
1718	!!----------------------------------------------------------------------
1719
1720	! East
1721	!-----
1722	DO iseg = 1, nbdysege
1723	ib_bdy = npckge(iseg)
1724	!
1725	! ------------ T points -------------
1726	igrd=1
1727	icount=0
1728	DO ir = 1, nn_rimwidth(ib_bdy)
1729	DO ij = jpjedt(iseg), jpjeft(iseg)
1730	icount = icount + 1
1731	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir + nn_hls
1732	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1733	nbrdta(icount, igrd, ib_bdy) = ir
1734	ENDDO
1735	ENDDO
1736	!
1737	! ------------ U points -------------
1738	igrd=2
1739	icount=0
1740	DO ir = 1, nn_rimwidth(ib_bdy)
1741	DO ij = jpjedt(iseg), jpjeft(iseg)
1742	icount = icount + 1
1743	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir + nn_hls
1744	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1745	nbrdta(icount, igrd, ib_bdy) = ir
1746	ENDDO
1747	ENDDO
1748	!
1749	! ------------ V points -------------
1750	igrd=3
1751	icount=0
1752	DO ir = 1, nn_rimwidth(ib_bdy)
1753	! DO ij = jpjedt(iseg), jpjeft(iseg) - 1
1754	DO ij = jpjedt(iseg), jpjeft(iseg)
1755	icount = icount + 1
1756	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir + nn_hls
1757	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1758	nbrdta(icount, igrd, ib_bdy) = ir
1759	ENDDO
1760	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1761	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1762	ENDDO
1763	ENDDO
1764	!
1765	! West
1766	!-----
1767	DO iseg = 1, nbdysegw
1768	ib_bdy = npckgw(iseg)
1769	!
1770	! ------------ T points -------------
1771	igrd=1
1772	icount=0
1773	DO ir = 1, nn_rimwidth(ib_bdy)
1774	DO ij = jpjwdt(iseg), jpjwft(iseg)
1775	icount = icount + 1
1776	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1777	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1778	nbrdta(icount, igrd, ib_bdy) = ir
1779	ENDDO
1780	ENDDO
1781	!
1782	! ------------ U points -------------
1783	igrd=2
1784	icount=0
1785	DO ir = 1, nn_rimwidth(ib_bdy)
1786	DO ij = jpjwdt(iseg), jpjwft(iseg)
1787	icount = icount + 1
1788	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1789	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1790	nbrdta(icount, igrd, ib_bdy) = ir
1791	ENDDO
1792	ENDDO
1793	!
1794	! ------------ V points -------------
1795	igrd=3
1796	icount=0
1797	DO ir = 1, nn_rimwidth(ib_bdy)
1798	! DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
1799	DO ij = jpjwdt(iseg), jpjwft(iseg)
1800	icount = icount + 1
1801	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1802	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1803	nbrdta(icount, igrd, ib_bdy) = ir
1804	ENDDO
1805	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1806	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1807	ENDDO
1808	ENDDO
1809	!
1810	! North
1811	!-----
1812	DO iseg = 1, nbdysegn
1813	ib_bdy = npckgn(iseg)
1814	!
1815	! ------------ T points -------------
1816	igrd=1
1817	icount=0
1818	DO ir = 1, nn_rimwidth(ib_bdy)
1819	DO ii = jpindt(iseg), jpinft(iseg)
1820	icount = icount + 1
1821	nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1822	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir + nn_hls
1823	nbrdta(icount, igrd, ib_bdy) = ir
1824	ENDDO
1825	ENDDO
1826	!
1827	! ------------ U points -------------
1828	igrd=2
1829	icount=0
1830	DO ir = 1, nn_rimwidth(ib_bdy)
1831	! DO ii = jpindt(iseg), jpinft(iseg) - 1
1832	DO ii = jpindt(iseg), jpinft(iseg)
1833	icount = icount + 1
1834	nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1835	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir + nn_hls
1836	nbrdta(icount, igrd, ib_bdy) = ir
1837	ENDDO
1838	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1839	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1840	ENDDO
1841	!
1842	! ------------ V points -------------
1843	igrd=3
1844	icount=0
1845	DO ir = 1, nn_rimwidth(ib_bdy)
1846	DO ii = jpindt(iseg), jpinft(iseg)
1847	icount = icount + 1
1848	nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1849	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir + nn_hls
1850	nbrdta(icount, igrd, ib_bdy) = ir
1851	ENDDO
1852	ENDDO
1853	ENDDO
1854	!
1855	! South
1856	!-----
1857	DO iseg = 1, nbdysegs
1858	ib_bdy = npckgs(iseg)
1859	!
1860	! ------------ T points -------------
1861	igrd=1
1862	icount=0
1863	DO ir = 1, nn_rimwidth(ib_bdy)
1864	DO ii = jpisdt(iseg), jpisft(iseg)
1865	icount = icount + 1
1866	nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1867	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1868	nbrdta(icount, igrd, ib_bdy) = ir
1869	ENDDO
1870	ENDDO
1871	!
1872	! ------------ U points -------------
1873	igrd=2
1874	icount=0
1875	DO ir = 1, nn_rimwidth(ib_bdy)
1876	! DO ii = jpisdt(iseg), jpisft(iseg) - 1
1877	DO ii = jpisdt(iseg), jpisft(iseg)
1878	icount = icount + 1
1879	nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1880	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1881	nbrdta(icount, igrd, ib_bdy) = ir
1882	ENDDO
1883	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1884	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1885	ENDDO
1886	!
1887	! ------------ V points -------------
1888	igrd=3
1889	icount=0
1890	DO ir = 1, nn_rimwidth(ib_bdy)
1891	DO ii = jpisdt(iseg), jpisft(iseg)
1892	icount = icount + 1
1893	nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1894	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1895	nbrdta(icount, igrd, ib_bdy) = ir
1896	ENDDO
1897	ENDDO
1898	ENDDO
1899
1900
1901	END SUBROUTINE bdy_coords_seg
1902
1903
1904	SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1905	!!----------------------------------------------------------------------
1906	!! * ROUTINE bdy_ctl_corn *
1907	!!
1908	!! ** Purpose : Check numerical schemes consistency between
1909	!! segments having a common corner
1910	!!
1911	!! ** Method :
1912	!!----------------------------------------------------------------------
1913	INTEGER, INTENT(in) :: ib1, ib2
1914	INTEGER :: itest
1915	!!----------------------------------------------------------------------
1916	itest = 0
1917
1918	IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) ) itest = itest + 1
1919	IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) ) itest = itest + 1
1920	IF( cn_tra (ib1) /= cn_tra (ib2) ) itest = itest + 1
1921	!
1922	IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) ) itest = itest + 1
1923	IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) ) itest = itest + 1
1924	IF( nn_tra_dta (ib1) /= nn_tra_dta (ib2) ) itest = itest + 1
1925	!
1926	IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) ) itest = itest + 1
1927	!
1928	IF( itest>0 ) THEN
1929	WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2
1930	CALL ctl_stop( ctmp1, ' have different open bdy schemes' )
1931	ENDIF
1932	!
1933	END SUBROUTINE bdy_ctl_corn
1934
1935
1936	SUBROUTINE bdy_meshwri()
1937	!!----------------------------------------------------------------------
1938	!! * ROUTINE bdy_meshwri *
1939	!!
1940	!! ** Purpose : write netcdf file with nbr, flagu, flagv, ntreat for T, U
1941	!! and V points in 2D arrays for easier visualisation/control
1942	!!
1943	!! ** Method : use iom_rstput as in domwri.F
1944	!!----------------------------------------------------------------------
1945	INTEGER :: ib_bdy, ii, ij, igrd, ib ! dummy loop indices
1946	INTEGER :: inum ! - -
1947	REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! pointer to 2D mask fields
1948	REAL(wp) , DIMENSION(jpi,jpj) :: ztmp
1949	CHARACTER(LEN=1) , DIMENSION(jpbgrd) :: cgrid
1950	!!----------------------------------------------------------------------
1951	cgrid = (/'t','u','v'/)
1952	CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. )
1953	DO igrd = 1, jpbgrd
1954	SELECT CASE( igrd )
1955	CASE( 1 ) ; zmask => tmask(:,:,1)
1956	CASE( 2 ) ; zmask => umask(:,:,1)
1957	CASE( 3 ) ; zmask => vmask(:,:,1)
1958	END SELECT
1959	ztmp(:,:) = zmask(:,:)
1960	DO ib_bdy = 1, nb_bdy
1961	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) ! nbr deined for all rims
1962	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1963	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1964	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10.
1965	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1966	END DO
1967	END DO
1968	CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1969	ztmp(:,:) = zmask(:,:)
1970	DO ib_bdy = 1, nb_bdy
1971	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! flagu defined only for rims 0 and 1
1972	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1973	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1974	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10.
1975	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1976	END DO
1977	END DO
1978	CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1979	ztmp(:,:) = zmask(:,:)
1980	DO ib_bdy = 1, nb_bdy
1981	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! flagv defined only for rims 0 and 1
1982	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1983	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1984	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10.
1985	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1986	END DO
1987	END DO
1988	CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1989	ztmp(:,:) = zmask(:,:)
1990	DO ib_bdy = 1, nb_bdy
1991	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! ntreat defined only for rims 0 and 1
1992	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1993	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1994	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10.
1995	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1996	END DO
1997	END DO
1998	CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1999	END DO
2000	CALL iom_close( inum )
2001
2002	END SUBROUTINE bdy_meshwri
2003
2004	!!=================================================================================
2005	END MODULE bdyini

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

Download in other formats: