Context Navigation

bdyini.F90

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

Download in other formats: