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

Context Navigation

source: NEMO/trunk/src/OCE/BDY/bdyini.F90 @ 10983

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

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

Download in other formats: