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/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/BDY – NEMO

Context Navigation

source: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/BDY/bdyini.F90 @ 14576

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