C     ******************************************************************

	FUNCTION YPARINTERP (XINT,X,Y,N,NTOT,lup)
C     ==========================
C
C     M. Cervino (Obs. Midi-Pyrenees, Toulouse, France 1999) 
C
C     Parabolic interpolation routine. Determines YINT = Y(XINT) from 
C     grid Y(X) with N points and dimension NTOT.
C     As far as 3 point are necessary, assuming that X(I) < X < X(I+1)
C       (or X(I) > X > X(I+1) ) the interpolation can be 
C       performed ussing:
C               a) Y(I), Y(I+1) and Y(I+2)  : lup=.TRUE.
C               b) Y(I-1), Y(I) and Y(I+1)  : lup=.FALSE.
C     Linear interpolations are always assumed in the top (low) boundaries 
C        of the table if lup=.TRUE. (.FALSE.)
C
C     CAUTIONS: a) Note the REAL*8 definitions
C               b) The use of lup=.true. or lup=.false.
C            may leads different results, so it is convenient
C            try to use a "physical" argument for the choice
C            example: assume x(i) < x(i+1) with x corresponds 
C                     to the initial mass, and y(i) to the life-time.
C              "physical" argument: lifetime of WR stars and normal 
C                     stars are, may be, different
C               so: it may be better  use the life-time of more massive stars
C                  for WR ones, i.e. lup=.TRUE. if X > Mwrlim
C                  and the life-time for low mases ones (lup=.false.) if
C                  X < Mwrlim, or use always lup=.true. better than always 
C                  lup=.false. 
C
C     Any question or problem related with this subroutine, 
C      please, contact me
C                            mcs@laeff.esa.es
C                            http://www.laeff.esa.es/users/mcs
C
C      
        IMPLICIT REAL*8    (A-H,O-Z)
        real*8 YPARINTERP
	DIMENSION X(NTOT),Y(NTOT)
	logical lup

C     bisection (see Numerical Recipes par 3.4 page 90)
c     obtain k where x is in the interval x(n) and x(n+1)

        if(xint.eq.x(1)) then
          yparinterp=y(1)
          return
        endif
        if(xint.eq.x(n)) then
          yparinterp=y(n)
          return
        endif
	JL = 0
	JU = N+1
 10	IF (JU-JL.GT.1) THEN
	  JM = (JU+JL)/2
	  IF ((X(N).GT.X(1)).EQV.(XINT.GT.X(JM))) THEN
            JL = JM
	  ELSE
            JU = JM
	  END IF
	  GO TO 10
	END IF
	J = JL
	IF (J.EQ.0.OR.J.EQ.N) CALL ERRPRI 
	1    ('YPARINTERP: EXTRAPOLATION')
C
	k=j
	if(j.eq.1) then
	  if(x(j).eq.x(j+1)) then
	    YPARINTERP = y(j)
	    return
	  else
	    if(x(j+1).eq.x(j+2)) then
	      YPARINTERP = y(j)+((y(j+1)-y(j))*(XINT-x(j))/(x(j+1)-x(j)))
	      return
	    else
	      k=2
	      lup=.false.
	    endif
	  endif
	endif
	
	if(j.eq.n-1) then
	  if(x(j).eq.x(j+1)) then
	    YPARINTERP = y(j)
	    return
	  else
	    if(x(j-1).eq.x(j+1)) then
	      YPARINTERP = y(j)+((y(j+1)-y(j))*(XINT-x(j))/(x(j+1)-x(j)))
	      return
	    else
	      k=j
	      lup=.false.
	    endif
	  endif
	endif

	if(x(k).eq.x(k+1)) then
	  YPARINTERP = y(k)
	  return
	else
	  if(lup) then
	    if(x(k+1).eq.x(k+2).or.x(k).eq.x(k+2)) then
	      YPARINTERP = y(k)+((y(k+1)-y(k))*(XINT-x(k))/(x(k+1)-x(k)))
	      return
	    else
	      help1=(y(k)-y(k+1))/((x(k)-x(k+1))*(x(k)-x(k+2)))
	      help2=(y(k+1)-y(k+2))/((x(k+1)-x(k+2))*(x(k)-x(k+2)))
	      A=help1-help2
	      B=(help2*(x(k)+x(k+1)))-(help1*(x(k+1)+x(k+2)))
	      C=y(k+1)-(A*(x(k+1)**2.d+00))-(b*x(k+1))
	      YPARINTERP = (A*(XINT**2.d+00))+(B*XINT)+C

	      ylintest = y(k)+((y(k+1)-y(k))*(XINT-x(k))/(x(k+1)-x(k)))
              ymax=max(y(k),y(k+1),y(k+2))
              ymin=min(y(k),y(k+1),y(k+2))
              if(yparinterp.gt.ymax.or.yparinterp.lt.ymin) 
     $             YPARINTERP=ylintest
	      return
	    endif
	  else
	    if(x(k-1).eq.x(k).or.x(k-1).eq.x(k+1)) then
	      YPARINTERP = y(k)+((y(k+1)-y(k))*(XINT-x(k))/(x(k+1)-x(k)))
	      return
	    else
	      help1=(y(k-1)-y(k))/((x(k-1)-x(k))*(x(k-1)-x(k+1)))
	      help2=(y(k)-y(k+1))/((x(k)-x(k+1))*(x(k-1)-x(k+1)))
	      A=help1-help2
	      B=(help2*(x(k-1)+x(k)))-(help1*(x(k)+x(k+1)))
	      C=y(k)-(A*(x(k)**2.d+00))-(b*x(k))
	      YPARINTERP = (A*(XINT**2.d+00))+(B*XINT)+C
	      ylintest = y(k)+((y(k+1)-y(k))*(XINT-x(k))/(x(k+1)-x(k)))
              ymax=max(y(k-1),y(k),y(k+1))
              ymin=min(y(k-1),y(k),y(k+1))
              if(yparinterp.gt.ymax.or.yparinterp.lt.ymin) 
     $             YPARINTERP=ylintest
	      return
	    endif
	  endif
	endif
C
	END

C
      SUBROUTINE ERRPRI (NAME)
C     =================
C   
C     From D. Schaerer
C     Prints name of subroutine where error occured.
C
      CHARACTER NAME*(*)
      WRITE (98,*) NAME
      STOP 1
      END
C
