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 branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 4333

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

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

Download in other formats: