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 @ 15360

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

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

Download in other formats: