Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

bdyini.F90 @ 12658

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

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.2_mirror/src/OCE/BDY/bdyini.F90 @ 12658

Download in other formats: