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_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/BDY – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/BDY/bdyini.F90 @ 11641

Last change on this file since 11641 was 11641, checked in by acc, 5 years ago
Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Further substantive changes to the BDY routines that read multiple copies of some namelists from namelist_cfg. Special treatment is required to replicate the original behaviour with internal files. These changes emulate previous behaviour but removal of the reliance on multiple, same-named namelist blocks would be preferable.
Property svn:keywords set to `Id`
File size: 89.6 KB

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

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

Download in other formats: