Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

bdyini.F90 @ 12094

Last change on this file since 12094 was 12094, checked in by deazer, 5 years ago
Allow ts_init to work properly for already interpolated scoords Allow bdy to be spread over arbitary rim points
File size: 89.4 KB

Line
1	MODULE bdyini
2	!!======================================================================
3	!! * MODULE bdyini *
4	!! Unstructured open boundaries : initialisation
5	!!======================================================================
6	!! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code
7	!! - ! 2007-01 (D. Storkey) Update to use IOM module
8	!! - ! 2007-01 (D. Storkey) Tidal forcing
9	!! 3.0 ! 2008-04 (NEMO team) add in the reference version
10	!! 3.3 ! 2010-09 (E.O'Dea) updates for Shelf configurations
11	!! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions
12	!! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge
13	!! 3.4 ! 2012 (J. Chanut) straight open boundary case update
14	!! 3.5 ! 2012 (S. Mocavero, I. Epicoco) optimization of BDY communications
15	!! 3.7 ! 2016 (T. Lovato) Remove bdy macro, call here init for dta and tides
16	!!----------------------------------------------------------------------
17	!! bdy_init : Initialization of unstructured open boundaries
18	!!----------------------------------------------------------------------
19	USE oce ! ocean dynamics and tracers variables
20	USE dom_oce ! ocean space and time domain
21	USE 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
401	! Calculate global boundary index arrays or read in from file
402	!------------------------------------------------------------
403	! 1. Read global index arrays from boundary coordinates file.
404	DO ib_bdy = 1, nb_bdy
405	!
406	IF( ln_coords_file(ib_bdy) ) THEN
407	!
408	ALLOCATE( zz_read( MAXVAL(nblendta), 1 ) )
409	CALL iom_open( cn_coords_file(ib_bdy), inum )
410	!
411	DO igrd = 1, jpbgrd
412	CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
413	DO ii = 1,nblendta(igrd,ib_bdy)
414	nbidta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
415	END DO
416	CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
417	DO ii = 1,nblendta(igrd,ib_bdy)
418	nbjdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
419	END DO
420	CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
421	DO ii = 1,nblendta(igrd,ib_bdy)
422	nbrdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
423	END DO
424	!
425	ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
426	IF(lwp) WRITE(numout,*)
427	IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
428	IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
429	IF (ibr_max < nn_rimwidth(ib_bdy)) &
430	CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
431	END DO
432	!
433	CALL iom_close( inum )
434	DEALLOCATE( zz_read )
435	!
436	ENDIF
437	!
438	END DO
439
440	! 2. Now fill indices corresponding to straight open boundary arrays:
441	CALL bdy_coords_seg( nbidta, nbjdta, nbrdta )
442
443	! Deal with duplicated points
444	!-----------------------------
445	! We assign negative indices to duplicated points (to remove them from bdy points to be updated)
446	! if their distance to the bdy is greater than the other
447	! If their distance are the same, just keep only one to avoid updating a point twice
448	DO igrd = 1, jpbgrd
449	DO ib_bdy1 = 1, nb_bdy
450	DO ib_bdy2 = 1, nb_bdy
451	IF (ib_bdy1/=ib_bdy2) THEN
452	DO ib1 = 1, nblendta(igrd,ib_bdy1)
453	DO ib2 = 1, nblendta(igrd,ib_bdy2)
454	IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
455	& (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
456	! IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &
457	! & nbidta(ib1, igrd, ib_bdy1), &
458	! & nbjdta(ib2, igrd, ib_bdy2)
459	! keep only points with the lowest distance to boundary:
460	IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN
461	nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
462	nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
463	ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN
464	nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
465	nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
466	! Arbitrary choice if distances are the same:
467	ELSE
468	nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
469	nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
470	ENDIF
471	END IF
472	END DO
473	END DO
474	ENDIF
475	END DO
476	END DO
477	END DO
478	!
479	! Find lenght of boundaries and rim on local mpi domain
480	!------------------------------------------------------
481	!
482	iwe = mig(1)
483	ies = mig(jpi)
484	iso = mjg(1)
485	ino = mjg(jpj)
486	!
487	DO ib_bdy = 1, nb_bdy
488	DO igrd = 1, jpbgrd
489	icount = 0 ! initialization of local bdy length
490	icountr = 0 ! initialization of local rim 0 and rim 1 bdy length
491	icountr0 = 0 ! initialization of local rim 0 bdy length
492	idx_bdy(ib_bdy)%nblen(igrd) = 0
493	idx_bdy(ib_bdy)%nblenrim(igrd) = 0
494	idx_bdy(ib_bdy)%nblenrim0(igrd) = 0
495	DO ib = 1, nblendta(igrd,ib_bdy)
496	! check that data is in correct order in file
497	IF( ib > 1 ) THEN
498	IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ib-1,igrd,ib_bdy) ) THEN
499	CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', &
500	& ' in order of distance from edge nbr A utility for re-ordering ', &
501	& ' boundary coordinates and data files exists in the TOOLS/OBC directory')
502	ENDIF
503	ENDIF
504	! check if point is in local domain
505	IF( nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND. &
506	& nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino ) THEN
507	!
508	icount = icount + 1
509	IF( nbrdta(ib,igrd,ib_bdy) == 1 .OR. nbrdta(ib,igrd,ib_bdy) == 0 ) icountr = icountr + 1
510	IF( nbrdta(ib,igrd,ib_bdy) == 0 ) icountr0 = icountr0 + 1
511	ENDIF
512	END DO
513	idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc
514	idx_bdy(ib_bdy)%nblenrim (igrd) = icountr !: length of rim 0 and rim 1 boundary data on each proc
515	idx_bdy(ib_bdy)%nblenrim0(igrd) = icountr0 !: length of rim 0 boundary data on each proc
516	END DO ! igrd
517
518	! Allocate index arrays for this boundary set
519	!--------------------------------------------
520	ilen1 = MAXVAL( idx_bdy(ib_bdy)%nblen(:) )
521	ALLOCATE( idx_bdy(ib_bdy)%nbi (ilen1,jpbgrd) , &
522	& idx_bdy(ib_bdy)%nbj (ilen1,jpbgrd) , &
523	& idx_bdy(ib_bdy)%nbr (ilen1,jpbgrd) , &
524	& idx_bdy(ib_bdy)%nbd (ilen1,jpbgrd) , &
525	& idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) , &
526	& idx_bdy(ib_bdy)%ntreat(ilen1,jpbgrd) , &
527	& idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) , &
528	& idx_bdy(ib_bdy)%nbw (ilen1,jpbgrd) , &
529	& idx_bdy(ib_bdy)%flagu (ilen1,jpbgrd) , &
530	& idx_bdy(ib_bdy)%flagv (ilen1,jpbgrd) )
531
532	! Dispatch mapping indices and discrete distances on each processor
533	! -----------------------------------------------------------------
534	DO igrd = 1, jpbgrd
535	icount = 0
536	! Outer loop on rimwidth to ensure outermost points come first in the local arrays.
537	DO ir = 0, nn_rimwidth(ib_bdy)
538	DO ib = 1, nblendta(igrd,ib_bdy)
539	! check if point is in local domain and equals ir
540	IF( nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND. &
541	& nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino .AND. &
542	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
543	!
544	icount = icount + 1
545	idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1 ! global to local indexes
546	idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 ! global to local indexes
547	idx_bdy(ib_bdy)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_bdy)
548	idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
549	ENDIF
550	END DO
551	END DO
552	END DO ! igrd
553
554	END DO ! ib_bdy
555
556	! Initialize array indicating communications in bdy
557	! -------------------------------------------------
558	ALLOCATE( lsend_bdy(nb_bdy,jpbgrd,4,0:1), lrecv_bdy(nb_bdy,jpbgrd,4,0:1) )
559	lsend_bdy(:,:,:,:) = .false.
560	lrecv_bdy(:,:,:,:) = .false.
561
562	DO ib_bdy = 1, nb_bdy
563	DO igrd = 1, jpbgrd
564	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! only the rim triggers communications, see bdy routines
565	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
566	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
567	IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN ; ir = 0
568	ELSE ; ir = 1
569	END IF
570	!
571	! check if point has to be sent to a neighbour
572	! W neighbour and on the inner left side
573	IF( ii == 2 .and. (nbondi == 0 .or. nbondi == 1) ) lsend_bdy(ib_bdy,igrd,1,ir) = .true.
574	! E neighbour and on the inner right side
575	IF( ii == jpi-1 .and. (nbondi == 0 .or. nbondi == -1) ) lsend_bdy(ib_bdy,igrd,2,ir) = .true.
576	! S neighbour and on the inner down side
577	IF( ij == 2 .and. (nbondj == 0 .or. nbondj == 1) ) lsend_bdy(ib_bdy,igrd,3,ir) = .true.
578	! N neighbour and on the inner up side
579	IF( ij == jpj-1 .and. (nbondj == 0 .or. nbondj == -1) ) lsend_bdy(ib_bdy,igrd,4,ir) = .true.
580	!
581	! check if point has to be received from a neighbour
582	! W neighbour and on the outter left side
583	IF( ii == 1 .and. (nbondi == 0 .or. nbondi == 1) ) lrecv_bdy(ib_bdy,igrd,1,ir) = .true.
584	! E neighbour and on the outter right side
585	IF( ii == jpi .and. (nbondi == 0 .or. nbondi == -1) ) lrecv_bdy(ib_bdy,igrd,2,ir) = .true.
586	! S neighbour and on the outter down side
587	IF( ij == 1 .and. (nbondj == 0 .or. nbondj == 1) ) lrecv_bdy(ib_bdy,igrd,3,ir) = .true.
588	! N neighbour and on the outter up side
589	IF( ij == jpj .and. (nbondj == 0 .or. nbondj == -1) ) lrecv_bdy(ib_bdy,igrd,4,ir) = .true.
590	!
591	END DO
592	END DO ! igrd
593
594	! Compute rim weights for FRS scheme
595	! ----------------------------------
596	DO igrd = 1, jpbgrd
597	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
598	ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) ) ! both rim 0 and rim 1 have the same weights
599	! idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 ) ! tanh formulation
600	idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( ir - 1 ) * 0.5 &
601	& *(10./FLOAT(nn_rimwidth(ib_bdy))) ) ! JGraham:modified for rim=15
602	! idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic
603	! idx_bdy(ib_bdy)%nbw(ib,igrd) = REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)) ! linear
604	END DO
605	END DO
606
607	! Compute damping coefficients
608	! ----------------------------
609	DO igrd = 1, jpbgrd
610	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
611	ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) ) ! both rim 0 and rim 1 have the same damping coefficients
612	idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &
613	& (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))*2. ! quadratic
614	idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) &
615	& (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))*2. ! quadratic
616	END DO
617	END DO
618
619	END DO ! ib_bdy
620
621	! ------------------------------------------------------
622	! Initialise masks and find normal/tangential directions
623	! ------------------------------------------------------
624
625	! ------------------------------------------
626	! handle rim0, do as if rim 1 was free ocean
627	! ------------------------------------------
628
629	ztmask(:,:) = tmask(:,:,1) ; zumask(:,:) = umask(:,:,1) ; zvmask(:,:) = vmask(:,:,1)
630	! For the flagu/flagv calculation below we require a version of fmask without
631	! the land boundary condition (shlat) included:
632	DO ij = 1, jpjm1
633	DO ii = 1, jpim1
634	zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) &
635	& * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
636	END DO
637	END DO
638	CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )
639
640	! Read global 2D mask at T-points: bdytmask
641	! -----------------------------------------
642	! bdytmask = 1 on the computational domain AND on open boundaries
643	! = 0 elsewhere
644
645	bdytmask(:,:) = ssmask(:,:)
646
647	! Derive mask on U and V grid from mask on T grid
648	DO ij = 1, jpjm1
649	DO ii = 1, jpim1
650	bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1,ij )
651	bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii ,ij+1)
652	END DO
653	END DO
654	CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1., bdyvmask, 'V', 1. ) ! Lateral boundary cond.
655
656	! bdy masks are now set to zero on rim 0 points:
657	DO ib_bdy = 1, nb_bdy
658	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1) ! extent of rim 0
659	bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
660	END DO
661	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2) ! extent of rim 0
662	bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
663	END DO
664	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3) ! extent of rim 0
665	bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
666	END DO
667	END DO
668
669	CALL bdy_rim_treat( zumask, zvmask, zfmask, .true. ) ! compute flagu, flagv, ntreat on rim 0
670
671	! ------------------------------------
672	! handle rim1, do as if rim 0 was land
673	! ------------------------------------
674
675	! z[tuv]mask are now set to zero on rim 0 points:
676	DO ib_bdy = 1, nb_bdy
677	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1) ! extent of rim 0
678	ztmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
679	END DO
680	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2) ! extent of rim 0
681	zumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
682	END DO
683	DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3) ! extent of rim 0
684	zvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
685	END DO
686	END DO
687
688	! Recompute zfmask
689	DO ij = 1, jpjm1
690	DO ii = 1, jpim1
691	zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) &
692	& * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
693	END DO
694	END DO
695	CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )
696
697	! bdy masks are now set to zero on rim1 points:
698	DO ib_bdy = 1, nb_bdy
699	DO ib = idx_bdy(ib_bdy)%nblenrim0(1) + 1, idx_bdy(ib_bdy)%nblenrim(1) ! extent of rim 1
700	bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
701	END DO
702	DO ib = idx_bdy(ib_bdy)%nblenrim0(2) + 1, idx_bdy(ib_bdy)%nblenrim(2) ! extent of rim 1
703	bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
704	END DO
705	DO ib = idx_bdy(ib_bdy)%nblenrim0(3) + 1, idx_bdy(ib_bdy)%nblenrim(3) ! extent of rim 1
706	bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
707	END DO
708	END DO
709
710	CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. ) ! compute flagu, flagv, ntreat on rim 1
711	!
712	! Check which boundaries might need communication
713	ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,4,0:1), lrecv_bdyint(nb_bdy,jpbgrd,4,0:1) )
714	lsend_bdyint(:,:,:,:) = .false.
715	lrecv_bdyint(:,:,:,:) = .false.
716	ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,4,0:1), lrecv_bdyext(nb_bdy,jpbgrd,4,0:1) )
717	lsend_bdyext(:,:,:,:) = .false.
718	lrecv_bdyext(:,:,:,:) = .false.
719	!
720	DO igrd = 1, jpbgrd
721	DO ib_bdy = 1, nb_bdy
722	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
723	IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE
724	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
725	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
726	ir = idx_bdy(ib_bdy)%nbr(ib,igrd)
727	flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd))
728	flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd))
729	iibe = ii - flagu ! neighbouring point towards the exterior of the computational domain
730	ijbe = ij - flagv
731	iibi = ii + flagu ! neighbouring point towards the interior of the computational domain
732	ijbi = ij + flagv
733	CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 ) ! free ocean neighbours
734	!
735	! search neighbour in the west/east direction
736	! Rim is on the halo and computed ocean is towards exterior of mpi domain
737	! <-- (o exterior) -->
738	! (1) o\|x OR (2) x\|o
739	! \|___ ___\|
740	IF( iibi == 0 .OR. ii1 == 0 .OR. ii2 == 0 .OR. ii3 == 0 ) lrecv_bdyint(ib_bdy,igrd,1,ir) = .true.
741	IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 ) lrecv_bdyint(ib_bdy,igrd,2,ir) = .true.
742	IF( iibe == 0 ) lrecv_bdyext(ib_bdy,igrd,1,ir) = .true.
743	IF( iibe == jpi+1 ) lrecv_bdyext(ib_bdy,igrd,2,ir) = .true.
744	! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo
745	! :¨¨¨¨¨\|¨¨--> \| \| <--¨¨\|¨¨¨¨¨:
746	! : \| x:o \| neighbour limited by ... would need o \| o:x \| :
747	! :.....\|_._:_____\| (1) W neighbour E neighbour (2) \|_____:_._\|.....:
748	IF( ii == 2 .AND. ( nbondi == 1 .OR. nbondi == 0 ) .AND. &
749	& ( iibi == 3 .OR. ii1 == 3 .OR. ii2 == 3 .OR. ii3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,1,ir)=.true.
750	IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. &
751	& ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) ) lsend_bdyint(ib_bdy,igrd,2,ir)=.true.
752	IF( ii == 2 .AND. ( nbondi == 1 .OR. nbondi == 0 ) .AND. iibe == 3 ) lsend_bdyext(ib_bdy,igrd,1,ir)=.true.
753	IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 ) lsend_bdyext(ib_bdy,igrd,2,ir)=.true.
754	!
755	! search neighbour in the north/south direction
756	! Rim is on the halo and computed ocean is towards exterior of mpi domain
757	!(3) \| \| ^ ___o___
758	! \| \|___x___\| OR \| \| x \|
759	! v o (4) \| \|
760	IF( ijbi == 0 .OR. ij1 == 0 .OR. ij2 == 0 .OR. ij3 == 0 ) lrecv_bdyint(ib_bdy,igrd,3,ir) = .true.
761	IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 ) lrecv_bdyint(ib_bdy,igrd,4,ir) = .true.
762	IF( ijbe == 0 ) lrecv_bdyext(ib_bdy,igrd,3,ir) = .true.
763	IF( ijbe == jpj+1 ) lrecv_bdyext(ib_bdy,igrd,4,ir) = .true.
764	! Check if neighbour has its rim parallel to its mpi subdomain _________ border and next to its halo
765	! ^ \| o \| : :
766	! \| \|¨¨¨¨x¨¨¨¨\| neighbour limited by ... would need o \| \|....x....\|
767	! :_________: (3) S neighbour N neighbour (4) v \| o \|
768	IF( ij == 2 .AND. ( nbondj == 1 .OR. nbondj == 0 ) .AND. &
769	& ( ijbi == 3 .OR. ij1 == 3 .OR. ij2 == 3 .OR. ij3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,3,ir)=.true.
770	IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. &
771	& ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) ) lsend_bdyint(ib_bdy,igrd,4,ir)=.true.
772	IF( ij == 2 .AND. ( nbondj == 1 .OR. nbondj == 0 ) .AND. ijbe == 3 ) lsend_bdyext(ib_bdy,igrd,3,ir)=.true.
773	IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 ) lsend_bdyext(ib_bdy,igrd,4,ir)=.true.
774	END DO
775	END DO
776	END DO
777
778	DO ib_bdy = 1,nb_bdy
779	IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &
780	& cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &
781	& cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo' ) THEN
782	DO igrd = 1, jpbgrd
783	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
784	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
785	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
786	IF( mig(ii) > 2 .AND. mig(ii) < jpiglo-2 .AND. mjg(ij) > 2 .AND. mjg(ij) < jpjglo-2 ) THEN
787	WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain'
788	CALL ctl_stop( ctmp1 )
789	END IF
790	END DO
791	END DO
792	END IF
793	END DO
794	!
795	DEALLOCATE( nbidta, nbjdta, nbrdta )
796	!
797	END SUBROUTINE bdy_def
798
799
800	SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 )
801	!!----------------------------------------------------------------------
802	!! * ROUTINE bdy_rim_treat *
803	!!
804	!! ** Purpose : Initialize structures ( flagu, flagv, ntreat ) indicating how rim points
805	!! are to be handled in the boundary condition treatment
806	!!
807	!! ** Method : - to handle rim 0 zmasks must indicate ocean points (set at one on rim 0 and rim 1 and interior)
808	!! and bdymasks must be set at 0 on rim 0 (set at one on rim 1 and interior)
809	!! (as if rim 1 was free ocean)
810	!! - to handle rim 1 zmasks must be set at 0 on rim 0 (set at one on rim 1 and interior)
811	!! and bdymasks must indicate free ocean points (set at one on interior)
812	!! (as if rim 0 was land)
813	!! - we can then check in which direction the interior of the computational domain is with the difference
814	!! mask array values on both sides to compute flagu and flagv
815	!! - and look at the ocean neighbours to compute ntreat
816	!!----------------------------------------------------------------------
817	REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: pfmask ! temporary fmask excluding coastal boundary condition (shlat)
818	REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: pumask, pvmask ! temporary t/u/v mask array
819	LOGICAL , INTENT (in ) :: lrim0 ! .true. -> rim 0 .false. -> rim 1
820	INTEGER :: ib_bdy, ii, ij, igrd, ib, icount ! dummy loop indices
821	INTEGER :: i_offset, j_offset, inn ! local integer
822	INTEGER :: ibeg, iend ! local integer
823	LOGICAL :: llnon, llson, llean, llwen ! local logicals indicating the presence of a ocean neighbour
824	REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! pointer to 2D mask fields
825	REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars
826	CHARACTER(LEN=1), DIMENSION(jpbgrd) :: cgrid
827	REAL(wp) , DIMENSION(jpi,jpj) :: ztmp
828	!!----------------------------------------------------------------------
829
830	cgrid = (/'t','u','v'/)
831
832	DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components
833
834	! Calculate relationship of U direction to the local orientation of the boundary
835	! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
836	! flagu = 0 : u is tangential
837	! flagu = 1 : u is normal to the boundary and is direction is inward
838	DO igrd = 1, jpbgrd
839	SELECT CASE( igrd )
840	CASE( 1 ) ; zmask => pumask ; i_offset = 0
841	CASE( 2 ) ; zmask => bdytmask ; i_offset = 1
842	CASE( 3 ) ; zmask => pfmask ; i_offset = 0
843	END SELECT
844	icount = 0
845	ztmp(:,:) = -999._wp
846	IF( lrim0 ) THEN ! extent of rim 0
847	ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
848	ELSE ! extent of rim 1
849	ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1 ; iend = idx_bdy(ib_bdy)%nblenrim(igrd)
850	END IF
851	DO ib = ibeg, iend
852	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
853	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
854	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE
855	zwfl = zmask(ii+i_offset-1,ij)
856	zefl = zmask(ii+i_offset ,ij)
857	! This error check only works if you are using the bdyXmask arrays
858	IF( i_offset == 1 .and. zefl + zwfl == 2. ) THEN
859	icount = icount + 1
860	IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
861	ELSE
862	ztmp(ii,ij) = -zwfl + zefl
863	ENDIF
864	END DO
865	IF( icount /= 0 ) THEN
866	WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,', &
867	' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
868	CALL ctl_stop( ctmp1 )
869	ENDIF
870	SELECT CASE( igrd )
871	CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )
872	CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. )
873	CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. )
874	END SELECT
875	DO ib = ibeg, iend
876	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
877	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
878	idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij)
879	END DO
880	END DO
881
882	! Calculate relationship of V direction to the local orientation of the boundary
883	! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
884	! flagv = 0 : v is tangential
885	! flagv = 1 : v is normal to the boundary and is direction is inward
886	DO igrd = 1, jpbgrd
887	SELECT CASE( igrd )
888	CASE( 1 ) ; zmask => pvmask ; j_offset = 0
889	CASE( 2 ) ; zmask => pfmask ; j_offset = 0
890	CASE( 3 ) ; zmask => bdytmask ; j_offset = 1
891	END SELECT
892	icount = 0
893	ztmp(:,:) = -999._wp
894	IF( lrim0 ) THEN ! extent of rim 0
895	ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
896	ELSE ! extent of rim 1
897	ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1 ; iend = idx_bdy(ib_bdy)%nblenrim(igrd)
898	END IF
899	DO ib = ibeg, iend
900	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
901	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
902	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE
903	zsfl = zmask(ii,ij+j_offset-1)
904	znfl = zmask(ii,ij+j_offset )
905	! This error check only works if you are using the bdyXmask arrays
906	IF( j_offset == 1 .and. znfl + zsfl == 2. ) THEN
907	IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
908	icount = icount + 1
909	ELSE
910	ztmp(ii,ij) = -zsfl + znfl
911	END IF
912	END DO
913	IF( icount /= 0 ) THEN
914	WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,', &
915	' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
916	CALL ctl_stop( ctmp1 )
917	ENDIF
918	SELECT CASE( igrd )
919	CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )
920	CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. )
921	CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. )
922	END SELECT
923	DO ib = ibeg, iend
924	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
925	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
926	idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij)
927	END DO
928	END DO
929	!
930	END DO ! ib_bdy
931
932	DO ib_bdy = 1, nb_bdy
933	DO igrd = 1, jpbgrd
934	SELECT CASE( igrd )
935	CASE( 1 ) ; zmask => bdytmask
936	CASE( 2 ) ; zmask => bdyumask
937	CASE( 3 ) ; zmask => bdyvmask
938	END SELECT
939	ztmp(:,:) = -999._wp
940	IF( lrim0 ) THEN ! extent of rim 0
941	ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
942	ELSE ! extent of rim 1
943	ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1 ; iend = idx_bdy(ib_bdy)%nblenrim(igrd)
944	END IF
945	DO ib = ibeg, iend
946	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
947	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
948	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE
949	llnon = zmask(ii ,ij+1) == 1.
950	llson = zmask(ii ,ij-1) == 1.
951	llean = zmask(ii+1,ij ) == 1.
952	llwen = zmask(ii-1,ij ) == 1.
953	inn = COUNT( (/ llnon, llson, llean, llwen /) )
954	IF( inn == 0 ) THEN ! no neighbours -> interior of a corner or cluster of rim points
955	! ! ! _____ ! _____ ! __ __
956	! 1 \| o ! 2 o \| ! 3 \| x ! 4 x \| ! \| \| -> error
957	! \|_x_ _ ! _ _x_\| ! \| o ! o \| ! \|x_x\|
958	IF( zmask(ii+1,ij+1) == 1. ) THEN ; ztmp(ii,ij) = 1.
959	ELSEIF( zmask(ii-1,ij+1) == 1. ) THEN ; ztmp(ii,ij) = 2.
960	ELSEIF( zmask(ii+1,ij-1) == 1. ) THEN ; ztmp(ii,ij) = 3.
961	ELSEIF( zmask(ii-1,ij-1) == 1. ) THEN ; ztmp(ii,ij) = 4.
962	ELSE ; ztmp(ii,ij) = -1.
963	WRITE(ctmp1,*) 'Problem with ',cgrid(igrd) ,' grid point', ii, ij, &
964	' on boundary set ', ib_bdy, ' has no free ocean neighbour'
965	IF( lrim0 ) THEN
966	WRITE(ctmp2,*) ' There seems to be a cluster of rim 0 points.'
967	ELSE
968	WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.'
969	END IF
970	CALL ctl_warn( ctmp1, ctmp2 )
971	END IF
972	END IF
973	IF( inn == 1 ) THEN ! middle of linear bdy or incomplete corner ! ___ o
974	! \| ! \| ! o ! ______ ! \|x___
975	! 5 \| x o ! 6 o x \| ! 7 __x__ ! 8 x
976	! \| ! \| ! ! o
977	IF( llean ) ztmp(ii,ij) = 5.
978	IF( llwen ) ztmp(ii,ij) = 6.
979	IF( llnon ) ztmp(ii,ij) = 7.
980	IF( llson ) ztmp(ii,ij) = 8.
981	END IF
982	IF( inn == 2 ) THEN ! exterior of a corner
983	! o ! o ! _____\| ! \|_____
984	! 9 ____x o ! 10 o x___ ! 11 x o ! 12 o x
985	! \| ! \| ! o ! o
986	IF( llnon .AND. llean ) ztmp(ii,ij) = 9.
987	IF( llnon .AND. llwen ) ztmp(ii,ij) = 10.
988	IF( llson .AND. llean ) ztmp(ii,ij) = 11.
989	IF( llson .AND. llwen ) ztmp(ii,ij) = 12.
990	END IF
991	IF( inn == 3 ) THEN ! 3 neighbours __ __
992	! \|_ o ! o _\| ! \|_\| ! o
993	! 13 _\| x o ! 14 o x \|_ ! 15 o x o ! 16 o x o
994	! \| o ! o \| ! o ! __\|¨\|__
995	IF( llnon .AND. llean .AND. llson ) ztmp(ii,ij) = 13.
996	IF( llnon .AND. llwen .AND. llson ) ztmp(ii,ij) = 14.
997	IF( llwen .AND. llson .AND. llean ) ztmp(ii,ij) = 15.
998	IF( llwen .AND. llnon .AND. llean ) ztmp(ii,ij) = 16.
999	END IF
1000	IF( inn == 4 ) THEN
1001	WRITE(ctmp1,*) 'Problem with ',cgrid(igrd) ,' grid point', ii, ij, &
1002	' on boundary set ', ib_bdy, ' have 4 neighbours'
1003	CALL ctl_stop( ctmp1 )
1004	END IF
1005	END DO
1006	SELECT CASE( igrd )
1007	CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )
1008	CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. )
1009	CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. )
1010	END SELECT
1011	DO ib = ibeg, iend
1012	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1013	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1014	idx_bdy(ib_bdy)%ntreat(ib,igrd) = NINT(ztmp(ii,ij))
1015	END DO
1016	END DO
1017	END DO
1018
1019	END SUBROUTINE bdy_rim_treat
1020
1021
1022	SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )
1023	!!----------------------------------------------------------------------
1024	!! * ROUTINE find_neib *
1025	!!
1026	!! ** Purpose : get ii1, ij1, ii2, ij2, ii3, ij3, the indices of
1027	!! the free ocean neighbours of (ii,ij) for bdy treatment
1028	!!
1029	!! ** Method : use itreat input to select a case
1030	!! N.B. ntreat is defined for all bdy points in routine bdy_rim_treat
1031	!!
1032	!!----------------------------------------------------------------------
1033	INTEGER, INTENT(in ) :: ii, ij, itreat
1034	INTEGER, INTENT( out) :: ii1, ij1, ii2, ij2, ii3, ij3
1035	!!----------------------------------------------------------------------
1036	SELECT CASE( itreat ) ! points that will be used by bdy routines, -1 will be discarded
1037	! ! ! _____ ! _____
1038	! 1 \| o ! 2 o \| ! 3 \| x ! 4 x \|
1039	! \|_x_ _ ! _ _x_\| ! \| o ! o \|
1040	CASE( 1 ) ; ii1 = ii+1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1041	CASE( 2 ) ; ii1 = ii-1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1042	CASE( 3 ) ; ii1 = ii+1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1043	CASE( 4 ) ; ii1 = ii-1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1044	! \| ! \| ! o ! ______ ! or incomplete corner
1045	! 5 \| x o ! 6 o x \| ! 7 __x__ ! 8 x ! 7 ____ o
1046	! \| ! \| ! ! o ! \|x___
1047	CASE( 5 ) ; ii1 = ii+1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1048	CASE( 6 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1049	CASE( 7 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1050	CASE( 8 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1
1051	! o ! o ! _____\| ! \|_____
1052	! 9 ____x o ! 10 o x___ ! 11 x o ! 12 o x
1053	! \| ! \| ! o ! o
1054	CASE( 9 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1055	CASE( 10 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1056	CASE( 11 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1057	CASE( 12 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1
1058	! \|_ o ! o _\| ! ¨¨\|_\|¨¨ ! o
1059	! 13 _\| x o ! 14 o x \|_ ! 15 o x o ! 16 o x o
1060	! \| o ! o \| ! o ! __\|¨\|__
1061	CASE( 13 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1
1062	CASE( 14 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1
1063	CASE( 15 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij-1 ; ii3 = ii+1 ; ij3 = ij
1064	CASE( 16 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij+1 ; ii3 = ii+1 ; ij3 = ij
1065	END SELECT
1066	END SUBROUTINE find_neib
1067
1068
1069	SUBROUTINE bdy_read_seg( kb_bdy, knblendta )
1070	!!----------------------------------------------------------------------
1071	!! * ROUTINE bdy_coords_seg *
1072	!!
1073	!! ** Purpose : build bdy coordinates with segments defined in namelist
1074	!!
1075	!! ** Method : read namelist nambdy_index blocks
1076	!!
1077	!!----------------------------------------------------------------------
1078	INTEGER , INTENT (in ) :: kb_bdy ! bdy number
1079	INTEGER, DIMENSION(jpbgrd), INTENT ( out) :: knblendta ! length of index arrays
1080	!!
1081	INTEGER :: ios ! Local integer output status for namelist read
1082	INTEGER :: nbdyind, nbdybeg, nbdyend
1083	CHARACTER(LEN=1) :: ctypebdy ! - -
1084	NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
1085	!!----------------------------------------------------------------------
1086
1087	! No REWIND here because may need to read more than one nambdy_index namelist.
1088	! Read only namelist_cfg to avoid unseccessfull overwrite
1089	! keep full control of the configuration namelist
1090	READ ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 )
1091	904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' )
1092	IF(lwm) WRITE ( numond, nambdy_index )
1093
1094	SELECT CASE ( TRIM(ctypebdy) )
1095	CASE( 'N' )
1096	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1097	nbdyind = jpjglo - 2 ! set boundary to whole side of model domain.
1098	nbdybeg = 2
1099	nbdyend = jpiglo - 1
1100	ENDIF
1101	nbdysegn = nbdysegn + 1
1102	npckgn(nbdysegn) = kb_bdy ! Save bdy package number
1103	jpjnob(nbdysegn) = nbdyind
1104	jpindt(nbdysegn) = nbdybeg
1105	jpinft(nbdysegn) = nbdyend
1106	!
1107	CASE( 'S' )
1108	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1109	nbdyind = 2 ! set boundary to whole side of model domain.
1110	nbdybeg = 2
1111	nbdyend = jpiglo - 1
1112	ENDIF
1113	nbdysegs = nbdysegs + 1
1114	npckgs(nbdysegs) = kb_bdy ! Save bdy package number
1115	jpjsob(nbdysegs) = nbdyind
1116	jpisdt(nbdysegs) = nbdybeg
1117	jpisft(nbdysegs) = nbdyend
1118	!
1119	CASE( 'E' )
1120	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1121	nbdyind = jpiglo - 2 ! set boundary to whole side of model domain.
1122	nbdybeg = 2
1123	nbdyend = jpjglo - 1
1124	ENDIF
1125	nbdysege = nbdysege + 1
1126	npckge(nbdysege) = kb_bdy ! Save bdy package number
1127	jpieob(nbdysege) = nbdyind
1128	jpjedt(nbdysege) = nbdybeg
1129	jpjeft(nbdysege) = nbdyend
1130	!
1131	CASE( 'W' )
1132	IF( nbdyind == -1 ) THEN ! Automatic boundary definition: if nbdysegX = -1
1133	nbdyind = 2 ! set boundary to whole side of model domain.
1134	nbdybeg = 2
1135	nbdyend = jpjglo - 1
1136	ENDIF
1137	nbdysegw = nbdysegw + 1
1138	npckgw(nbdysegw) = kb_bdy ! Save bdy package number
1139	jpiwob(nbdysegw) = nbdyind
1140	jpjwdt(nbdysegw) = nbdybeg
1141	jpjwft(nbdysegw) = nbdyend
1142	!
1143	CASE DEFAULT ; CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
1144	END SELECT
1145
1146	! For simplicity we assume that in case of straight bdy, arrays have the same length
1147	! (even if it is true that last tangential velocity points
1148	! are useless). This simplifies a little bit boundary data format (and agrees with format
1149	! used so far in obc package)
1150
1151	knblendta(1:jpbgrd) = (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy)
1152
1153	END SUBROUTINE bdy_read_seg
1154
1155
1156	SUBROUTINE bdy_ctl_seg
1157	!!----------------------------------------------------------------------
1158	!! * ROUTINE bdy_ctl_seg *
1159	!!
1160	!! ** Purpose : Check straight open boundary segments location
1161	!!
1162	!! ** Method : - Look for open boundary corners
1163	!! - Check that segments start or end on land
1164	!!----------------------------------------------------------------------
1165	INTEGER :: ib, ib1, ib2, ji ,jj, itest
1166	INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns
1167	REAL(wp), DIMENSION(2) :: ztestmask
1168	!!----------------------------------------------------------------------
1169	!
1170	IF (lwp) WRITE(numout,*) ' '
1171	IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1172	IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1173	!
1174	IF(lwp) WRITE(numout,*) 'Number of east segments : ', nbdysege
1175	IF(lwp) WRITE(numout,*) 'Number of west segments : ', nbdysegw
1176	IF(lwp) WRITE(numout,*) 'Number of north segments : ', nbdysegn
1177	IF(lwp) WRITE(numout,*) 'Number of south segments : ', nbdysegs
1178	! 1. Check bounds
1179	!----------------
1180	DO ib = 1, nbdysegn
1181	IF (lwp) WRITE(numout,) '*check north seg bounds pckg: ', npckgn(ib)
1182	IF ((jpjnob(ib).ge.jpjglo-1).or.&
1183	&(jpjnob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1184	IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1185	IF (jpindt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1186	IF (jpinft(ib).gt.jpiglo) CALL ctl_stop( 'End index out of domain' )
1187	END DO
1188	!
1189	DO ib = 1, nbdysegs
1190	IF (lwp) WRITE(numout,) '*check south seg bounds pckg: ', npckgs(ib)
1191	IF ((jpjsob(ib).ge.jpjglo-1).or.&
1192	&(jpjsob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1193	IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1194	IF (jpisdt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1195	IF (jpisft(ib).gt.jpiglo) CALL ctl_stop( 'End index out of domain' )
1196	END DO
1197	!
1198	DO ib = 1, nbdysege
1199	IF (lwp) WRITE(numout,) '*check east seg bounds pckg: ', npckge(ib)
1200	IF ((jpieob(ib).ge.jpiglo-1).or.&
1201	&(jpieob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1202	IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1203	IF (jpjedt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1204	IF (jpjeft(ib).gt.jpjglo) CALL ctl_stop( 'End index out of domain' )
1205	END DO
1206	!
1207	DO ib = 1, nbdysegw
1208	IF (lwp) WRITE(numout,) '*check west seg bounds pckg: ', npckgw(ib)
1209	IF ((jpiwob(ib).ge.jpiglo-1).or.&
1210	&(jpiwob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
1211	IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1212	IF (jpjwdt(ib).lt.1 ) CALL ctl_stop( 'Start index out of domain' )
1213	IF (jpjwft(ib).gt.jpjglo) CALL ctl_stop( 'End index out of domain' )
1214	ENDDO
1215	!
1216	!
1217	! 2. Look for segment crossings
1218	!------------------------------
1219	IF (lwp) WRITE(numout,) '*Look for segments corners :'
1220	!
1221	itest = 0 ! corner number
1222	!
1223	! flag to detect if start or end of open boundary belongs to a corner
1224	! if not (=0), it must be on land.
1225	! if a corner is detected, save bdy package number for further tests
1226	icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1227	! South/West crossings
1228	IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1229	DO ib1 = 1, nbdysegw
1230	DO ib2 = 1, nbdysegs
1231	IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1232	& ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1233	& ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1234	& ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1235	IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN
1236	! We have a possible South-West corner
1237	! WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1238	! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1239	icornw(ib1,1) = npckgs(ib2)
1240	icorns(ib2,1) = npckgw(ib1)
1241	ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1242	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1243	& jpisft(ib2), jpjwft(ib1)
1244	WRITE(ctmp2,*) ' Not allowed yet'
1245	WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1246	& ' and South segment: ',npckgs(ib2)
1247	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1248	ELSE
1249	WRITE(ctmp1,*) ' Check South and West Open boundary indices'
1250	WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , &
1251	& ' and South segment: ',npckgs(ib2)
1252	CALL ctl_stop( ctmp1, ctmp2 )
1253	END IF
1254	END IF
1255	END DO
1256	END DO
1257	END IF
1258	!
1259	! South/East crossings
1260	IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1261	DO ib1 = 1, nbdysege
1262	DO ib2 = 1, nbdysegs
1263	IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1264	& ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1265	& ( jpjedt(ib1)<=jpjsob(ib2) ).AND. &
1266	& ( jpjeft(ib1)>=jpjsob(ib2) )) THEN
1267	IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1268	! We have a possible South-East corner
1269	! WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1270	! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1271	icorne(ib1,1) = npckgs(ib2)
1272	icorns(ib2,2) = npckge(ib1)
1273	ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1274	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1275	& jpisdt(ib2), jpjeft(ib1)
1276	WRITE(ctmp2,*) ' Not allowed yet'
1277	WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1278	& ' and South segment: ',npckgs(ib2)
1279	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1280	ELSE
1281	WRITE(ctmp1,*) ' Check South and East Open boundary indices'
1282	WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1283	& ' and South segment: ',npckgs(ib2)
1284	CALL ctl_stop( ctmp1, ctmp2 )
1285	END IF
1286	END IF
1287	END DO
1288	END DO
1289	END IF
1290	!
1291	! North/West crossings
1292	IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1293	DO ib1 = 1, nbdysegw
1294	DO ib2 = 1, nbdysegn
1295	IF (( jpindt(ib2)<=jpiwob(ib1) ).AND. &
1296	& ( jpinft(ib2)>=jpiwob(ib1) ).AND. &
1297	& ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1298	& ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1299	IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1300	! We have a possible North-West corner
1301	! WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1302	! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1303	icornw(ib1,2) = npckgn(ib2)
1304	icornn(ib2,1) = npckgw(ib1)
1305	ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1306	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1307	& jpinft(ib2), jpjwdt(ib1)
1308	WRITE(ctmp2,*) ' Not allowed yet'
1309	WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1310	& ' and North segment: ',npckgn(ib2)
1311	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1312	ELSE
1313	WRITE(ctmp1,*) ' Check North and West Open boundary indices'
1314	WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1315	& ' and North segment: ',npckgn(ib2)
1316	CALL ctl_stop( ctmp1, ctmp2 )
1317	END IF
1318	END IF
1319	END DO
1320	END DO
1321	END IF
1322	!
1323	! North/East crossings
1324	IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1325	DO ib1 = 1, nbdysege
1326	DO ib2 = 1, nbdysegn
1327	IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1328	& ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1329	& ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1330	& ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1331	IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1332	! We have a possible North-East corner
1333	! WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1334	! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1335	icorne(ib1,2) = npckgn(ib2)
1336	icornn(ib2,2) = npckge(ib1)
1337	ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1338	WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1339	& jpindt(ib2), jpjedt(ib1)
1340	WRITE(ctmp2,*) ' Not allowed yet'
1341	WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1342	& ' and North segment: ',npckgn(ib2)
1343	CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1344	ELSE
1345	WRITE(ctmp1,*) ' Check North and East Open boundary indices'
1346	WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1347	& ' and North segment: ',npckgn(ib2)
1348	CALL ctl_stop( ctmp1, ctmp2 )
1349	END IF
1350	END IF
1351	END DO
1352	END DO
1353	END IF
1354	!
1355	! 3. Check if segment extremities are on land
1356	!--------------------------------------------
1357	!
1358	! West segments
1359	DO ib = 1, nbdysegw
1360	! get mask at boundary extremities:
1361	ztestmask(1:2)=0.
1362	DO ji = 1, jpi
1363	DO jj = 1, jpj
1364	IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &
1365	& ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1366	IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &
1367	& ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1)
1368	END DO
1369	END DO
1370	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1371
1372	IF (ztestmask(1)==1) THEN
1373	IF (icornw(ib,1)==0) THEN
1374	WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1375	CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1376	ELSE
1377	! This is a corner
1378	IF(lwp) WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1379	CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1380	itest=itest+1
1381	ENDIF
1382	ENDIF
1383	IF (ztestmask(2)==1) THEN
1384	IF (icornw(ib,2)==0) THEN
1385	WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1386	CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' )
1387	ELSE
1388	! This is a corner
1389	IF(lwp) WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1390	CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1391	itest=itest+1
1392	ENDIF
1393	ENDIF
1394	END DO
1395	!
1396	! East segments
1397	DO ib = 1, nbdysege
1398	! get mask at boundary extremities:
1399	ztestmask(1:2)=0.
1400	DO ji = 1, jpi
1401	DO jj = 1, jpj
1402	IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &
1403	& ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
1404	IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &
1405	& ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1)
1406	END DO
1407	END DO
1408	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1409
1410	IF (ztestmask(1)==1) THEN
1411	IF (icorne(ib,1)==0) THEN
1412	WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1413	CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1414	ELSE
1415	! This is a corner
1416	IF(lwp) WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1417	CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1418	itest=itest+1
1419	ENDIF
1420	ENDIF
1421	IF (ztestmask(2)==1) THEN
1422	IF (icorne(ib,2)==0) THEN
1423	WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1424	CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1425	ELSE
1426	! This is a corner
1427	IF(lwp) WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1428	CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1429	itest=itest+1
1430	ENDIF
1431	ENDIF
1432	END DO
1433	!
1434	! South segments
1435	DO ib = 1, nbdysegs
1436	! get mask at boundary extremities:
1437	ztestmask(1:2)=0.
1438	DO ji = 1, jpi
1439	DO jj = 1, jpj
1440	IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &
1441	& ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1442	IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &
1443	& ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1)
1444	END DO
1445	END DO
1446	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1447
1448	IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1449	WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1450	CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1451	ENDIF
1452	IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1453	WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1454	CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1455	ENDIF
1456	END DO
1457	!
1458	! North segments
1459	DO ib = 1, nbdysegn
1460	! get mask at boundary extremities:
1461	ztestmask(1:2)=0.
1462	DO ji = 1, jpi
1463	DO jj = 1, jpj
1464	IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &
1465	& ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
1466	IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &
1467	& ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1)
1468	END DO
1469	END DO
1470	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
1471
1472	IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1473	WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1474	CALL ctl_stop( ctmp1, ' does not start on land' )
1475	ENDIF
1476	IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1477	WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1478	CALL ctl_stop( ctmp1, ' does not end on land' )
1479	ENDIF
1480	END DO
1481	!
1482	IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1483	!
1484	! Other tests TBD:
1485	! segments completly on land
1486	! optimized open boundary array length according to landmask
1487	! Nudging layers that overlap with interior domain
1488	!
1489	END SUBROUTINE bdy_ctl_seg
1490
1491
1492	SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta )
1493	!!----------------------------------------------------------------------
1494	!! * ROUTINE bdy_coords_seg *
1495	!!
1496	!! ** Purpose : build nbidta, nbidta, nbrdta for bdy built with segments
1497	!!
1498	!! ** Method :
1499	!!
1500	!!----------------------------------------------------------------------
1501	INTEGER, DIMENSION(:,:,:), intent( out) :: nbidta, nbjdta, nbrdta ! Index arrays: i and j indices of bdy dta
1502	!!
1503	INTEGER :: ii, ij, ir, iseg
1504	INTEGER :: igrd ! grid type (t=1, u=2, v=3)
1505	INTEGER :: icount !
1506	INTEGER :: ib_bdy ! bdy number
1507	!!----------------------------------------------------------------------
1508
1509	! East
1510	!-----
1511	DO iseg = 1, nbdysege
1512	ib_bdy = npckge(iseg)
1513	!
1514	! ------------ T points -------------
1515	igrd=1
1516	icount=0
1517	DO ir = 1, nn_rimwidth(ib_bdy)
1518	DO ij = jpjedt(iseg), jpjeft(iseg)
1519	icount = icount + 1
1520	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
1521	nbjdta(icount, igrd, ib_bdy) = ij
1522	nbrdta(icount, igrd, ib_bdy) = ir
1523	ENDDO
1524	ENDDO
1525	!
1526	! ------------ U points -------------
1527	igrd=2
1528	icount=0
1529	DO ir = 1, nn_rimwidth(ib_bdy)
1530	DO ij = jpjedt(iseg), jpjeft(iseg)
1531	icount = icount + 1
1532	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir
1533	nbjdta(icount, igrd, ib_bdy) = ij
1534	nbrdta(icount, igrd, ib_bdy) = ir
1535	ENDDO
1536	ENDDO
1537	!
1538	! ------------ V points -------------
1539	igrd=3
1540	icount=0
1541	DO ir = 1, nn_rimwidth(ib_bdy)
1542	! DO ij = jpjedt(iseg), jpjeft(iseg) - 1
1543	DO ij = jpjedt(iseg), jpjeft(iseg)
1544	icount = icount + 1
1545	nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
1546	nbjdta(icount, igrd, ib_bdy) = ij
1547	nbrdta(icount, igrd, ib_bdy) = ir
1548	ENDDO
1549	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1550	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1551	ENDDO
1552	ENDDO
1553	!
1554	! West
1555	!-----
1556	DO iseg = 1, nbdysegw
1557	ib_bdy = npckgw(iseg)
1558	!
1559	! ------------ T points -------------
1560	igrd=1
1561	icount=0
1562	DO ir = 1, nn_rimwidth(ib_bdy)
1563	DO ij = jpjwdt(iseg), jpjwft(iseg)
1564	icount = icount + 1
1565	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
1566	nbjdta(icount, igrd, ib_bdy) = ij
1567	nbrdta(icount, igrd, ib_bdy) = ir
1568	ENDDO
1569	ENDDO
1570	!
1571	! ------------ U points -------------
1572	igrd=2
1573	icount=0
1574	DO ir = 1, nn_rimwidth(ib_bdy)
1575	DO ij = jpjwdt(iseg), jpjwft(iseg)
1576	icount = icount + 1
1577	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
1578	nbjdta(icount, igrd, ib_bdy) = ij
1579	nbrdta(icount, igrd, ib_bdy) = ir
1580	ENDDO
1581	ENDDO
1582	!
1583	! ------------ V points -------------
1584	igrd=3
1585	icount=0
1586	DO ir = 1, nn_rimwidth(ib_bdy)
1587	! DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
1588	DO ij = jpjwdt(iseg), jpjwft(iseg)
1589	icount = icount + 1
1590	nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
1591	nbjdta(icount, igrd, ib_bdy) = ij
1592	nbrdta(icount, igrd, ib_bdy) = ir
1593	ENDDO
1594	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1595	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1596	ENDDO
1597	ENDDO
1598	!
1599	! North
1600	!-----
1601	DO iseg = 1, nbdysegn
1602	ib_bdy = npckgn(iseg)
1603	!
1604	! ------------ T points -------------
1605	igrd=1
1606	icount=0
1607	DO ir = 1, nn_rimwidth(ib_bdy)
1608	DO ii = jpindt(iseg), jpinft(iseg)
1609	icount = icount + 1
1610	nbidta(icount, igrd, ib_bdy) = ii
1611	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
1612	nbrdta(icount, igrd, ib_bdy) = ir
1613	ENDDO
1614	ENDDO
1615	!
1616	! ------------ U points -------------
1617	igrd=2
1618	icount=0
1619	DO ir = 1, nn_rimwidth(ib_bdy)
1620	! DO ii = jpindt(iseg), jpinft(iseg) - 1
1621	DO ii = jpindt(iseg), jpinft(iseg)
1622	icount = icount + 1
1623	nbidta(icount, igrd, ib_bdy) = ii
1624	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
1625	nbrdta(icount, igrd, ib_bdy) = ir
1626	ENDDO
1627	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1628	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1629	ENDDO
1630	!
1631	! ------------ V points -------------
1632	igrd=3
1633	icount=0
1634	DO ir = 1, nn_rimwidth(ib_bdy)
1635	DO ii = jpindt(iseg), jpinft(iseg)
1636	icount = icount + 1
1637	nbidta(icount, igrd, ib_bdy) = ii
1638	nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir
1639	nbrdta(icount, igrd, ib_bdy) = ir
1640	ENDDO
1641	ENDDO
1642	ENDDO
1643	!
1644	! South
1645	!-----
1646	DO iseg = 1, nbdysegs
1647	ib_bdy = npckgs(iseg)
1648	!
1649	! ------------ T points -------------
1650	igrd=1
1651	icount=0
1652	DO ir = 1, nn_rimwidth(ib_bdy)
1653	DO ii = jpisdt(iseg), jpisft(iseg)
1654	icount = icount + 1
1655	nbidta(icount, igrd, ib_bdy) = ii
1656	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
1657	nbrdta(icount, igrd, ib_bdy) = ir
1658	ENDDO
1659	ENDDO
1660	!
1661	! ------------ U points -------------
1662	igrd=2
1663	icount=0
1664	DO ir = 1, nn_rimwidth(ib_bdy)
1665	! DO ii = jpisdt(iseg), jpisft(iseg) - 1
1666	DO ii = jpisdt(iseg), jpisft(iseg)
1667	icount = icount + 1
1668	nbidta(icount, igrd, ib_bdy) = ii
1669	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
1670	nbrdta(icount, igrd, ib_bdy) = ir
1671	ENDDO
1672	nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1673	nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1674	ENDDO
1675	!
1676	! ------------ V points -------------
1677	igrd=3
1678	icount=0
1679	DO ir = 1, nn_rimwidth(ib_bdy)
1680	DO ii = jpisdt(iseg), jpisft(iseg)
1681	icount = icount + 1
1682	nbidta(icount, igrd, ib_bdy) = ii
1683	nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
1684	nbrdta(icount, igrd, ib_bdy) = ir
1685	ENDDO
1686	ENDDO
1687	ENDDO
1688
1689
1690	END SUBROUTINE bdy_coords_seg
1691
1692
1693	SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1694	!!----------------------------------------------------------------------
1695	!! * ROUTINE bdy_ctl_corn *
1696	!!
1697	!! ** Purpose : Check numerical schemes consistency between
1698	!! segments having a common corner
1699	!!
1700	!! ** Method :
1701	!!----------------------------------------------------------------------
1702	INTEGER, INTENT(in) :: ib1, ib2
1703	INTEGER :: itest
1704	!!----------------------------------------------------------------------
1705	itest = 0
1706
1707	IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) ) itest = itest + 1
1708	IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) ) itest = itest + 1
1709	IF( cn_tra (ib1) /= cn_tra (ib2) ) itest = itest + 1
1710	!
1711	IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) ) itest = itest + 1
1712	IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) ) itest = itest + 1
1713	IF( nn_tra_dta (ib1) /= nn_tra_dta (ib2) ) itest = itest + 1
1714	!
1715	IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) ) itest = itest + 1
1716	!
1717	IF( itest>0 ) THEN
1718	WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2
1719	CALL ctl_stop( ctmp1, ' have different open bdy schemes' )
1720	ENDIF
1721	!
1722	END SUBROUTINE bdy_ctl_corn
1723
1724
1725	SUBROUTINE bdy_meshwri()
1726	!!----------------------------------------------------------------------
1727	!! * ROUTINE bdy_meshwri *
1728	!!
1729	!! ** Purpose : write netcdf file with nbr, flagu, flagv, ntreat for T, U
1730	!! and V points in 2D arrays for easier visualisation/control
1731	!!
1732	!! ** Method : use iom_rstput as in domwri.F
1733	!!----------------------------------------------------------------------
1734	INTEGER :: ib_bdy, ii, ij, igrd, ib ! dummy loop indices
1735	INTEGER :: inum ! - -
1736	REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! pointer to 2D mask fields
1737	REAL(wp) , DIMENSION(jpi,jpj) :: ztmp
1738	CHARACTER(LEN=1) , DIMENSION(jpbgrd) :: cgrid
1739	!!----------------------------------------------------------------------
1740	cgrid = (/'t','u','v'/)
1741	CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. )
1742	DO igrd = 1, jpbgrd
1743	SELECT CASE( igrd )
1744	CASE( 1 ) ; zmask => tmask(:,:,1)
1745	CASE( 2 ) ; zmask => umask(:,:,1)
1746	CASE( 3 ) ; zmask => vmask(:,:,1)
1747	END SELECT
1748	ztmp(:,:) = zmask(:,:)
1749	DO ib_bdy = 1, nb_bdy
1750	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) ! nbr deined for all rims
1751	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1752	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1753	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10.
1754	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1755	END DO
1756	END DO
1757	CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1758	ztmp(:,:) = zmask(:,:)
1759	DO ib_bdy = 1, nb_bdy
1760	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! flagu defined only for rims 0 and 1
1761	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1762	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1763	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10.
1764	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1765	END DO
1766	END DO
1767	CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1768	ztmp(:,:) = zmask(:,:)
1769	DO ib_bdy = 1, nb_bdy
1770	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! flagv defined only for rims 0 and 1
1771	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1772	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1773	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10.
1774	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1775	END DO
1776	END DO
1777	CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1778	ztmp(:,:) = zmask(:,:)
1779	DO ib_bdy = 1, nb_bdy
1780	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! ntreat defined only for rims 0 and 1
1781	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1782	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1783	ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10.
1784	IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1785	END DO
1786	END DO
1787	CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1788	END DO
1789	CALL iom_close( inum )
1790
1791	END SUBROUTINE bdy_meshwri
1792
1793	!!=================================================================================
1794	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.1_MIRROR_CO8_package/src/OCE/BDY/bdyini.F90 @ 12094

Download in other formats: