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

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

bdyini.F90 in NEMO/trunk/src/OCE/BDY – NEMO

Context Navigation

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

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

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

Download in other formats: