New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

bdyini.F90 in NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90 @ 11258

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

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

Download in other formats: