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/ticket2680_C1D_PAPA/src/OCE/BDY – NEMO

Context Navigation

source: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/BDY/bdyini.F90 @ 15020

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

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

Download in other formats: