	subroutine  interv(   xt   ,  lxt   ,  x    ,  left   , mflag  )
        implicit double precision (a-h,o-z)
c
c
c	computes
c		left  :=  max(i : 1.le.i.le.lxt .and. xt(i).le.x ) .
c
c	**** I N P U T ****
c
c	xt ...	a real sequence, of length lxt, assumed to be
c		nondecresing.
c	lxt...	number of terms in the sequnce xt.
c	x  ...	the point whose location with respect to the sequence xt
c		is to be determined.
c
c      **** O U T P U T ****
c
c	left,mflag     both integers,whose value is
c	 1    -1       if  x .lt. xt(1)
c	 i     0       if  xt(i)   .le. x .lt. xt(i+1)
c	 lxt   1       if  xt(lxt) .le. x
c		In particular, mflag = 0 is the 'usual' case.
c		Mflag .ne. 0 indicates	thet x lies outside the halfopen
c		interval  xt(1) .le. y .lt. xt(lxt). the asymmetric
c		treatment of the interval is due to the decision to
c		make all  pp functions continuos from the right.
c
c     **** M E T H O D ****
c
c     The program is designed to be efficient in the common situation that
c     it is called repeatedly with taken from an increasing or decreasing
c     sequence. This will happen, e.g. ,when a pp function is to be
c     graphed.The first guess for left is therefore taken to be the value
c     returned at the previsiuos call and stored in the L O C A L variable
c     ilo  and are done after just three comparisons.
c	 Otherwise, we repeatedly double the difference istep=ihi-ilo
c     while also moving ilo and ihi in the direction of x , until
c		      xt(ilo).le. x .lt. xt(ihi) ,
c     after which we use bisection to get,in addition, ilo+1= ihi .
c     Left= ilo is then returned .
c
      integer	left  , lxt   , mflag  ,  ihi	,
     *		ilo   , istep  ,  middle
c
      double precision	 x     , xt    (  lxt	)
c
      data     ilo   / 1  /
c
c     save  ilo (a valid fortran statement in the new 1977 standard)
c
	if  ( x . ne . xt ( lxt ) )  go to  5
	left = lxt - 1
	mflag = 0
	return
 5	ihi  =	ilo + 1
	if  ( ihi .lt. lxt )  go  to  20
c
	  if  ( x . ge . xt ( lxt ) )  go  to  110
	    if	( lxt . le . 1 )  go to  90
	      ilo = lxt - 1
	      ihi  =  lxt
c
20			continue
			if  (x .ge. xt(ihi))  go to  40
				if  (x .ge. xt(ilo))  go to  100
c
c					now  x.lt.xt(ilo).
c					decrease  ilo  to capture  x.
c
					istep  =  1
31					continue
				    	    ihi  =  ilo
					    ilo  =  ihi - istep
					    if	(ilo .le. 1)
     *						go to  35
						if  (x.ge.xt(ilo))
     *						    go to  50
						    istep  =  istep*2
					go  to	31
35					continue
					ilo  =	1
					if  (x .lt. xt(1))  go to  90
					    go	to  50
c
40	  		continue
c
c			now  x .ge. xt(ihi).
c 			increase  ihi  to capture  x.
c
			istep  =  1
41			continue
				ilo  =	ihi
				ihi  =	ilo + istep
				if  (ihi .ge. lxt)  go to  45
					if  (x .lt. xt(ihi))  go to  50
					    istep  =  istep*2
			go  to	41
45			continue
			if  (x .ge. xt(lxt))  go to  110
 				ihi  =	lxt
c
c				now  xt(ilo) .le. x .lt. xt(ihi).
c				narrow the interval.
50			continue
c
				middle	=  (ilo + ihi)/2
				if  (middle .eq. ilo)  go to  100
c
c					note. it is assumed that
c					middle = ilo in case ihi = ilo+1
c
					if  (x .lt. xt(middle))
     *					    go	to  53
c
					    ilo  =  middle
					go to  50
53					continue
					ihi  =	middle
			go to  50
c
90		continue
c
c		set output flag and return
c
		mflag = -1
		left = 1
	return
100	continue
	mflag  =  0
	left   =  ilo
	return
c
110	continue
	mflag  =  1
	left   =  lxt
 	return
c
	end
 
 

