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

Last change on this file since 11049 was 11049, checked in by girrmann, 5 years ago
dev_r10984_HPC-13 : CYCLE instruction is not systematic anymore, computation is done on the halo whenever possible and overwritten by lbc_bdy instruction, see #2285
Property svn:keywords set to `Id`
File size: 82.6 KB

Rev	Line
[1125]	1	MODULE bdyini
	2	!!======================================================================
[911]	3	!! * MODULE bdyini *
[1125]	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
[2528]	10	!! 3.3 ! 2010-09 (E.O'Dea) updates for Shelf configurations
	11	!! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions
[3294]	12	!! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge
[3651]	13	!! 3.4 ! 2012 (J. Chanut) straight open boundary case update
[6140]	14	!! 3.5 ! 2012 (S. Mocavero, I. Epicoco) optimization of BDY communications
[7646]	15	!! 3.7 ! 2016 (T. Lovato) Remove bdy macro, call here init for dta and tides
[1125]	16	!!----------------------------------------------------------------------
[6140]	17	!! bdy_init : Initialization of unstructured open boundaries
[1125]	18	!!----------------------------------------------------------------------
[6140]	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
[7646]	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
[6140]	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
[911]	31
	32	IMPLICIT NONE
	33	PRIVATE
	34
[3294]	35	PUBLIC bdy_init ! routine called in nemo_init
[911]	36
[6140]	37	INTEGER, PARAMETER :: jp_nseg = 100 !
	38	INTEGER, PARAMETER :: nrimmax = 20 ! maximum rimwidth in structured
[11048]	39	INTEGER :: nde = 1 ! domain extended in the halo to deal with bondaries
[3651]	40	! open boundary data files
	41	! Straight open boundary segment parameters:
[6140]	42	INTEGER :: nbdysege, nbdysegw, nbdysegn, nbdysegs
	43	INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge !
	44	INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw !
	45	INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn !
	46	INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs !
[1125]	47	!!----------------------------------------------------------------------
[9598]	48	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1146]	49	!! $Id$
[10068]	50	!! Software governed by the CeCILL license (see ./LICENSE)
[2528]	51	!!----------------------------------------------------------------------
[911]	52	CONTAINS
[7646]	53
[911]	54	SUBROUTINE bdy_init
	55	!!----------------------------------------------------------------------
	56	!! * ROUTINE bdy_init *
[7646]	57	!!
	58	!! ** Purpose : Initialization of the dynamics and tracer fields with
[2715]	59	!! unstructured open boundaries.
[911]	60	!!
[7646]	61	!! ** Method : Read initialization arrays (mask, indices) to identify
	62	!! an unstructured open boundary
	63	!!
	64	!! ** Input : bdy_init.nc, input file for unstructured open boundaries
	65	!!----------------------------------------------------------------------
	66	NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file, &
	67	& ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, &
	68	& cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, &
	69	& ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
[9657]	70	& cn_ice, nn_ice_dta, &
[7646]	71	& rn_ice_tem, rn_ice_sal, rn_ice_age, &
	72	& ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
	73	!
	74	INTEGER :: ios ! Local integer output status for namelist read
	75	!!----------------------------------------------------------------------
	76
	77	! ------------------------
	78	! Read namelist parameters
	79	! ------------------------
	80	REWIND( numnam_ref ) ! Namelist nambdy in reference namelist :Unstructured open boundaries
	81	READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901)
	82	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
	83	REWIND( numnam_cfg ) ! Namelist nambdy in configuration namelist :Unstructured open boundaries
	84	READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
[9168]	85	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
[7646]	86	IF(lwm) WRITE ( numond, nambdy )
	87
[9449]	88	IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE. ! forced for Agrif children
	89
[7646]	90	! -----------------------------------------
	91	! unstructured open boundaries use control
	92	! -----------------------------------------
	93	IF ( ln_bdy ) THEN
	94	IF(lwp) WRITE(numout,*)
	95	IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
	96	IF(lwp) WRITE(numout,*) '~~~~~~~~'
	97	!
	98	! Open boundaries definition (arrays and masks)
	99	CALL bdy_segs
	100	!
	101	! Open boundaries initialisation of external data arrays
	102	CALL bdy_dta_init
	103	!
	104	! Open boundaries initialisation of tidal harmonic forcing
	105	IF( ln_tide ) CALL bdytide_init
	106	!
	107	ELSE
	108	IF(lwp) WRITE(numout,*)
	109	IF(lwp) WRITE(numout,*) 'bdy_init : open boundaries not used (ln_bdy = F)'
	110	IF(lwp) WRITE(numout,*) '~~~~~~~~'
	111	!
	112	ENDIF
	113	!
	114	END SUBROUTINE bdy_init
[9019]	115
	116
[7646]	117	SUBROUTINE bdy_segs
	118	!!----------------------------------------------------------------------
	119	!! * ROUTINE bdy_init *
	120	!!
	121	!! ** Purpose : Definition of unstructured open boundaries.
	122	!!
[2715]	123	!! ** Method : Read initialization arrays (mask, indices) to identify
	124	!! an unstructured open boundary
[911]	125	!!
	126	!! ** Input : bdy_init.nc, input file for unstructured open boundaries
	127	!!----------------------------------------------------------------------
[3294]	128	INTEGER :: ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
	129	INTEGER :: icount, icountr, ibr_max, ilen1, ibm1 ! local integers
[5836]	130	INTEGER :: iwe, ies, iso, ino, inum, id_dummy ! - -
[3294]	131	INTEGER :: igrd_start, igrd_end, jpbdta ! - -
[3651]	132	INTEGER :: jpbdtau, jpbdtas ! - -
	133	INTEGER :: ib_bdy1, ib_bdy2, ib1, ib2 ! - -
[11049]	134	INTEGER :: i_offset, j_offset, inbdy, itreat ! - -
[6140]	135	INTEGER , POINTER :: nbi, nbj, nbr ! short cuts
[4292]	136	REAL(wp), POINTER, DIMENSION(:,:) :: pmask ! pointer to 2D mask fields
[3294]	137	REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars
[3651]	138	INTEGER, DIMENSION (2) :: kdimsz
[3294]	139	INTEGER, DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays
	140	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta
	141	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points
[3651]	142	CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid
[10934]	143	INTEGER :: com_east, com_west, com_south, com_north, jpk_max ! Flags for boundaries sending
[3680]	144	INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b ! Flags for boundaries receiving
	145	INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4) ! Arrays for neighbours coordinates
[9019]	146	REAL(wp), TARGET, DIMENSION(jpi,jpj) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat)
[11048]	147	REAL(wp) , DIMENSION(jpi,jpj) :: ztmp
[11049]	148	REAL(wp), POINTER :: flagu, flagv ! short cuts
[11024]	149	LOGICAL :: llnobdy, llsobdy, lleabdy, llwebdy ! local logicals
[1125]	150	!!
[6140]	151	CHARACTER(LEN=1) :: ctypebdy ! - -
	152	INTEGER :: nbdyind, nbdybeg, nbdyend
	153	!!
[3651]	154	NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
[4147]	155	INTEGER :: ios ! Local integer output status for namelist read
[911]	156	!!----------------------------------------------------------------------
[6140]	157	!
	158	cgrid = (/'t','u','v'/)
[911]	159
[3294]	160	! -----------------------------------------
	161	! Check and write out namelist parameters
	162	! -----------------------------------------
[7646]	163	IF( jperio /= 0 ) CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,', &
	164	& ' and general open boundary condition are not compatible' )
[911]	165
[6140]	166	IF( nb_bdy == 0 ) THEN
[3294]	167	IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
[911]	168	ELSE
[6140]	169	IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy
[911]	170	ENDIF
	171
[3294]	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
[2528]	181	IF(lwp) WRITE(numout,*)
[1125]	182
[3294]	183	IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: '
[4292]	184	SELECT CASE( cn_dyn2d(ib_bdy) )
[6140]	185	CASE( 'none' )
[4292]	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.
[6140]	190	CASE( 'frs' )
[4292]	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.
[6140]	195	CASE( 'flather' )
[4292]	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.
[6140]	200	CASE( 'orlanski' )
[4292]	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.
[6140]	205	CASE( 'orlanski_npo' )
[4292]	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' )
[3294]	211	END SELECT
[4292]	212	IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN
[3294]	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
[7646]	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' )
[3651]	222	ENDIF
[3294]	223	ENDIF
[2528]	224	IF(lwp) WRITE(numout,*)
[911]	225
[3294]	226	IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: '
[4292]	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.
[7646]	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.
[4292]	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' )
[3294]	261	END SELECT
[4292]	262	IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
[3294]	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
[3651]	269
	270	IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
[4292]	271	IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN
[3651]	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.
[4292]	274	ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN
[3651]	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' )
[4292]	280	dta_bdy(ib_bdy)%ll_u3d = .true.
	281	dta_bdy(ib_bdy)%ll_v3d = .true.
[3651]	282	ENDIF
	283	ELSE
	284	IF(lwp) WRITE(numout,*) ' NO relaxation on baroclinic velocities'
	285	ENDIF
[2528]	286	IF(lwp) WRITE(numout,*)
[1125]	287
[3294]	288	IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: '
[4292]	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' )
[3294]	319	END SELECT
[4292]	320	IF( cn_tra(ib_bdy) /= 'none' ) THEN
[3294]	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
[3651]	327
	328	IF ( ln_tra_dmp(ib_bdy) ) THEN
[4292]	329	IF ( cn_tra(ib_bdy) == 'none' ) THEN
[3651]	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.
[4292]	332	ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN
[3651]	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'
[4292]	337	IF(lwp) WRITE(numout,*) ' Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'
[3651]	338	IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
[4292]	339	dta_bdy(ib_bdy)%ll_tem = .true.
	340	dta_bdy(ib_bdy)%ll_sal = .true.
[3651]	341	ENDIF
	342	ELSE
	343	IF(lwp) WRITE(numout,*) ' NO T/S relaxation'
	344	ENDIF
[2528]	345	IF(lwp) WRITE(numout,*)
[1125]	346
[9570]	347	#if defined key_si3
[9019]	348	IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: '
[9657]	349	SELECT CASE( cn_ice(ib_bdy) )
[9019]	350	CASE('none')
[4292]	351	IF(lwp) WRITE(numout,*) ' no open boundary condition'
[9019]	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')
[4292]	356	IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme'
[9019]	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.
[9657]	360	CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_ice' )
[9019]	361	END SELECT
[9657]	362	IF( cn_ice(ib_bdy) /= 'none' ) THEN
	363	SELECT CASE( nn_ice_dta(ib_bdy) ) !
[3294]	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'
[9657]	366	CASE DEFAULT ; CALL ctl_stop( 'nn_ice_dta must be 0 or 1' )
[3294]	367	END SELECT
	368	ENDIF
[2528]	369	IF(lwp) WRITE(numout,*)
[4689]	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)
[3294]	373	#endif
[911]	374
[3651]	375	IF(lwp) WRITE(numout,*) ' Width of relaxation zone = ', nn_rimwidth(ib_bdy)
[3294]	376	IF(lwp) WRITE(numout,*)
[9019]	377	!
	378	END DO
[2528]	379
[9019]	380	IF( nb_bdy > 0 ) THEN
[3651]	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,*)
[10481]	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
[3651]	399	ELSE
	400	IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
	401	IF(lwp) WRITE(numout,*)
	402	ENDIF
[10934]	403	IF( nb_jpk_bdy(ib_bdy) > 0 ) THEN
[7646]	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
[3294]	408	ENDIF
	409
[1125]	410	! -------------------------------------------------
[3294]	411	! Initialise indices arrays for open boundaries
	412	! -------------------------------------------------
[911]	413
[3294]	414	! Work out global dimensions of boundary data
	415	! ---------------------------------------------
[4147]	416	REWIND( numnam_cfg )
	417
[3651]	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
[3294]	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
[3651]	432	icount = icount + 1
[3294]	433	! No REWIND here because may need to read more than one nambdy_index namelist.
[9168]	434	! Read only namelist_cfg to avoid unseccessfull overwrite
	435	! keep full control of the configuration namelist
[4147]	436	READ ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 )
[9169]	437	904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp )
[4624]	438	IF(lwm) WRITE ( numond, nambdy_index )
[4147]	439
[3651]	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
[3294]	491
[3651]	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)
[3294]	496
[3651]	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' )
[3294]	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 )
[4333]	506	nblendta(igrd,ib_bdy) = MAXVAL(kdimsz)
	507	jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz))
[6140]	508	END DO
[3651]	509	CALL iom_close( inum )
[6140]	510	!
[3294]	511	ENDIF
[6140]	512	!
	513	END DO ! ib_bdy
[3294]	514
[3651]	515	IF (nb_bdy>0) THEN
	516	jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
[3294]	517
[3651]	518	! Allocate arrays
	519	!---------------
	520	ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), &
	521	& nbrdta(jpbdta, jpbgrd, nb_bdy) )
[10934]	522
	523	jpk_max = MAXVAL(nb_jpk_bdy)
	524	jpk_max = MAX(jpk_max, jpk)
[3294]	525
[10934]	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
[7646]	529
	530	IF ( icount>0 ) THEN
[10934]	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
[7646]	534	ENDIF
[3651]	535	!
	536	ENDIF
	537
	538	! Now look for crossings in user (namelist) defined open boundary segments:
	539	!--------------------------------------------------------------------------
[6140]	540	IF( icount>0 ) CALL bdy_ctl_seg
[3651]	541
[3294]	542	! Calculate global boundary index arrays or read in from file
[3651]	543	!------------------------------------------------------------
	544	! 1. Read global index arrays from boundary coordinates file.
[3294]	545	DO ib_bdy = 1, nb_bdy
[6140]	546	!
[3651]	547	IF( ln_coords_file(ib_bdy) ) THEN
[6140]	548	!
[3651]	549	CALL iom_open( cn_coords_file(ib_bdy), inum )
[3294]	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
[6140]	563	!
[3294]	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 )
[6140]	572	!
[3294]	573	ENDIF
[6140]	574	!
	575	END DO
[3651]	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
[3294]	757
[3651]	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
[3294]	794	! Work out dimensions of boundary data on each processor
	795	! ------------------------------------------------------
	796
[3651]	797	! Rather assume that boundary data indices are given on global domain
	798	! TO BE DISCUSSED ?
	799	! iw = mig(1) + 1 ! if monotasking and no zoom, iw=2
	800	! ie = mig(1) + nlci-1 - 1 ! if monotasking and no zoom, ie=jpim1
	801	! is = mjg(1) + 1 ! if monotasking and no zoom, is=2
	802	! in = mjg(1) + nlcj-1 - 1 ! if monotasking and no zoom, in=jpjm1
[11048]	803	iwe = mig(1) - 1 + 2 - nde ! if monotasking and no zoom, iw=2
	804	ies = mig(1) + nlci-1 - 1 + nde ! if monotasking and no zoom, ie=jpim1
[11049]	805	iso = mjg(1) - 1 + 2 - nde ! if monotasking and no zoom, is=2
[11048]	806	ino = mjg(1) + nlcj-1 - 1 + nde ! if monotasking and no zoom, in=jpjm1
[3651]	807
[3680]	808	ALLOCATE( nbondi_bdy(nb_bdy))
	809	ALLOCATE( nbondj_bdy(nb_bdy))
	810	nbondi_bdy(:)=2
	811	nbondj_bdy(:)=2
	812	ALLOCATE( nbondi_bdy_b(nb_bdy))
	813	ALLOCATE( nbondj_bdy_b(nb_bdy))
	814	nbondi_bdy_b(:)=2
	815	nbondj_bdy_b(:)=2
	816
	817	! Work out dimensions of boundary data on each neighbour process
[6140]	818	IF(nbondi == 0) THEN
[7646]	819	iw_b(1) = 1 + nimppt(nowe+1)
	820	ie_b(1) = 1 + nimppt(nowe+1)+nlcit(nowe+1)-3
	821	is_b(1) = 1 + njmppt(nowe+1)
	822	in_b(1) = 1 + njmppt(nowe+1)+nlcjt(nowe+1)-3
[3680]	823
[7646]	824	iw_b(2) = 1 + nimppt(noea+1)
	825	ie_b(2) = 1 + nimppt(noea+1)+nlcit(noea+1)-3
	826	is_b(2) = 1 + njmppt(noea+1)
	827	in_b(2) = 1 + njmppt(noea+1)+nlcjt(noea+1)-3
[6140]	828	ELSEIF(nbondi == 1) THEN
[7646]	829	iw_b(1) = 1 + nimppt(nowe+1)
	830	ie_b(1) = 1 + nimppt(nowe+1)+nlcit(nowe+1)-3
	831	is_b(1) = 1 + njmppt(nowe+1)
	832	in_b(1) = 1 + njmppt(nowe+1)+nlcjt(nowe+1)-3
[6140]	833	ELSEIF(nbondi == -1) THEN
[7646]	834	iw_b(2) = 1 + nimppt(noea+1)
	835	ie_b(2) = 1 + nimppt(noea+1)+nlcit(noea+1)-3
	836	is_b(2) = 1 + njmppt(noea+1)
	837	in_b(2) = 1 + njmppt(noea+1)+nlcjt(noea+1)-3
[3680]	838	ENDIF
	839
[6140]	840	IF(nbondj == 0) THEN
[7646]	841	iw_b(3) = 1 + nimppt(noso+1)
	842	ie_b(3) = 1 + nimppt(noso+1)+nlcit(noso+1)-3
	843	is_b(3) = 1 + njmppt(noso+1)
	844	in_b(3) = 1 + njmppt(noso+1)+nlcjt(noso+1)-3
[3680]	845
[7646]	846	iw_b(4) = 1 + nimppt(nono+1)
	847	ie_b(4) = 1 + nimppt(nono+1)+nlcit(nono+1)-3
	848	is_b(4) = 1 + njmppt(nono+1)
	849	in_b(4) = 1 + njmppt(nono+1)+nlcjt(nono+1)-3
[6140]	850	ELSEIF(nbondj == 1) THEN
[7646]	851	iw_b(3) = 1 + nimppt(noso+1)
	852	ie_b(3) = 1 + nimppt(noso+1)+nlcit(noso+1)-3
	853	is_b(3) = 1 + njmppt(noso+1)
	854	in_b(3) = 1 + njmppt(noso+1)+nlcjt(noso+1)-3
[6140]	855	ELSEIF(nbondj == -1) THEN
[7646]	856	iw_b(4) = 1 + nimppt(nono+1)
	857	ie_b(4) = 1 + nimppt(nono+1)+nlcit(nono+1)-3
	858	is_b(4) = 1 + njmppt(nono+1)
	859	in_b(4) = 1 + njmppt(nono+1)+nlcjt(nono+1)-3
[3680]	860	ENDIF
	861
[3294]	862	DO ib_bdy = 1, nb_bdy
	863	DO igrd = 1, jpbgrd
	864	icount = 0
	865	icountr = 0
	866	idx_bdy(ib_bdy)%nblen(igrd) = 0
	867	idx_bdy(ib_bdy)%nblenrim(igrd) = 0
	868	DO ib = 1, nblendta(igrd,ib_bdy)
	869	! check that data is in correct order in file
	870	ibm1 = MAX(1,ib-1)
	871	IF(lwp) THEN ! Since all procs read global data only need to do this check on one proc...
	872	IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
[7646]	873	CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', &
[4294]	874	& ' in order of distance from edge nbr A utility for re-ordering ', &
	875	& ' boundary coordinates and data files exists in the TOOLS/OBC directory')
[3294]	876	ENDIF
	877	ENDIF
	878	! check if point is in local domain
[5656]	879	IF( nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND. &
	880	& nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino ) THEN
[3294]	881	!
	882	icount = icount + 1
	883	!
	884	IF( nbrdta(ib,igrd,ib_bdy) == 1 ) icountr = icountr+1
	885	ENDIF
[9019]	886	END DO
[3294]	887	idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
	888	idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc
[9019]	889	END DO ! igrd
[3294]	890
	891	! Allocate index arrays for this boundary set
	892	!--------------------------------------------
[6140]	893	ilen1 = MAXVAL( idx_bdy(ib_bdy)%nblen(:) )
[9019]	894	ALLOCATE( idx_bdy(ib_bdy)%nbi (ilen1,jpbgrd) , &
	895	& idx_bdy(ib_bdy)%nbj (ilen1,jpbgrd) , &
	896	& idx_bdy(ib_bdy)%nbr (ilen1,jpbgrd) , &
	897	& idx_bdy(ib_bdy)%nbd (ilen1,jpbgrd) , &
	898	& idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) , &
[11024]	899	& idx_bdy(ib_bdy)%ntreat(ilen1,jpbgrd) , &
[9019]	900	& idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) , &
	901	& idx_bdy(ib_bdy)%nbw (ilen1,jpbgrd) , &
	902	& idx_bdy(ib_bdy)%flagu (ilen1,jpbgrd) , &
	903	& idx_bdy(ib_bdy)%flagv (ilen1,jpbgrd) )
[3294]	904
	905	! Dispatch mapping indices and discrete distances on each processor
	906	! -----------------------------------------------------------------
	907
[6140]	908	com_east = 0
	909	com_west = 0
[3680]	910	com_south = 0
	911	com_north = 0
	912
[6140]	913	com_east_b = 0
	914	com_west_b = 0
[3680]	915	com_south_b = 0
	916	com_north_b = 0
[5656]	917
[3294]	918	DO igrd = 1, jpbgrd
	919	icount = 0
	920	! Loop on rimwidth to ensure outermost points come first in the local arrays.
	921	DO ir=1, nn_rimwidth(ib_bdy)
	922	DO ib = 1, nblendta(igrd,ib_bdy)
	923	! check if point is in local domain and equals ir
[5656]	924	IF( nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND. &
	925	& nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino .AND. &
[3294]	926	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
	927	!
	928	icount = icount + 1
[3651]	929
	930	! Rather assume that boundary data indices are given on global domain
	931	! TO BE DISCUSSED ?
	932	! idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1
	933	! idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
[7646]	934	idx_bdy(ib_bdy)%nbi(icount,igrd) = nbidta(ib,igrd,ib_bdy)- mig(1)+1
	935	idx_bdy(ib_bdy)%nbj(icount,igrd) = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
[3680]	936	! check if point has to be sent
	937	ii = idx_bdy(ib_bdy)%nbi(icount,igrd)
	938	ij = idx_bdy(ib_bdy)%nbj(icount,igrd)
[10629]	939	if((com_east .ne. 1) .and. (ii == (nlci-1)) .and. (nbondi .le. 0)) then
[3680]	940	com_east = 1
[6140]	941	elseif((com_west .ne. 1) .and. (ii == 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then
[3680]	942	com_west = 1
	943	endif
[6140]	944	if((com_south .ne. 1) .and. (ij == 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then
[3680]	945	com_south = 1
[10629]	946	elseif((com_north .ne. 1) .and. (ij == (nlcj-1)) .and. (nbondj .le. 0)) then
[3680]	947	com_north = 1
	948	endif
[3294]	949	idx_bdy(ib_bdy)%nbr(icount,igrd) = nbrdta(ib,igrd,ib_bdy)
	950	idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
	951	ENDIF
[3680]	952	! check if point has to be received from a neighbour
[6140]	953	IF(nbondi == 0) THEN
[3680]	954	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND. &
	955	& nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND. &
	956	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
	957	ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
[10983]	958	if( ii == (nlcit(nowe+1)-1) ) then
[3680]	959	ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
[6140]	960	if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
[3680]	961	com_south = 1
[6140]	962	elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
[3680]	963	com_north = 1
	964	endif
	965	com_west_b = 1
	966	endif
	967	ENDIF
	968	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND. &
	969	& nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND. &
	970	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
	971	ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
[10983]	972	if( ii == 2 ) then
[3680]	973	ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
[6140]	974	if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
[3680]	975	com_south = 1
[6140]	976	elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
[3680]	977	com_north = 1
	978	endif
	979	com_east_b = 1
	980	endif
	981	ENDIF
[6140]	982	ELSEIF(nbondi == 1) THEN
[3680]	983	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND. &
	984	& nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND. &
	985	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
	986	ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
[10983]	987	if( ii == (nlcit(nowe+1)-1) ) then
[3680]	988	ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
[6140]	989	if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
[3680]	990	com_south = 1
[6140]	991	elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
[3680]	992	com_north = 1
	993	endif
	994	com_west_b = 1
	995	endif
	996	ENDIF
[6140]	997	ELSEIF(nbondi == -1) THEN
[3680]	998	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND. &
	999	& nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND. &
	1000	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
	1001	ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
[10983]	1002	if( ii == 2 ) then
[3680]	1003	ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
[6140]	1004	if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
[3680]	1005	com_south = 1
[6140]	1006	elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
[3680]	1007	com_north = 1
	1008	endif
	1009	com_east_b = 1
	1010	endif
	1011	ENDIF
	1012	ENDIF
[6140]	1013	IF(nbondj == 0) THEN
[3703]	1014	IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 &
	1015	& .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
[3680]	1016	& nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
	1017	com_north_b = 1
	1018	ENDIF
[3703]	1019	IF(com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 &
	1020	&.OR. nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
[3680]	1021	& nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
	1022	com_south_b = 1
	1023	ENDIF
	1024	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND. &
	1025	& nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND. &
	1026	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
	1027	ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
[10629]	1028	if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then
[3680]	1029	com_south_b = 1
	1030	endif
	1031	ENDIF
	1032	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND. &
	1033	& nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND. &
	1034	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
	1035	ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
[6140]	1036	if((com_north_b .ne. 1) .and. (ij == 2)) then
[3680]	1037	com_north_b = 1
	1038	endif
	1039	ENDIF
[6140]	1040	ELSEIF(nbondj == 1) THEN
[3703]	1041	IF( com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. &
	1042	& nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
[3680]	1043	& nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
	1044	com_south_b = 1
	1045	ENDIF
	1046	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND. &
	1047	& nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND. &
	1048	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
	1049	ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
[10629]	1050	if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then
[3680]	1051	com_south_b = 1
	1052	endif
	1053	ENDIF
[6140]	1054	ELSEIF(nbondj == -1) THEN
[3703]	1055	IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 &
	1056	& .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
[3680]	1057	& nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
	1058	com_north_b = 1
	1059	ENDIF
	1060	IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND. &
	1061	& nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND. &
	1062	& nbrdta(ib,igrd,ib_bdy) == ir ) THEN
	1063	ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
[6140]	1064	if((com_north_b .ne. 1) .and. (ij == 2)) then
[3680]	1065	com_north_b = 1
	1066	endif
	1067	ENDIF
	1068	ENDIF
[3294]	1069	ENDDO
	1070	ENDDO
	1071	ENDDO
[4292]	1072
[6140]	1073	! definition of the i- and j- direction local boundaries arrays used for sending the boundaries
	1074	IF( (com_east == 1) .and. (com_west == 1) ) THEN ; nbondi_bdy(ib_bdy) = 0
	1075	ELSEIF( (com_east == 1) .and. (com_west == 0) ) THEN ; nbondi_bdy(ib_bdy) = -1
	1076	ELSEIF( (com_east == 0) .and. (com_west == 1) ) THEN ; nbondi_bdy(ib_bdy) = 1
[3680]	1077	ENDIF
[6140]	1078	IF( (com_north == 1) .and. (com_south == 1) ) THEN ; nbondj_bdy(ib_bdy) = 0
	1079	ELSEIF( (com_north == 1) .and. (com_south == 0) ) THEN ; nbondj_bdy(ib_bdy) = -1
	1080	ELSEIF( (com_north == 0) .and. (com_south == 1) ) THEN ; nbondj_bdy(ib_bdy) = 1
[3680]	1081	ENDIF
	1082
[6140]	1083	! definition of the i- and j- direction local boundaries arrays used for receiving the boundaries
	1084	IF( (com_east_b == 1) .and. (com_west_b == 1) ) THEN ; nbondi_bdy_b(ib_bdy) = 0
	1085	ELSEIF( (com_east_b == 1) .and. (com_west_b == 0) ) THEN ; nbondi_bdy_b(ib_bdy) = -1
	1086	ELSEIF( (com_east_b == 0) .and. (com_west_b == 1) ) THEN ; nbondi_bdy_b(ib_bdy) = 1
[3680]	1087	ENDIF
[6140]	1088	IF( (com_north_b == 1) .and. (com_south_b == 1) ) THEN ; nbondj_bdy_b(ib_bdy) = 0
	1089	ELSEIF( (com_north_b == 1) .and. (com_south_b == 0) ) THEN ; nbondj_bdy_b(ib_bdy) = -1
	1090	ELSEIF( (com_north_b == 0) .and. (com_south_b == 1) ) THEN ; nbondj_bdy_b(ib_bdy) = 1
[3680]	1091	ENDIF
	1092
[3294]	1093	! Compute rim weights for FRS scheme
	1094	! ----------------------------------
	1095	DO igrd = 1, jpbgrd
	1096	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
	1097	nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
[6140]	1098	idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( nbr - 1 ) *0.5 ) ! tanh formulation
	1099	! idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic
	1100	! idx_bdy(ib_bdy)%nbw(ib,igrd) = REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)) ! linear
[3294]	1101	END DO
	1102	END DO
	1103
[3651]	1104	! Compute damping coefficients
	1105	! ----------------------------
	1106	DO igrd = 1, jpbgrd
	1107	DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
	1108	nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
	1109	idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &
[6140]	1110	& (REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))*2. ! quadratic
[4292]	1111	idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) &
[6140]	1112	& (REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))*2. ! quadratic
[3651]	1113	END DO
	1114	END DO
	1115
[9019]	1116	END DO
[3294]	1117
	1118	! ------------------------------------------------------
	1119	! Initialise masks and find normal/tangential directions
	1120	! ------------------------------------------------------
	1121
[1125]	1122	! Read global 2D mask at T-points: bdytmask
[3294]	1123	! -----------------------------------------
[1125]	1124	! bdytmask = 1 on the computational domain AND on open boundaries
	1125	! = 0 elsewhere
[911]	1126
[7646]	1127	bdytmask(:,:) = ssmask(:,:)
	1128
[9600]	1129	! Derive mask on U and V grid from mask on T grid
[911]	1130
[9600]	1131	bdyumask(:,:) = 0._wp
	1132	bdyvmask(:,:) = 0._wp
	1133	DO ij = 1, jpjm1
	1134	DO ii = 1, jpim1
	1135	bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1, ij )
	1136	bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii ,ij+1)
[1125]	1137	END DO
[9600]	1138	END DO
[10425]	1139	CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1. , bdyvmask, 'V', 1. ) ! Lateral boundary cond.
[911]	1140
[6125]	1141	! bdy masks are now set to zero on boundary points:
[1125]	1142	!
	1143	igrd = 1
[3294]	1144	DO ib_bdy = 1, nb_bdy
	1145	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
[6140]	1146	bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
	1147	END DO
	1148	END DO
[1125]	1149	!
	1150	igrd = 2
[3294]	1151	DO ib_bdy = 1, nb_bdy
	1152	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
[6140]	1153	bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
	1154	END DO
	1155	END DO
[1125]	1156	!
	1157	igrd = 3
[3294]	1158	DO ib_bdy = 1, nb_bdy
	1159	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
[6140]	1160	bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
[9019]	1161	END DO
	1162	END DO
[911]	1163
[4292]	1164	! For the flagu/flagv calculation below we require a version of fmask without
	1165	! the land boundary condition (shlat) included:
[9807]	1166	zfmask(:,:) = 0
[4292]	1167	DO ij = 2, jpjm1
	1168	DO ii = 2, jpim1
	1169	zfmask(ii,ij) = tmask(ii,ij ,1) * tmask(ii+1,ij ,1) &
	1170	& * tmask(ii,ij+1,1) * tmask(ii+1,ij+1,1)
	1171	END DO
	1172	END DO
	1173
[1125]	1174	! Lateral boundary conditions
[10425]	1175	CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )
	1176	CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1. , bdyvmask, 'V', 1., bdytmask, 'T', 1. )
[3294]	1177	DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components
	1178	icount = 0
	1179
[4292]	1180	! Calculate relationship of U direction to the local orientation of the boundary
	1181	! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
	1182	! flagu = 0 : u is tangential
	1183	! flagu = 1 : u is normal to the boundary and is direction is inward
[3294]	1184
[9019]	1185	DO igrd = 1, jpbgrd
[4292]	1186	SELECT CASE( igrd )
[6140]	1187	CASE( 1 ) ; pmask => umask (:,:,1) ; i_offset = 0
	1188	CASE( 2 ) ; pmask => bdytmask(:,:) ; i_offset = 1
	1189	CASE( 3 ) ; pmask => zfmask (:,:) ; i_offset = 0
[4292]	1190	END SELECT
	1191	icount = 0
[11048]	1192	ztmp(:,:) = 0._wp
[4292]	1193	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
[11048]	1194	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	1195	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
	1196	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE
	1197	zefl = pmask(ii+i_offset-1,ij)
	1198	zwfl = pmask(ii+i_offset ,ij)
[4292]	1199	! This error check only works if you are using the bdyXmask arrays
[11048]	1200	IF( i_offset == 1 .and. zefl + zwfl == 2. ) THEN
[4292]	1201	icount = icount + 1
[11048]	1202	IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
[4292]	1203	ELSE
[11048]	1204	ztmp(ii,ij) = -zefl + zwfl
[4292]	1205	ENDIF
	1206	END DO
	1207	IF( icount /= 0 ) THEN
[10425]	1208	WRITE(ctmp1,*) ' E R R O R : Some ',cgrid(igrd),' grid points,', &
[4292]	1209	' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
[10425]	1210	WRITE(ctmp2,*) ' ========== '
	1211	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[4292]	1212	ENDIF
[11048]	1213	CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )
	1214	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
	1215	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	1216	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
	1217	idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij)
	1218	END DO
[1125]	1219	END DO
[911]	1220
[4292]	1221	! Calculate relationship of V direction to the local orientation of the boundary
	1222	! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
	1223	! flagv = 0 : v is tangential
	1224	! flagv = 1 : v is normal to the boundary and is direction is inward
[911]	1225
[6140]	1226	DO igrd = 1, jpbgrd
[4292]	1227	SELECT CASE( igrd )
[6140]	1228	CASE( 1 ) ; pmask => vmask (:,:,1) ; j_offset = 0
	1229	CASE( 2 ) ; pmask => zfmask(:,:) ; j_offset = 0
	1230	CASE( 3 ) ; pmask => bdytmask ; j_offset = 1
[4292]	1231	END SELECT
	1232	icount = 0
[11048]	1233	ztmp(:,:) = 0._wp
[4292]	1234	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
[11048]	1235	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	1236	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
	1237	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE
	1238	znfl = pmask(ii,ij+j_offset-1)
	1239	zsfl = pmask(ii,ij+j_offset )
[4292]	1240	! This error check only works if you are using the bdyXmask arrays
[11048]	1241	IF( j_offset == 1 .and. znfl + zsfl == 2. ) THEN
	1242	IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
[4292]	1243	icount = icount + 1
	1244	ELSE
[11048]	1245	ztmp(ii,ij) = -znfl + zsfl
[4292]	1246	END IF
	1247	END DO
	1248	IF( icount /= 0 ) THEN
[10425]	1249	WRITE(ctmp1,*) ' E R R O R : Some ',cgrid(igrd),' grid points,', &
[4292]	1250	' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
[10425]	1251	WRITE(ctmp2,*) ' ========== '
	1252	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[11048]	1253	ENDIF
	1254	CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )
	1255	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
	1256	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	1257	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
	1258	idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij)
	1259	END DO
[1125]	1260	END DO
[6140]	1261	!
[4292]	1262	END DO
[11024]	1263
[11044]	1264	! detect corner interior and its orientation index 1 to 4 depending on the orientation
	1265	! detect corner exterior and its orientation index 5 to 8 depending on the orientation
	1266	! detect geometries with 3 neighbours index 9 to 12 depending on the orientation
	1267	! else index 0
[11024]	1268	DO ib_bdy = 1, nb_bdy
	1269	DO igrd = 1, jpbgrd
	1270	SELECT CASE( igrd )
	1271	CASE( 1 ) ; pmask => bdytmask
	1272	CASE( 2 ) ; pmask => bdyumask
	1273	CASE( 3 ) ; pmask => bdyvmask
	1274	END SELECT
[11048]	1275	ztmp(:,:) = 0._wp
[11024]	1276	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
[11049]	1277	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	1278	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
[11048]	1279	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE
[11024]	1280	llnobdy = pmask(ii ,ij+1) == 1.
	1281	llsobdy = pmask(ii ,ij-1) == 1.
	1282	lleabdy = pmask(ii+1,ij ) == 1.
	1283	llwebdy = pmask(ii-1,ij ) == 1.
	1284	inbdy = COUNT( (/ llnobdy, llsobdy, lleabdy, llwebdy /) )
	1285	IF( inbdy == 0 ) THEN ! no neighbours -> interior of a corner
	1286	! ! ! _____ ! _____
	1287	! 1 \| o ! 2 o \| ! 3 \| x ! 4 x \|
[11048]	1288	! \|_x_ _ ! _ _x_\| ! \| o ! o \|
	1289	IF( pmask(ii+1,ij+1) == 1. ) ztmp(ii,ij) = 1
	1290	IF( pmask(ii-1,ij+1) == 1. ) ztmp(ii,ij) = 2
	1291	IF( pmask(ii+1,ij-1) == 1. ) ztmp(ii,ij) = 3
	1292	IF( pmask(ii-1,ij-1) == 1. ) ztmp(ii,ij) = 4
[11024]	1293	END IF
	1294	IF( inbdy == 1 ) THEN ! middle of linear bdy
[11048]	1295	ztmp(ii,ij) = 0 ! regular treatment with flags
[11024]	1296	END IF
	1297	IF( inbdy == 2 ) THEN ! exterior of a corner
[11044]	1298	! o ! o ! _____\| ! \|_____
	1299	! 5 ____x o ! 6 o x___ ! 7 x o ! 8 o x
	1300	! \| ! \| ! o ! o
[11048]	1301	IF( llnobdy .AND. lleabdy ) ztmp(ii,ij) = 5
	1302	IF( llnobdy .AND. llwebdy ) ztmp(ii,ij) = 6
	1303	IF( llsobdy .AND. lleabdy ) ztmp(ii,ij) = 7
	1304	IF( llsobdy .AND. llwebdy ) ztmp(ii,ij) = 8
[11024]	1305	END IF
	1306	IF( inbdy == 3 ) THEN ! 3 neighbours __ __
	1307	! \|_ o ! o _\| ! \|_\| ! o
[11044]	1308	! 9 _\| x o ! 10 o x \|_ ! 11 o x o ! 12 o x o
	1309	! \| o ! o \| ! o ! __\|¨\|__
[11048]	1310	IF( llnobdy .AND. lleabdy .AND. llsobdy ) ztmp(ii,ij) = 9
	1311	IF( llnobdy .AND. llwebdy .AND. llsobdy ) ztmp(ii,ij) = 10
	1312	IF( llwebdy .AND. llsobdy .AND. lleabdy ) ztmp(ii,ij) = 11
	1313	IF( llwebdy .AND. llnobdy .AND. lleabdy ) ztmp(ii,ij) = 12
[11024]	1314	END IF
	1315	IF( inbdy == 4 ) THEN
	1316	WRITE(ctmp1,*) ' E R R O R : Problem with ',cgrid(igrd) ,' grid points,', &
	1317	' some points on boundary set ', ib_bdy, ' have 4 neighbours'
	1318	WRITE(ctmp2,*) ' ========== '
	1319	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
	1320	END IF
	1321	END DO
[11048]	1322	CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )
	1323	DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
	1324	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	1325	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
	1326	idx_bdy(ib_bdy)%ntreat(ib,igrd) = ztmp(ii,ij)
	1327	END DO
[11024]	1328	END DO
	1329	END DO
[11049]	1330
[3294]	1331	!
	1332	! Tidy up
	1333	!--------
[6140]	1334	IF( nb_bdy>0 ) DEALLOCATE( nbidta, nbjdta, nbrdta )
	1335	!
[7646]	1336	END SUBROUTINE bdy_segs
[911]	1337
[3651]	1338	SUBROUTINE bdy_ctl_seg
	1339	!!----------------------------------------------------------------------
	1340	!! * ROUTINE bdy_ctl_seg *
	1341	!!
	1342	!! ** Purpose : Check straight open boundary segments location
	1343	!!
	1344	!! ** Method : - Look for open boundary corners
	1345	!! - Check that segments start or end on land
	1346	!!----------------------------------------------------------------------
	1347	INTEGER :: ib, ib1, ib2, ji ,jj, itest
	1348	INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns
	1349	REAL(wp), DIMENSION(2) :: ztestmask
	1350	!!----------------------------------------------------------------------
	1351	!
	1352	IF (lwp) WRITE(numout,*) ' '
	1353	IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
	1354	IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
	1355	!
	1356	IF(lwp) WRITE(numout,*) 'Number of east segments : ', nbdysege
	1357	IF(lwp) WRITE(numout,*) 'Number of west segments : ', nbdysegw
	1358	IF(lwp) WRITE(numout,*) 'Number of north segments : ', nbdysegn
	1359	IF(lwp) WRITE(numout,*) 'Number of south segments : ', nbdysegs
	1360	! 1. Check bounds
	1361	!----------------
	1362	DO ib = 1, nbdysegn
	1363	IF (lwp) WRITE(numout,) '*check north seg bounds pckg: ', npckgn(ib)
	1364	IF ((jpjnob(ib).ge.jpjglo-1).or.&
	1365	&(jpjnob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
	1366	IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
	1367	IF (jpindt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
	1368	IF (jpinft(ib).ge.jpiglo) CALL ctl_stop( 'End index out of domain' )
	1369	END DO
	1370	!
	1371	DO ib = 1, nbdysegs
	1372	IF (lwp) WRITE(numout,) '*check south seg bounds pckg: ', npckgs(ib)
	1373	IF ((jpjsob(ib).ge.jpjglo-1).or.&
	1374	&(jpjsob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
	1375	IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
	1376	IF (jpisdt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
	1377	IF (jpisft(ib).ge.jpiglo) CALL ctl_stop( 'End index out of domain' )
	1378	END DO
	1379	!
	1380	DO ib = 1, nbdysege
	1381	IF (lwp) WRITE(numout,) '*check east seg bounds pckg: ', npckge(ib)
	1382	IF ((jpieob(ib).ge.jpiglo-1).or.&
	1383	&(jpieob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
	1384	IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
	1385	IF (jpjedt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
	1386	IF (jpjeft(ib).ge.jpjglo) CALL ctl_stop( 'End index out of domain' )
	1387	END DO
	1388	!
	1389	DO ib = 1, nbdysegw
	1390	IF (lwp) WRITE(numout,) '*check west seg bounds pckg: ', npckgw(ib)
	1391	IF ((jpiwob(ib).ge.jpiglo-1).or.&
	1392	&(jpiwob(ib).le.1)) CALL ctl_stop( 'nbdyind out of domain' )
	1393	IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
	1394	IF (jpjwdt(ib).le.1 ) CALL ctl_stop( 'Start index out of domain' )
	1395	IF (jpjwft(ib).ge.jpjglo) CALL ctl_stop( 'End index out of domain' )
	1396	ENDDO
	1397	!
	1398	!
	1399	! 2. Look for segment crossings
	1400	!------------------------------
	1401	IF (lwp) WRITE(numout,) '*Look for segments corners :'
	1402	!
	1403	itest = 0 ! corner number
	1404	!
	1405	! flag to detect if start or end of open boundary belongs to a corner
	1406	! if not (=0), it must be on land.
	1407	! if a corner is detected, save bdy package number for further tests
	1408	icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
	1409	! South/West crossings
	1410	IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
	1411	DO ib1 = 1, nbdysegw
	1412	DO ib2 = 1, nbdysegs
	1413	IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
	1414	& ( jpisft(ib2)>=jpiwob(ib1)).AND. &
	1415	& ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
	1416	& ( jpjwft(ib1)>=jpjsob(ib2))) THEN
	1417	IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN
	1418	! We have a possible South-West corner
	1419	! WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
	1420	! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
	1421	icornw(ib1,1) = npckgs(ib2)
	1422	icorns(ib2,1) = npckgw(ib1)
	1423	ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
[10425]	1424	WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
	1425	& jpisft(ib2), jpjwft(ib1)
	1426	WRITE(ctmp2,*) ' ========== Not allowed yet'
	1427	WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
	1428	& ' and South segment: ',npckgs(ib2)
	1429	CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' )
[3651]	1430	ELSE
[10425]	1431	WRITE(ctmp1,*) ' E R R O R : Check South and West Open boundary indices'
	1432	WRITE(ctmp2,*) ' ========== Crossing problem with West segment: ',npckgw(ib1) , &
	1433	& ' and South segment: ',npckgs(ib2)
	1434	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1435	END IF
	1436	END IF
	1437	END DO
	1438	END DO
	1439	END IF
	1440	!
	1441	! South/East crossings
	1442	IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
	1443	DO ib1 = 1, nbdysege
	1444	DO ib2 = 1, nbdysegs
	1445	IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
	1446	& ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
	1447	& ( jpjedt(ib1)<=jpjsob(ib2) ).AND. &
	1448	& ( jpjeft(ib1)>=jpjsob(ib2) )) THEN
	1449	IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
	1450	! We have a possible South-East corner
	1451	! WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
	1452	! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
	1453	icorne(ib1,1) = npckgs(ib2)
	1454	icorns(ib2,2) = npckge(ib1)
	1455	ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
[10425]	1456	WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
	1457	& jpisdt(ib2), jpjeft(ib1)
	1458	WRITE(ctmp2,*) ' ========== Not allowed yet'
	1459	WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
	1460	& ' and South segment: ',npckgs(ib2)
	1461	CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' )
[3651]	1462	ELSE
[10425]	1463	WRITE(ctmp1,*) ' E R R O R : Check South and East Open boundary indices'
	1464	WRITE(ctmp2,*) ' ========== Crossing problem with East segment: ',npckge(ib1), &
	1465	& ' and South segment: ',npckgs(ib2)
	1466	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1467	END IF
	1468	END IF
	1469	END DO
	1470	END DO
	1471	END IF
	1472	!
	1473	! North/West crossings
	1474	IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
	1475	DO ib1 = 1, nbdysegw
	1476	DO ib2 = 1, nbdysegn
	1477	IF (( jpindt(ib2)<=jpiwob(ib1) ).AND. &
	1478	& ( jpinft(ib2)>=jpiwob(ib1) ).AND. &
	1479	& ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
	1480	& ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
	1481	IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
	1482	! We have a possible North-West corner
	1483	! WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
	1484	! WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
	1485	icornw(ib1,2) = npckgn(ib2)
	1486	icornn(ib2,1) = npckgw(ib1)
	1487	ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
[10425]	1488	WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
[3651]	1489	& jpinft(ib2), jpjwdt(ib1)
[10425]	1490	WRITE(ctmp2,*) ' ========== Not allowed yet'
	1491	WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
[3651]	1492	& ' and North segment: ',npckgn(ib2)
[10425]	1493	CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' )
[3651]	1494	ELSE
[10425]	1495	WRITE(ctmp1,*) ' E R R O R : Check North and West Open boundary indices'
	1496	WRITE(ctmp2,*) ' ========== Crossing problem with West segment: ',npckgw(ib1), &
[3651]	1497	& ' and North segment: ',npckgn(ib2)
[10425]	1498	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1499	END IF
	1500	END IF
	1501	END DO
	1502	END DO
	1503	END IF
	1504	!
	1505	! North/East crossings
	1506	IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
	1507	DO ib1 = 1, nbdysege
	1508	DO ib2 = 1, nbdysegn
	1509	IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
	1510	& ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
	1511	& ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
	1512	& ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
	1513	IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
	1514	! We have a possible North-East corner
	1515	! WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
	1516	! WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
	1517	icorne(ib1,2) = npckgn(ib2)
	1518	icornn(ib2,2) = npckge(ib1)
	1519	ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
[10425]	1520	WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
[3651]	1521	& jpindt(ib2), jpjedt(ib1)
[10425]	1522	WRITE(ctmp2,*) ' ========== Not allowed yet'
	1523	WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
	1524	& ' and North segment: ',npckgn(ib2)
	1525	CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' )
[3651]	1526	ELSE
[10425]	1527	WRITE(ctmp1,*) ' E R R O R : Check North and East Open boundary indices'
	1528	WRITE(ctmp2,*) ' ========== Crossing problem with East segment: ',npckge(ib1), &
	1529	& ' and North segment: ',npckgn(ib2)
	1530	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1531	END IF
	1532	END IF
	1533	END DO
	1534	END DO
	1535	END IF
	1536	!
	1537	! 3. Check if segment extremities are on land
	1538	!--------------------------------------------
	1539	!
	1540	! West segments
	1541	DO ib = 1, nbdysegw
	1542	! get mask at boundary extremities:
	1543	ztestmask(1:2)=0.
	1544	DO ji = 1, jpi
	1545	DO jj = 1, jpj
	1546	IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &
	1547	& ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
	1548	IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &
	1549	& ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1)
	1550	END DO
	1551	END DO
[10425]	1552	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
[3651]	1553
	1554	IF (ztestmask(1)==1) THEN
	1555	IF (icornw(ib,1)==0) THEN
[10425]	1556	WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib)
	1557	WRITE(ctmp2,*) ' ========== does not start on land or on a corner'
	1558	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1559	ELSE
	1560	! This is a corner
[5656]	1561	IF(lwp) WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
[3651]	1562	CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
	1563	itest=itest+1
	1564	ENDIF
	1565	ENDIF
	1566	IF (ztestmask(2)==1) THEN
	1567	IF (icornw(ib,2)==0) THEN
[10425]	1568	WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib)
	1569	WRITE(ctmp2,*) ' ========== does not end on land or on a corner'
	1570	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1571	ELSE
	1572	! This is a corner
[5656]	1573	IF(lwp) WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
[3651]	1574	CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
	1575	itest=itest+1
	1576	ENDIF
	1577	ENDIF
	1578	END DO
	1579	!
	1580	! East segments
	1581	DO ib = 1, nbdysege
	1582	! get mask at boundary extremities:
	1583	ztestmask(1:2)=0.
	1584	DO ji = 1, jpi
	1585	DO jj = 1, jpj
	1586	IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &
	1587	& ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
	1588	IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &
	1589	& ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1)
	1590	END DO
	1591	END DO
[10425]	1592	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
[3651]	1593
	1594	IF (ztestmask(1)==1) THEN
	1595	IF (icorne(ib,1)==0) THEN
[10425]	1596	WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib)
	1597	WRITE(ctmp2,*) ' ========== does not start on land or on a corner'
	1598	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1599	ELSE
	1600	! This is a corner
[5656]	1601	IF(lwp) WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
[3651]	1602	CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
	1603	itest=itest+1
	1604	ENDIF
	1605	ENDIF
	1606	IF (ztestmask(2)==1) THEN
	1607	IF (icorne(ib,2)==0) THEN
[10425]	1608	WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib)
	1609	WRITE(ctmp2,*) ' ========== does not end on land or on a corner'
	1610	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1611	ELSE
	1612	! This is a corner
[5656]	1613	IF(lwp) WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
[3651]	1614	CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
	1615	itest=itest+1
	1616	ENDIF
	1617	ENDIF
	1618	END DO
	1619	!
	1620	! South segments
	1621	DO ib = 1, nbdysegs
	1622	! get mask at boundary extremities:
	1623	ztestmask(1:2)=0.
	1624	DO ji = 1, jpi
	1625	DO jj = 1, jpj
	1626	IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &
	1627	& ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
	1628	IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &
	1629	& ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1)
	1630	END DO
	1631	END DO
[10425]	1632	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
[3651]	1633
	1634	IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
[10425]	1635	WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib)
	1636	WRITE(ctmp2,*) ' ========== does not start on land or on a corner'
	1637	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1638	ENDIF
	1639	IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
[10425]	1640	WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib)
	1641	WRITE(ctmp2,*) ' ========== does not end on land or on a corner'
	1642	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1643	ENDIF
	1644	END DO
	1645	!
	1646	! North segments
	1647	DO ib = 1, nbdysegn
	1648	! get mask at boundary extremities:
	1649	ztestmask(1:2)=0.
	1650	DO ji = 1, jpi
	1651	DO jj = 1, jpj
	1652	IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &
	1653	& ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
	1654	IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &
	1655	& ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1)
	1656	END DO
	1657	END DO
[10425]	1658	CALL mpp_sum( 'bdyini', ztestmask, 2 ) ! sum over the global domain
[3651]	1659
	1660	IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
[10425]	1661	WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib)
	1662	WRITE(ctmp2,*) ' ========== does not start on land'
	1663	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1664	ENDIF
	1665	IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
[10425]	1666	WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib)
	1667	WRITE(ctmp2,*) ' ========== does not end on land'
	1668	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1669	ENDIF
	1670	END DO
	1671	!
	1672	IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
	1673	!
	1674	! Other tests TBD:
	1675	! segments completly on land
	1676	! optimized open boundary array length according to landmask
	1677	! Nudging layers that overlap with interior domain
	1678	!
	1679	END SUBROUTINE bdy_ctl_seg
	1680
[9019]	1681
[3651]	1682	SUBROUTINE bdy_ctl_corn( ib1, ib2 )
	1683	!!----------------------------------------------------------------------
	1684	!! * ROUTINE bdy_ctl_corn *
	1685	!!
	1686	!! ** Purpose : Check numerical schemes consistency between
	1687	!! segments having a common corner
	1688	!!
	1689	!! ** Method :
	1690	!!----------------------------------------------------------------------
	1691	INTEGER, INTENT(in) :: ib1, ib2
	1692	INTEGER :: itest
	1693	!!----------------------------------------------------------------------
	1694	itest = 0
	1695
[6140]	1696	IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) ) itest = itest + 1
	1697	IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) ) itest = itest + 1
	1698	IF( cn_tra (ib1) /= cn_tra (ib2) ) itest = itest + 1
[3651]	1699	!
[6140]	1700	IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) ) itest = itest + 1
	1701	IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) ) itest = itest + 1
	1702	IF( nn_tra_dta (ib1) /= nn_tra_dta (ib2) ) itest = itest + 1
[3651]	1703	!
[6140]	1704	IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) ) itest = itest + 1
[3651]	1705	!
[6140]	1706	IF( itest>0 ) THEN
[10425]	1707	WRITE(ctmp1,*) ' E R R O R : Segments ', ib1, 'and ', ib2
	1708	WRITE(ctmp2,*) ' ========== have different open bdy schemes'
	1709	CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
[3651]	1710	ENDIF
	1711	!
	1712	END SUBROUTINE bdy_ctl_corn
	1713
[911]	1714	!!=================================================================================
	1715	END MODULE bdyini

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

Download in other formats: