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 @ 11071

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

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

Download in other formats: