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/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/BDY – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/BDY/bdyini.F90 @ 11609

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

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

Download in other formats: