Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

bdyini.F90 @ 15328

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

Download in other formats: