Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

bdyini.F90

Last change on this file was 15368, checked in by smasson, 3 years ago
trunk: final version (hopefully) for ticket #2731
Property svn:keywords set to `Id`
File size: 109.4 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, iiSstdiag, iiSnddiag, iiSstsono, iiSndsono
150	INTEGER :: ijRst, ijRnd, ijSst, ijSnd, ijSstdiag, ijSnddiag, ijSstsono, ijSndsono
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 = Njs0 ! Rcv we-side starting point, excluding sw-corner
810	iiRnd = nn_hls ; ijRnd = Nje0 ! Rcv we-side ending point, excluding nw-corner
811	iiSst = Nie0-nn_hls+1 ; ijSst = Njs0 ! Snd ea-side starting point, excluding se-corner
812	iiSnd = Nie0 ; ijSnd = Nje0 ! 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 = Nie0+1 ; ijRst = Njs0 ! Rcv ea-side starting point, excluding se-corner
820	iiRnd = jpi ; ijRnd = Nje0 ! Rcv ea-side ending point, excluding ne-corner
821	iiSst = Nis0 ; ijSst = Njs0 ! Snd we-side starting point, excluding sw-corner
822	iiSnd = Nis0+nn_hls-1 ; ijSnd = Nje0 ! 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 = Nis0 ; ijRst = 1 ! Rcv so-side starting point, excluding sw-corner
832	iiRnd = Nie0 ; ijRnd = nn_hls ! Rcv so-side ending point, excluding se-corner
833	iiSst = Nis0 ; ijSst = Nje0-nn_hls+1 ! Snd no-side starting point, excluding nw-corner
834	iiSnd = Nie0 ; ijSnd = Nje0 ! 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 = Nis0 ; ijRst = Nje0+1 ! Rcv no-side starting point, excluding nw-corner
843	iiRnd = Nie0 ; ijRnd = jpj ! Rcv no-side ending point, excluding ne-corner
844	iiSst = Nis0 ; ijSst = Njs0 ! Snd so-side starting point, excluding sw-corner
845	iiSnd = Nie0 ; ijSnd = Njs0+nn_hls-1 ! 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	iiRst = 1 ; ijRst = 1 ! Rcv sw-corner starting point
889	iiRnd = nn_hls ; ijRnd = nn_hls ! Rcv sw-corner ending point
890	iiSstdiag = Nie0-nn_hls+1 ; ijSstdiag = Nje0-nn_hls+1 ! send to sw-corner of ne neighbourg
891	iiSnddiag = Nie0 ; ijSnddiag = Nje0 ! send to sw-corner of ne neighbourg
892	iiSstsono = 1 ; ijSstsono = Nje0-nn_hls+1 ! send to sw-corner of no neighbourg
893	iiSndsono = nn_hls ; ijSndsono = Nje0 ! send to sw-corner of no neighbourg
894	iioutdir = -1 ; ijoutdir = -1 ! outside MPI domain: westward or southward
895	! ....\|...
896	CASE( 2 ) ! x: rim on se-corner : \|o
897	! o: potential neighbour(s) __:__x\|o
898	! outside of the MPI domain :o o o
899	! :
900	iRdiag = jpse ; iRsono = jpso ! Recv: for se or so
901	iSdiag = jpnw ; iSsono = jpno ! Send: to nw or no
902	iiRst = Nie0+1 ; ijRst = 1 ! Rcv se-corner starting point
903	iiRnd = jpi ; ijRnd = nn_hls ! Rcv se-corner ending point
904	iiSstdiag = Nis0 ; ijSstdiag = Nje0-nn_hls+1 ! send to se-corner of nw neighbourg
905	iiSnddiag = Nis0+nn_hls-1 ; ijSnddiag = Nje0 ! send to se-corner of nw neighbourg
906	iiSstsono = Nie0+1 ; ijSstsono = Nje0-nn_hls+1 ! send to se-corner of no neighbourg
907	iiSndsono = jpi ; ijSndsono = Nje0 ! send to se-corner of no neighbourg
908	iioutdir = 1 ; ijoutdir = -1 ! outside MPI domain: eastward or southward
909	! :
910	! o o_o:___
911	CASE( 3 ) ! x: rim on nw-corner o\|x :
912	! o: potential neighbour(s) ..o\|...:
913	! outside of the MPI domain \|
914	iRdiag = jpnw ; iRsono = jpno ! Recv: for nw or no
915	iSdiag = jpse ; iSsono = jpso ! Send: to se or so
916	iiRst = 1 ; ijRst = Nje0+1 ! Rcv nw-corner starting point
917	iiRnd = nn_hls ; ijRnd = jpj ! Rcv nw-corner ending point
918	iiSstdiag = Nie0-nn_hls+1 ; ijSstdiag = Njs0 ! send to nw-corner of se neighbourg
919	iiSnddiag = Nie0 ; ijSnddiag = Njs0+nn_hls-1 ! send to nw-corner of se neighbourg
920	iiSstsono = 1 ; ijSstsono = Njs0 ! send to nw-corner of so neighbourg
921	iiSndsono = nn_hls ; ijSndsono = Njs0+nn_hls-1 ! send to nw-corner of so neighbourg
922	iioutdir = -1 ; ijoutdir = 1 ! outside MPI domain: westward or northward
923	! :
924	! ___:o_o o
925	CASE( 4 ) ! x: rim on ne-corner : x\|o
926	! o: potential neighbour(s) :...\|o...
927	! outside of the MPI domain \|
928	iRdiag = jpne ; iRsono = jpno ! Recv: for ne or no
929	iSdiag = jpsw ; iSsono = jpso ! Send: to sw or so
930	iiRst = Nie0+1 ; ijRst = Nje0+1 ! Rcv ne-corner starting point
931	iiRnd = jpi ; ijRnd = jpj ! Rcv ne-corner ending point
932	iiSstdiag = Nis0 ; ijSstdiag = Njs0 ! send to ne-corner of sw neighbourg
933	iiSnddiag = Nis0+nn_hls-1 ; ijSnddiag = Njs0+nn_hls-1 ! send to ne-corner of sw neighbourg
934	iiSstsono = Nie0+1 ; ijSstsono = Njs0 ! send to ne-corner of so neighbourg
935	iiSndsono = jpi ; ijSndsono = Njs0+nn_hls-1 ! send to ne-corner of so neighbourg
936	iioutdir = 1 ; ijoutdir = 1 ! outside MPI domain: eastward or southward
937	END SELECT
938	!
939	! Check if we need to receive data for this rim point
940	IF( ii >= iiRst .AND. ii <= iiRnd .AND. ij >= ijRst .AND. ij <= ijRnd ) THEN ! rim point on the corner
941	iiout = ii+iioutdir ; ijout = ij+ijoutdir ! in which direction do we go outside of the MPI domain?
942	! take care of neighbourg(s) in the interior of the computational domain
943	IF( iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR. & ! Neib outside of the MPI domain
944	& ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) THEN ! -> I cannot compute it -> recv it
945	IF( mpiRnei(nn_hls,iRdiag) > -1 ) lrecv_bdyint(ib_bdy,igrd,iRdiag,ir) = .TRUE. ! Receive directly from diagonal neighbourg
946	IF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 ) lrecv_bdyint(ib_bdy,igrd,iRsono,ir) = .TRUE. ! Receive through the South/North neighbourg
947	ENDIF
948	! take care of neighbourg in the exterior of the computational domain
949	IF( iibe==iiout .OR. ijbe==ijout ) THEN ! Neib outside of the MPI domain -> I cannot compute it -> recv it
950	IF( mpiRnei(nn_hls,iRdiag) > -1 ) lrecv_bdyext(ib_bdy,igrd,iRdiag,ir) = .TRUE. ! Receive directly from diagonal neighbourg
951	IF( mpiRnei(nn_hls,iRsono) > -1 .AND. nn_comm == 1 ) lrecv_bdyext(ib_bdy,igrd,iRsono,ir) = .TRUE. ! Receive through the South/North neighbourg
952	ENDIF
953	ENDIF
954	!
955	! Check if this rim point corresponds to the corner of one neighbourg. if yes, do we need to send data?
956	! Direct send to diag: Is this rim point the corner point of a diag neighbour with which we communicate?
957	IF( ii >= iiSstdiag .AND. ii <= iiSnddiag .AND. ij >= ijSstdiag .AND. ij <= ijSnddiag &
958	& .AND. mpiSnei(nn_hls,iSdiag) > -1 ) THEN
959	iiout = ii+iioutdir ; ijout = ij+ijoutdir ! in which direction do we go outside of the nei MPI domain?
960	! take care of neighbourg(s) in the interior of the computational domain
961	IF( iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR. & ! Neib outside of diag nei MPI
962	& ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) & ! domain -> nei cannot compute it
963	& lsend_bdyint(ib_bdy,igrd,iSdiag,ir) = .TRUE. ! send rim point data to diag nei
964	! take care of neighbourg in the exterior of the computational domain
965	IF( iibe==iiout .OR. ijbe==ijout ) &
966	& lsend_bdyext(ib_bdy,igrd,iSdiag,ir) = .TRUE.
967	ENDIF
968	! Indirect send to diag (through so/no): rim point is the corner point of a so/no nei with which we communicate
969	IF( ii >= iiSstsono .AND. ii <= iiSndsono .AND. ij >= ijSstsono .AND. ij <= ijSndsono &
970	& .AND. mpiSnei(nn_hls,iSsono) > -1 .AND. nn_comm == 1 ) THEN
971	iiout = ii+iioutdir ; ijout = ij+ijoutdir ! in which direction do we go outside of the nei MPI domain?
972	! take care of neighbourg(s) in the interior of the computational domain
973	IF( iibi==iiout .OR. ii1==iiout .OR. ii2==iiout .OR. ii3==iiout .OR. & ! Neib outside of so/no nei MPI
974	& ijbi==ijout .OR. ij1==ijout .OR. ij2==ijout .OR. ij3==ijout ) & ! domain -> nei cannot compute it
975	& lsend_bdyint(ib_bdy,igrd,iSsono,ir) = .TRUE. ! send rim point data to so/no nei
976	! take care of neighbourg in the exterior of the computational domain
977	IF( iibe==iiout .OR. ijbe==ijout ) &
978	& lsend_bdyext(ib_bdy,igrd,iSsono,ir) = .TRUE.
979	ENDIF
980	!
981	END DO ! 4 corners
982	END DO ! ib
983	END DO ! igrd
984
985	! Comment out for debug
986	!!$ DO ir = 0,1
987	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing, &
988	!!$ & lsend = lsend_bdyint(ib_bdy,1,:,ir), lrecv = lrecv_bdyint(ib_bdy,1,:,ir) )
989	!!$ IF(lwp) WRITE(numout,*) ' bdy debug int T', ir ; CALL FLUSH(numout)
990	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing, &
991	!!$ & lsend = lsend_bdyint(ib_bdy,2,:,ir), lrecv = lrecv_bdyint(ib_bdy,2,:,ir) )
992	!!$ IF(lwp) WRITE(numout,*) ' bdy debug int U', ir ; CALL FLUSH(numout)
993	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing, &
994	!!$ & lsend = lsend_bdyint(ib_bdy,3,:,ir), lrecv = lrecv_bdyint(ib_bdy,3,:,ir) )
995	!!$ IF(lwp) WRITE(numout,*) ' bdy debug int V', ir ; CALL FLUSH(numout)
996	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'T', 1._wp, kfillmode = jpfillnothing, &
997	!!$ & lsend = lsend_bdyext(ib_bdy,1,:,ir), lrecv = lrecv_bdyext(ib_bdy,1,:,ir) )
998	!!$ IF(lwp) WRITE(numout,*) ' bdy debug ext T', ir ; CALL FLUSH(numout)
999	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'U', 1._wp, kfillmode = jpfillnothing, &
1000	!!$ & lsend = lsend_bdyext(ib_bdy,2,:,ir), lrecv = lrecv_bdyext(ib_bdy,2,:,ir) )
1001	!!$ IF(lwp) WRITE(numout,*) ' bdy debug ext U', ir ; CALL FLUSH(numout)
1002	!!$ zzbdy(:,:) = narea ; CALL lbc_lnk('bdy debug', zzbdy, 'V', 1._wp, kfillmode = jpfillnothing, &
1003	!!$ & lsend = lsend_bdyext(ib_bdy,3,:,ir), lrecv = lrecv_bdyext(ib_bdy,3,:,ir) )
1004	!!$ IF(lwp) WRITE(numout,*) ' bdy debug ext V', ir ; CALL FLUSH(numout)
1005	!!$ END DO
1006
1007	END DO ! ib_bdy
1008
1009	DO ib_bdy = 1,nb_bdy
1010	IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &
1011	& cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &
1012	& cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo' ) THEN
1013	DO igrd = 1, jpbgrd
1014	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1015	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1016	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1017	IF( mig0(ii) > 2 .AND. mig0(ii) < Ni0glo-2 .AND. mjg0(ij) > 2 .AND. mjg0(ij) < Nj0glo-2 ) THEN
1018	WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain'
1019	CALL ctl_stop( ctmp1 )
1020	END IF
1021	END DO
1022	END DO
1023	END IF
1024	END DO
1025	!
1026	DEALLOCATE( nbidta, nbjdta, nbrdta )
1027	!
1028	END SUBROUTINE bdy_def
1029
1030
1031	SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 )
1032	!!----------------------------------------------------------------------
1033	!! * ROUTINE bdy_rim_treat *
1034	!!
1035	!! ** Purpose : Initialize structures ( flagu, flagv, ntreat ) indicating how rim points
1036	!! are to be handled in the boundary condition treatment
1037	!!
1038	!! ** Method : - to handle rim 0 zmasks must indicate ocean points (set at one on rim 0 and rim 1 and interior)
1039	!! and bdymasks must be set at 0 on rim 0 (set at one on rim 1 and interior)
1040	!! (as if rim 1 was free ocean)
1041	!! - to handle rim 1 zmasks must be set at 0 on rim 0 (set at one on rim 1 and interior)
1042	!! and bdymasks must indicate free ocean points (set at one on interior)
1043	!! (as if rim 0 was land)
1044	!! - we can then check in which direction the interior of the computational domain is with the difference
1045	!! mask array values on both sides to compute flagu and flagv
1046	!! - and look at the ocean neighbours to compute ntreat
1047	!!----------------------------------------------------------------------
1048	REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: pumask, pvmask ! temporary u/v mask array
1049	REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: pfmask ! temporary fmask excluding coastal boundary condition (shlat)
1050	LOGICAL , INTENT (in ) :: lrim0 ! .true. -> rim 0 .false. -> rim 1
1051	INTEGER :: ib_bdy, ii, ij, igrd, ib, icount ! dummy loop indices
1052	INTEGER :: i_offset, j_offset, inn ! local integer
1053	INTEGER :: ibeg, iend ! local integer
1054	LOGICAL :: llnon, llson, llean, llwen ! local logicals indicating the presence of a ocean neighbour
1055	REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! pointer to 2D mask fields
1056	REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars
1057	CHARACTER(LEN=1), DIMENSION(jpbgrd) :: cgrid
1058	REAL(wp) , DIMENSION(jpi,jpj) :: ztmp
1059	!!----------------------------------------------------------------------
1060
1061	cgrid = (/'t','u','v'/)
1062
1063	DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components
1064
1065	DO igrd = 1, jpbgrd
1066
1067	IF( lrim0 ) THEN ! extent of rim 0
1068	ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
1069	ELSE ! extent of rim 1
1070	ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1 ; iend = idx_bdy(ib_bdy)%nblenrim(igrd)
1071	END IF
1072
1073	! Calculate relationship of U direction to the local orientation of the boundary
1074	! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
1075	! flagu = 0 : u is tangential
1076	! flagu = 1 : u is normal to the boundary and is direction is inward
1077	SELECT CASE( igrd )
1078	CASE( 1 ) ; zmask => pumask ; i_offset = 0 ! U(i-1) T(i) U(i )
1079	CASE( 2 ) ; zmask => bdytmask ; i_offset = 1 ! T(i ) U(i) T(i+1)
1080	CASE( 3 ) ; zmask => pfmask ; i_offset = 0 ! F(i-1) V(i) F(i )
1081	END SELECT
1082	icount = 0
1083	ztmp(:,:) = -999._wp
1084	DO ib = ibeg, iend
1085	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1086	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1087	IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 ) CYCLE ! call lbc_lnk -> no need to compute these pts
1088	zwfl = zmask(ii+i_offset-1,ij)
1089	zefl = zmask(ii+i_offset ,ij)
1090	! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
1091	IF( i_offset == 1 .and. zefl + zwfl == 2._wp ) THEN
1092	icount = icount + 1
1093	IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
1094	ELSE
1095	ztmp(ii,ij) = -zwfl + zefl
1096	ENDIF
1097	END DO
1098	IF( icount /= 0 ) THEN
1099	WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,', &
1100	' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1101	CALL ctl_stop( ctmp1 )
1102	ENDIF
1103	SELECT CASE( igrd )
1104	CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1105	CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1106	CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1107	END SELECT
1108	DO ib = ibeg, iend
1109	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1110	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1111	idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij)
1112	END DO
1113
1114	! Calculate relationship of V direction to the local orientation of the boundary
1115	! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
1116	! flagv = 0 : v is tangential
1117	! flagv = 1 : v is normal to the boundary and is direction is inward
1118	SELECT CASE( igrd )
1119	CASE( 1 ) ; zmask => pvmask ; j_offset = 0
1120	CASE( 2 ) ; zmask => pfmask ; j_offset = 0
1121	CASE( 3 ) ; zmask => bdytmask ; j_offset = 1
1122	END SELECT
1123	icount = 0
1124	ztmp(:,:) = -999._wp
1125	DO ib = ibeg, iend
1126	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1127	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1128	IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 ) CYCLE ! call lbc_lnk -> no need to compute these pts
1129	zsfl = zmask(ii,ij+j_offset-1)
1130	znfl = zmask(ii,ij+j_offset )
1131	! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
1132	IF( j_offset == 1 .and. znfl + zsfl == 2._wp ) THEN
1133	IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
1134	icount = icount + 1
1135	ELSE
1136	ztmp(ii,ij) = -zsfl + znfl
1137	END IF
1138	END DO
1139	IF( icount /= 0 ) THEN
1140	WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,', &
1141	' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
1142	CALL ctl_stop( ctmp1 )
1143	ENDIF
1144	SELECT CASE( igrd )
1145	CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1146	CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1147	CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1148	END SELECT
1149	DO ib = ibeg, iend
1150	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1151	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1152	idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij)
1153	END DO
1154
1155	! Calculate ntreat
1156	SELECT CASE( igrd )
1157	CASE( 1 ) ; zmask => bdytmask
1158	CASE( 2 ) ; zmask => bdyumask
1159	CASE( 3 ) ; zmask => bdyvmask
1160	END SELECT
1161	ztmp(:,:) = -999._wp
1162	DO ib = ibeg, iend
1163	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1164	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1165	IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 ) CYCLE ! call lbc_lnk -> no need to compute these pts
1166	llnon = zmask(ii ,ij+1) == 1._wp
1167	llson = zmask(ii ,ij-1) == 1._wp
1168	llean = zmask(ii+1,ij ) == 1._wp
1169	llwen = zmask(ii-1,ij ) == 1._wp
1170	inn = COUNT( (/ llnon, llson, llean, llwen /) )
1171	IF( inn == 0 ) THEN ! no neighbours -> interior of a corner or cluster of rim points
1172	! ! ! _____ ! _____ ! __ __
1173	! 1 \| o ! 2 o \| ! 3 \| x ! 4 x \| ! \| \| -> error
1174	! \|_x_ _ ! _ _x_\| ! \| o ! o \| ! \|x_x\|
1175	IF( zmask(ii+1,ij+1) == 1._wp ) THEN ; ztmp(ii,ij) = 1._wp
1176	ELSEIF( zmask(ii-1,ij+1) == 1._wp ) THEN ; ztmp(ii,ij) = 2._wp
1177	ELSEIF( zmask(ii+1,ij-1) == 1._wp ) THEN ; ztmp(ii,ij) = 3._wp
1178	ELSEIF( zmask(ii-1,ij-1) == 1._wp ) THEN ; ztmp(ii,ij) = 4._wp
1179	ELSE ; ztmp(ii,ij) = -1._wp
1180	WRITE(ctmp1,*) 'Problem with ',cgrid(igrd) ,' grid point', ii, ij, &
1181	' on boundary set ', ib_bdy, ' has no free ocean neighbour'
1182	IF( lrim0 ) THEN
1183	WRITE(ctmp2,*) ' There seems to be a cluster of rim 0 points.'
1184	ELSE
1185	WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.'
1186	END IF
1187	CALL ctl_warn( ctmp1, ctmp2 )
1188	END IF
1189	END IF
1190	IF( inn == 1 ) THEN ! middle of linear bdy or incomplete corner ! ___ o
1191	! \| ! \| ! o ! ______ ! \|x___
1192	! 5 \| x o ! 6 o x \| ! 7 __x__ ! 8 x
1193	! \| ! \| ! ! o
1194	IF( llean ) ztmp(ii,ij) = 5._wp
1195	IF( llwen ) ztmp(ii,ij) = 6._wp
1196	IF( llnon ) ztmp(ii,ij) = 7._wp
1197	IF( llson ) ztmp(ii,ij) = 8._wp
1198	END IF
1199	IF( inn == 2 ) THEN ! exterior of a corner
1200	! o ! o ! _____\| ! \|_____
1201	! 9 ____x o ! 10 o x___ ! 11 x o ! 12 o x
1202	! \| ! \| ! o ! o
1203	IF( llnon .AND. llean ) ztmp(ii,ij) = 9._wp
1204	IF( llnon .AND. llwen ) ztmp(ii,ij) = 10._wp
1205	IF( llson .AND. llean ) ztmp(ii,ij) = 11._wp
1206	IF( llson .AND. llwen ) ztmp(ii,ij) = 12._wp
1207	END IF
1208	IF( inn == 3 ) THEN ! 3 neighbours __ __
1209	! \|_ o ! o _\| ! \|_\| ! o
1210	! 13 _\| x o ! 14 o x \|_ ! 15 o x o ! 16 o x o
1211	! \| o ! o \| ! o ! __\|¨\|__
1212	IF( llnon .AND. llean .AND. llson ) ztmp(ii,ij) = 13._wp
1213	IF( llnon .AND. llwen .AND. llson ) ztmp(ii,ij) = 14._wp
1214	IF( llwen .AND. llson .AND. llean ) ztmp(ii,ij) = 15._wp
1215	IF( llwen .AND. llnon .AND. llean ) ztmp(ii,ij) = 16._wp
1216	END IF
1217	IF( inn == 4 ) THEN
1218	WRITE(ctmp1,*) 'Problem with ',cgrid(igrd) ,' grid point', ii, ij, &
1219	' on boundary set ', ib_bdy, ' have 4 neighbours'
1220	CALL ctl_stop( ctmp1 )
1221	END IF
1222	END DO
1223	SELECT CASE( igrd )
1224	CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1.0_wp )
1225	CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1.0_wp )
1226	CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1.0_wp )
1227	END SELECT
1228	DO ib = ibeg, iend
1229	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1230	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1231	idx_bdy(ib_bdy)%ntreat(ib,igrd) = NINT(ztmp(ii,ij))
1232	END DO
1233	!
1234	END DO ! jpbgrd
1235	!
1236	END DO ! ib_bdy
1237
1238	END SUBROUTINE bdy_rim_treat
1239
1240
1241	SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )
1242	!!----------------------------------------------------------------------
1243	!! * ROUTINE find_neib *
1244	!!
1245	!! ** Purpose : get ii1, ij1, ii2, ij2, ii3, ij3, the indices of
1246	!! the free ocean neighbours of (ii,ij) for bdy treatment
1247	!!
1248	!! ** Method : use itreat input to select a case
1249	!! N.B. ntreat is defined for all bdy points in routine bdy_rim_treat
1250	!!
1251	!!----------------------------------------------------------------------
1252	INTEGER, INTENT(in ) :: ii, ij, itreat
1253	INTEGER, INTENT( out) :: ii1, ij1, ii2, ij2, ii3, ij3
1254	!!----------------------------------------------------------------------
1255	SELECT CASE( itreat ) ! points that will be used by bdy routines, -1 will be discarded
1256	! ! ! _____ ! _____
1257	! 1 \| o ! 2 o \| ! 3 \| x ! 4 x \|
1258	! \|_x_ _ ! _ _x_\| ! \| o ! o \|
1259	CASE( 1 ) ; ii1 = ii+1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1260	CASE( 2 ) ; ii1 = ii-1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1261	CASE( 3 ) ; ii1 = ii+1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1262	CASE( 4 ) ; ii1 = ii-1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1263	! \| ! \| ! o ! ______ ! or incomplete corner
1264	! 5 \| x o ! 6 o x \| ! 7 __x__ ! 8 x ! 7 ____ o
1265	! \| ! \| ! ! o ! \|x___
1266	CASE( 5 ) ; ii1 = ii+1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1267	CASE( 6 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1268	CASE( 7 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1269	CASE( 8 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1270	! o ! o ! _____\| ! \|_____
1271	! 9 ____x o ! 10 o x___ ! 11 x o ! 12 o x
1272	! \| ! \| ! o ! o
1273	CASE( 9 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1274	CASE( 10 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1275	CASE( 11 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1276	CASE( 12 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1277	! \|_ o ! o _\| ! ¨¨\|_\|¨¨ ! o
1278	! 13 _\| x o ! 14 o x \|_ ! 15 o x o ! 16 o x o
1279	! \| o ! o \| ! o ! __\|¨\|__
1280	CASE( 13 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1
1281	CASE( 14 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1
1282	CASE( 15 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij-1 ; ii3 = ii+1 ; ij3 = ij
1283	CASE( 16 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij+1 ; ii3 = ii+1 ; ij3 = ij
1284	END SELECT
1285	END SUBROUTINE find_neib
1286
1287
1288	SUBROUTINE bdy_read_seg( kb_bdy, knblendta )
1289	!!----------------------------------------------------------------------
1290	!! * ROUTINE bdy_read_seg *
1291	!!
1292	!! ** Purpose : build bdy coordinates with segments defined in namelist
1293	!!
1294	!! ** Method : read namelist nambdy_index blocks
1295	!!
1296	!!----------------------------------------------------------------------
1297	INTEGER , INTENT (in ) :: kb_bdy ! bdy number
1298	INTEGER, DIMENSION(jpbgrd), INTENT ( out) :: knblendta ! length of index arrays
1299	!!
1300	INTEGER :: ios ! Local integer output status for namelist read
1301	INTEGER :: nbdyind, nbdybeg, nbdyend
1302	INTEGER :: nbdy_count, nbdy_rdstart, nbdy_loc
1303	CHARACTER(LEN=1) :: ctypebdy ! - -
1304	CHARACTER(LEN=50):: cerrmsg ! - -
1305	NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
1306	!!----------------------------------------------------------------------
1307	! Need to support possibility of reading more than one nambdy_index from
1308	! the namelist_cfg internal file.
1309	! Do this by finding the kb_bdy'th occurence of nambdy_index in the
1310	! character buffer as the starting point.
1311	nbdy_rdstart = 1
1312	DO nbdy_count = 1, kb_bdy
1313	nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_index' )
1314	IF( nbdy_loc .GT. 0 ) THEN
1315	nbdy_rdstart = nbdy_rdstart + nbdy_loc
1316	ELSE
1317	WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',kb_bdy,' of nambdy_index not found'
1318	ios = -1
1319	CALL ctl_nam ( ios , cerrmsg )
1320	ENDIF
1321	END DO
1322	nbdy_rdstart = MAX( 1, nbdy_rdstart - 2 )
1323	READ ( numnam_cfg( nbdy_rdstart: ), nambdy_index, IOSTAT = ios, ERR = 904)
1324	904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' )
1325	IF(lwm) WRITE ( numond, nambdy_index )
1326
1327	SELECT CASE ( TRIM(ctypebdy) )
1328	CASE( 'N' )
1329	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1330	nbdyind = Nj0glo - 2 ! set boundary to whole side of model domain.
1331	nbdybeg = 2
1332	nbdyend = Ni0glo - 1
1333	ENDIF
1334	nbdysegn = nbdysegn + 1
1335	npckgn(nbdysegn) = kb_bdy ! Save bdy package number
1336	jpjnob(nbdysegn) = nbdyind
1337	jpindt(nbdysegn) = nbdybeg
1338	jpinft(nbdysegn) = nbdyend
1339	!
1340	CASE( 'S' )
1341	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1342	nbdyind = 2 ! set boundary to whole side of model domain.
1343	nbdybeg = 2
1344	nbdyend = Ni0glo - 1
1345	ENDIF
1346	nbdysegs = nbdysegs + 1
1347	npckgs(nbdysegs) = kb_bdy ! Save bdy package number
1348	jpjsob(nbdysegs) = nbdyind
1349	jpisdt(nbdysegs) = nbdybeg
1350	jpisft(nbdysegs) = nbdyend
1351	!
1352	CASE( 'E' )
1353	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1354	nbdyind = Ni0glo - 2 ! set boundary to whole side of model domain.
1355	nbdybeg = 2
1356	nbdyend = Nj0glo - 1
1357	ENDIF
1358	nbdysege = nbdysege + 1
1359	npckge(nbdysege) = kb_bdy ! Save bdy package number
1360	jpieob(nbdysege) = nbdyind
1361	jpjedt(nbdysege) = nbdybeg
1362	jpjeft(nbdysege) = nbdyend
1363	!
1364	CASE( 'W' )
1365	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1366	nbdyind = 2 ! set boundary to whole side of model domain.
1367	nbdybeg = 2
1368	nbdyend = Nj0glo - 1
1369	ENDIF
1370	nbdysegw = nbdysegw + 1
1371	npckgw(nbdysegw) = kb_bdy ! Save bdy package number
1372	jpiwob(nbdysegw) = nbdyind
1373	jpjwdt(nbdysegw) = nbdybeg
1374	jpjwft(nbdysegw) = nbdyend
1375	!
1376	CASE DEFAULT ; CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
1377	END SELECT
1378
1379	! For simplicity we assume that in case of straight bdy, arrays have the same length
1380	! (even if it is true that last tangential velocity points
1381	! are useless). This simplifies a little bit boundary data format (and agrees with format
1382	! used so far in obc package)
1383
1384	knblendta(1:jpbgrd) = (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy)
1385
1386	END SUBROUTINE bdy_read_seg
1387
1388
1389	SUBROUTINE bdy_ctl_seg
1390	!!----------------------------------------------------------------------
1391	!! * ROUTINE bdy_ctl_seg *
1392	!!
1393	!! ** Purpose : Check straight open boundary segments location
1394	!!
1395	!! ** Method : - Look for open boundary corners
1396	!! - Check that segments start or end on land
1397	!!----------------------------------------------------------------------
1398	INTEGER :: ib, ib1, ib2, ji ,jj, itest
1399	INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns
1400	REAL(wp), DIMENSION(2) :: ztestmask
1401	!!----------------------------------------------------------------------
1402	!
1403	IF (lwp) WRITE(numout,*) ' '
1404	IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1405	IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1406	!
1407	IF(lwp) WRITE(numout,*) 'Number of east segments : ', nbdysege
1408	IF(lwp) WRITE(numout,*) 'Number of west segments : ', nbdysegw
1409	IF(lwp) WRITE(numout,*) 'Number of north segments : ', nbdysegn
1410	IF(lwp) WRITE(numout,*) 'Number of south segments : ', nbdysegs
1411	!
1412	! 1. Check bounds
1413	!----------------
1414	DO ib = 1, nbdysegn
1415	IF (lwp) WRITE(numout,) '*check north seg bounds pckg: ', npckgn(ib)
1416	IF ((jpjnob(ib).ge.Nj0glo-1).or.&
1417	&(jpjnob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1418	IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1419	IF (jpindt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1420	IF (jpinft(ib).gt.Ni0glo) CALL ctl_stop( 'End index out of domain' )
1421	END DO
1422	!
1423	DO ib = 1, nbdysegs
1424	IF (lwp) WRITE(numout,) '*check south seg bounds pckg: ', npckgs(ib)
1425	IF ((jpjsob(ib).ge.Nj0glo-1).or.&
1426	&(jpjsob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1427	IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1428	IF (jpisdt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1429	IF (jpisft(ib).gt.Ni0glo) CALL ctl_stop( 'End index out of domain' )
1430	END DO
1431	!
1432	DO ib = 1, nbdysege
1433	IF (lwp) WRITE(numout,) '*check east seg bounds pckg: ', npckge(ib)
1434	IF ((jpieob(ib).ge.Ni0glo-1).or.&
1435	&(jpieob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1436	IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1437	IF (jpjedt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1438	IF (jpjeft(ib).gt.Nj0glo) CALL ctl_stop( 'End index out of domain' )
1439	END DO
1440	!
1441	DO ib = 1, nbdysegw
1442	IF (lwp) WRITE(numout,) '*check west seg bounds pckg: ', npckgw(ib)
1443	IF ((jpiwob(ib).ge.Ni0glo-1).or.&
1444	&(jpiwob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1445	IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1446	IF (jpjwdt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1447	IF (jpjwft(ib).gt.Nj0glo) CALL ctl_stop( 'End index out of domain' )
1448	ENDDO
1449	!
1450	! 2. Look for segment crossings
1451	!------------------------------
1452	IF (lwp) WRITE(numout,) '*Look for segments corners :'
1453	!
1454	itest = 0 ! corner number
1455	!
1456	! flag to detect if start or end of open boundary belongs to a corner
1457	! if not (=0), it must be on land.
1458	! if a corner is detected, save bdy package number for further tests
1459	icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1460	! South/West crossings
1461	IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1462	DO ib1 = 1, nbdysegw
1463	DO ib2 = 1, nbdysegs
1464	IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1465	& ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1466	& ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1467	& ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1468	IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN
1469	! We have a possible South-West corner
1470	! WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1471	! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1472	icornw(ib1,1) = npckgs(ib2)
1473	icorns(ib2,1) = npckgw(ib1)
1474	ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1475	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1476	& jpisft(ib2), jpjwft(ib1)
1477	WRITE(ctmp2,*) ' Not allowed yet'
1478	WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1479	& ' and South segment: ',npckgs(ib2)
1480	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1481	ELSE
1482	WRITE(ctmp1,*) ' Check South and West Open boundary indices'
1483	WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , &
1484	& ' and South segment: ',npckgs(ib2)
1485	CALL ctl_stop( ctmp1, ctmp2 )
1486	END IF
1487	END IF
1488	END DO
1489	END DO
1490	END IF
1491	!
1492	! South/East crossings
1493	IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1494	DO ib1 = 1, nbdysege
1495	DO ib2 = 1, nbdysegs
1496	IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1497	& ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1498	& ( jpjedt(ib1)<=jpjsob(ib2) ).AND. &
1499	& ( jpjeft(ib1)>=jpjsob(ib2) )) THEN
1500	IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1501	! We have a possible South-East corner
1502	! WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1503	! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1504	icorne(ib1,1) = npckgs(ib2)
1505	icorns(ib2,2) = npckge(ib1)
1506	ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1507	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1508	& jpisdt(ib2), jpjeft(ib1)
1509	WRITE(ctmp2,*) ' Not allowed yet'
1510	WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1511	& ' and South segment: ',npckgs(ib2)
1512	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1513	ELSE
1514	WRITE(ctmp1,*) ' Check South and East Open boundary indices'
1515	WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1516	& ' and South segment: ',npckgs(ib2)
1517	CALL ctl_stop( ctmp1, ctmp2 )
1518	END IF
1519	END IF
1520	END DO
1521	END DO
1522	END IF
1523	!
1524	! North/West crossings
1525	IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1526	DO ib1 = 1, nbdysegw
1527	DO ib2 = 1, nbdysegn
1528	IF (( jpindt(ib2)<=jpiwob(ib1) ).AND. &
1529	& ( jpinft(ib2)>=jpiwob(ib1) ).AND. &
1530	& ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1531	& ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1532	IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1533	! We have a possible North-West corner
1534	! WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1535	! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1536	icornw(ib1,2) = npckgn(ib2)
1537	icornn(ib2,1) = npckgw(ib1)
1538	ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1539	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1540	& jpinft(ib2), jpjwdt(ib1)
1541	WRITE(ctmp2,*) ' Not allowed yet'
1542	WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1543	& ' and North segment: ',npckgn(ib2)
1544	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1545	ELSE
1546	WRITE(ctmp1,*) ' Check North and West Open boundary indices'
1547	WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1548	& ' and North segment: ',npckgn(ib2)
1549	CALL ctl_stop( ctmp1, ctmp2 )
1550	END IF
1551	END IF
1552	END DO
1553	END DO
1554	END IF
1555	!
1556	! North/East crossings
1557	IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1558	DO ib1 = 1, nbdysege
1559	DO ib2 = 1, nbdysegn
1560	IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1561	& ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1562	& ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1563	& ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1564	IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1565	! We have a possible North-East corner
1566	! WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1567	! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1568	icorne(ib1,2) = npckgn(ib2)
1569	icornn(ib2,2) = npckge(ib1)
1570	ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1571	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1572	& jpindt(ib2), jpjedt(ib1)
1573	WRITE(ctmp2,*) ' Not allowed yet'
1574	WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1575	& ' and North segment: ',npckgn(ib2)
1576	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1577	ELSE
1578	WRITE(ctmp1,*) ' Check North and East Open boundary indices'
1579	WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1580	& ' and North segment: ',npckgn(ib2)
1581	CALL ctl_stop( ctmp1, ctmp2 )
1582	END IF
1583	END IF
1584	END DO
1585	END DO
1586	END IF
1587	!
1588	! 3. Check if segment extremities are on land
1589	!--------------------------------------------
1590	!
1591	! West segments
1592	DO ib = 1, nbdysegw
1593	! get mask at boundary extremities:
1594	ztestmask(1:2)=0.
1595	DO ji = 1, jpi
1596	DO jj = 1, jpj
1597	IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwdt(ib) ) ztestmask(1) = tmask(ji,jj,1)
1598	IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwft(ib) ) ztestmask(2) = tmask(ji,jj,1)
1599	END DO
1600	END DO
1601	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1602
1603	IF (ztestmask(1)==1) THEN
1604	IF (icornw(ib,1)==0) THEN
1605	WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1606	CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1607	ELSE
1608	! This is a corner
1609	IF(lwp) WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1610	CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1611	itest=itest+1
1612	ENDIF
1613	ENDIF
1614	IF (ztestmask(2)==1) THEN
1615	IF (icornw(ib,2)==0) THEN
1616	WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1617	CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' )
1618	ELSE
1619	! This is a corner
1620	IF(lwp) WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1621	CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1622	itest=itest+1
1623	ENDIF
1624	ENDIF
1625	END DO
1626	!
1627	! East segments
1628	DO ib = 1, nbdysege
1629	! get mask at boundary extremities:
1630	ztestmask(1:2)=0.
1631	DO ji = 1, jpi
1632	DO jj = 1, jpj
1633	IF( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjedt(ib) ) ztestmask(1) = tmask(ji,jj,1)
1634	IF( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjeft(ib) ) ztestmask(2) = tmask(ji,jj,1)
1635	END DO
1636	END DO
1637	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1638
1639	IF (ztestmask(1)==1) THEN
1640	IF (icorne(ib,1)==0) THEN
1641	WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1642	CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1643	ELSE
1644	! This is a corner
1645	IF(lwp) WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1646	CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1647	itest=itest+1
1648	ENDIF
1649	ENDIF
1650	IF (ztestmask(2)==1) THEN
1651	IF (icorne(ib,2)==0) THEN
1652	WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1653	CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1654	ELSE
1655	! This is a corner
1656	IF(lwp) WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1657	CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1658	itest=itest+1
1659	ENDIF
1660	ENDIF
1661	END DO
1662	!
1663	! South segments
1664	DO ib = 1, nbdysegs
1665	! get mask at boundary extremities:
1666	ztestmask(1:2)=0.
1667	DO ji = 1, jpi
1668	DO jj = 1, jpj
1669	IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisdt(ib) ) ztestmask(1) = tmask(ji,jj,1)
1670	IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisft(ib) ) ztestmask(2) = tmask(ji,jj,1)
1671	END DO
1672	END DO
1673	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1674
1675	IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1676	WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1677	CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1678	ENDIF
1679	IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1680	WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1681	CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1682	ENDIF
1683	END DO
1684	!
1685	! North segments
1686	DO ib = 1, nbdysegn
1687	! get mask at boundary extremities:
1688	ztestmask(1:2)=0.
1689	DO ji = 1, jpi
1690	DO jj = 1, jpj
1691	IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpindt(ib) ) ztestmask(1) = tmask(ji,jj,1)
1692	IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpinft(ib) ) ztestmask(2) = tmask(ji,jj,1)
1693	END DO
1694	END DO
1695	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1696
1697	IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1698	WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1699	CALL ctl_stop( ctmp1, ' does not start on land' )
1700	ENDIF
1701	IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1702	WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1703	CALL ctl_stop( ctmp1, ' does not end on land' )
1704	ENDIF
1705	END DO
1706	!
1707	IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1708	!
1709	! Other tests TBD:
1710	! segments completly on land
1711	! optimized open boundary array length according to landmask
1712	! Nudging layers that overlap with interior domain
1713	!
1714	END SUBROUTINE bdy_ctl_seg
1715
1716
1717	SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta )
1718	!!----------------------------------------------------------------------
1719	!! * ROUTINE bdy_coords_seg *
1720	!!
1721	!! ** Purpose : build nbidta, nbidta, nbrdta for bdy built with segments
1722	!!
1723	!! ** Method :
1724	!!
1725	!!----------------------------------------------------------------------
1726	INTEGER, DIMENSION(:,:,:), intent( out) :: nbidta, nbjdta, nbrdta ! Index arrays: i and j indices of bdy dta
1727	!!
1728	INTEGER :: ii, ij, ir, iseg
1729	INTEGER :: igrd ! grid type (t=1, u=2, v=3)
1730	INTEGER :: icount !
1731	INTEGER :: ib_bdy ! bdy number
1732	!!----------------------------------------------------------------------
1733
1734	! East
1735	!-----
1736	DO iseg = 1, nbdysege
1737	ib_bdy = npckge(iseg)
1738	!
1739	! ------------ T points -------------
1740	igrd=1
1741	icount=0
1742	DO ir = 1, nn_rimwidth(ib_bdy)
1743	DO ij = jpjedt(iseg), jpjeft(iseg)
1744	icount = icount + 1
1745	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir + nn_hls
1746	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1747	nbrdta(icount, igrd, ib_bdy) = ir
1748	ENDDO
1749	ENDDO
1750	!
1751	! ------------ U points -------------
1752	igrd=2
1753	icount=0
1754	DO ir = 1, nn_rimwidth(ib_bdy)
1755	DO ij = jpjedt(iseg), jpjeft(iseg)
1756	icount = icount + 1
1757	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir + nn_hls
1758	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1759	nbrdta(icount, igrd, ib_bdy) = ir
1760	ENDDO
1761	ENDDO
1762	!
1763	! ------------ V points -------------
1764	igrd=3
1765	icount=0
1766	DO ir = 1, nn_rimwidth(ib_bdy)
1767	! DO ij = jpjedt(iseg), jpjeft(iseg) - 1
1768	DO ij = jpjedt(iseg), jpjeft(iseg)
1769	icount = icount + 1
1770	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir + nn_hls
1771	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1772	nbrdta(icount, igrd, ib_bdy) = ir
1773	ENDDO
1774	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1775	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1776	ENDDO
1777	ENDDO
1778	!
1779	! West
1780	!-----
1781	DO iseg = 1, nbdysegw
1782	ib_bdy = npckgw(iseg)
1783	!
1784	! ------------ T points -------------
1785	igrd=1
1786	icount=0
1787	DO ir = 1, nn_rimwidth(ib_bdy)
1788	DO ij = jpjwdt(iseg), jpjwft(iseg)
1789	icount = icount + 1
1790	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1791	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1792	nbrdta(icount, igrd, ib_bdy) = ir
1793	ENDDO
1794	ENDDO
1795	!
1796	! ------------ U points -------------
1797	igrd=2
1798	icount=0
1799	DO ir = 1, nn_rimwidth(ib_bdy)
1800	DO ij = jpjwdt(iseg), jpjwft(iseg)
1801	icount = icount + 1
1802	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1803	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1804	nbrdta(icount, igrd, ib_bdy) = ir
1805	ENDDO
1806	ENDDO
1807	!
1808	! ------------ V points -------------
1809	igrd=3
1810	icount=0
1811	DO ir = 1, nn_rimwidth(ib_bdy)
1812	! DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
1813	DO ij = jpjwdt(iseg), jpjwft(iseg)
1814	icount = icount + 1
1815	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 + nn_hls
1816	nbjdta(icount, igrd, ib_bdy) = ij + nn_hls
1817	nbrdta(icount, igrd, ib_bdy) = ir
1818	ENDDO
1819	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1820	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1821	ENDDO
1822	ENDDO
1823	!
1824	! North
1825	!-----
1826	DO iseg = 1, nbdysegn
1827	ib_bdy = npckgn(iseg)
1828	!
1829	! ------------ T points -------------
1830	igrd=1
1831	icount=0
1832	DO ir = 1, nn_rimwidth(ib_bdy)
1833	DO ii = jpindt(iseg), jpinft(iseg)
1834	icount = icount + 1
1835	nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1836	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir + nn_hls
1837	nbrdta(icount, igrd, ib_bdy) = ir
1838	ENDDO
1839	ENDDO
1840	!
1841	! ------------ U points -------------
1842	igrd=2
1843	icount=0
1844	DO ir = 1, nn_rimwidth(ib_bdy)
1845	! DO ii = jpindt(iseg), jpinft(iseg) - 1
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) + 2 - ir + nn_hls
1850	nbrdta(icount, igrd, ib_bdy) = ir
1851	ENDDO
1852	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1853	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1854	ENDDO
1855	!
1856	! ------------ V points -------------
1857	igrd=3
1858	icount=0
1859	DO ir = 1, nn_rimwidth(ib_bdy)
1860	DO ii = jpindt(iseg), jpinft(iseg)
1861	icount = icount + 1
1862	nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1863	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir + nn_hls
1864	nbrdta(icount, igrd, ib_bdy) = ir
1865	ENDDO
1866	ENDDO
1867	ENDDO
1868	!
1869	! South
1870	!-----
1871	DO iseg = 1, nbdysegs
1872	ib_bdy = npckgs(iseg)
1873	!
1874	! ------------ T points -------------
1875	igrd=1
1876	icount=0
1877	DO ir = 1, nn_rimwidth(ib_bdy)
1878	DO ii = jpisdt(iseg), jpisft(iseg)
1879	icount = icount + 1
1880	nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1881	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1882	nbrdta(icount, igrd, ib_bdy) = ir
1883	ENDDO
1884	ENDDO
1885	!
1886	! ------------ U points -------------
1887	igrd=2
1888	icount=0
1889	DO ir = 1, nn_rimwidth(ib_bdy)
1890	! DO ii = jpisdt(iseg), jpisft(iseg) - 1
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	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1898	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1899	ENDDO
1900	!
1901	! ------------ V points -------------
1902	igrd=3
1903	icount=0
1904	DO ir = 1, nn_rimwidth(ib_bdy)
1905	DO ii = jpisdt(iseg), jpisft(iseg)
1906	icount = icount + 1
1907	nbidta(icount, igrd, ib_bdy) = ii + nn_hls
1908	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 + nn_hls
1909	nbrdta(icount, igrd, ib_bdy) = ir
1910	ENDDO
1911	ENDDO
1912	ENDDO
1913
1914
1915	END SUBROUTINE bdy_coords_seg
1916
1917
1918	SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1919	!!----------------------------------------------------------------------
1920	!! * ROUTINE bdy_ctl_corn *
1921	!!
1922	!! ** Purpose : Check numerical schemes consistency between
1923	!! segments having a common corner
1924	!!
1925	!! ** Method :
1926	!!----------------------------------------------------------------------
1927	INTEGER, INTENT(in) :: ib1, ib2
1928	INTEGER :: itest
1929	!!----------------------------------------------------------------------
1930	itest = 0
1931
1932	IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) ) itest = itest + 1
1933	IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) ) itest = itest + 1
1934	IF( cn_tra (ib1) /= cn_tra (ib2) ) itest = itest + 1
1935	!
1936	IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) ) itest = itest + 1
1937	IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) ) itest = itest + 1
1938	IF( nn_tra_dta (ib1) /= nn_tra_dta (ib2) ) itest = itest + 1
1939	!
1940	IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) ) itest = itest + 1
1941	!
1942	IF( itest>0 ) THEN
1943	WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2
1944	CALL ctl_stop( ctmp1, ' have different open bdy schemes' )
1945	ENDIF
1946	!
1947	END SUBROUTINE bdy_ctl_corn
1948
1949
1950	SUBROUTINE bdy_meshwri()
1951	!!----------------------------------------------------------------------
1952	!! * ROUTINE bdy_meshwri *
1953	!!
1954	!! ** Purpose : write netcdf file with nbr, flagu, flagv, ntreat for T, U
1955	!! and V points in 2D arrays for easier visualisation/control
1956	!!
1957	!! ** Method : use iom_rstput as in domwri.F
1958	!!----------------------------------------------------------------------
1959	INTEGER :: ib_bdy, ii, ij, igrd, ib ! dummy loop indices
1960	INTEGER :: inum ! - -
1961	REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! pointer to 2D mask fields
1962	REAL(wp) , DIMENSION(jpi,jpj) :: ztmp
1963	CHARACTER(LEN=1) , DIMENSION(jpbgrd) :: cgrid
1964	!!----------------------------------------------------------------------
1965	cgrid = (/'t','u','v'/)
1966	CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. )
1967	DO igrd = 1, jpbgrd
1968	SELECT CASE( igrd )
1969	CASE( 1 ) ; zmask => tmask(:,:,1)
1970	CASE( 2 ) ; zmask => umask(:,:,1)
1971	CASE( 3 ) ; zmask => vmask(:,:,1)
1972	END SELECT
1973	ztmp(:,:) = zmask(:,:)
1974	DO ib_bdy = 1, nb_bdy
1975	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) ! nbr deined for all rims
1976	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1977	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1978	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10.
1979	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1980	END DO
1981	END DO
1982	CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1983	ztmp(:,:) = zmask(:,:)
1984	DO ib_bdy = 1, nb_bdy
1985	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! flagu defined only for rims 0 and 1
1986	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1987	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1988	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10.
1989	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1990	END DO
1991	END DO
1992	CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1993	ztmp(:,:) = zmask(:,:)
1994	DO ib_bdy = 1, nb_bdy
1995	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! flagv defined only for rims 0 and 1
1996	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1997	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1998	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10.
1999	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
2000	END DO
2001	END DO
2002	CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
2003	ztmp(:,:) = zmask(:,:)
2004	DO ib_bdy = 1, nb_bdy
2005	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! ntreat defined only for rims 0 and 1
2006	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
2007	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
2008	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10.
2009	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
2010	END DO
2011	END DO
2012	CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
2013	END DO
2014	CALL iom_close( inum )
2015
2016	END SUBROUTINE bdy_meshwri
2017
2018	!!=================================================================================
2019	END MODULE bdyini

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

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

Context Navigation

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

Download in other formats: