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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

bdyini.F90 in NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY – NEMO

Context Navigation

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

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

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

Download in other formats: