FUNCTION YINT(XL,YL,XL0)
C
C     Quadratic interpolation routine
C
C     Input:  XL - array of x
C             YL - array of f(x)
C             XL0 - the point x(0) to which one interpolates
C
      INCLUDE 'PARAMS.FOR'
      DIMENSION XL(3),YL(3)
      A0=(XL(2)-XL(1))*(XL(3)-XL(2))*(XL(3)-XL(1))
      A1=(XL0-XL(2))*(XL0-XL(3))*(XL(3)-XL(2))
      A2=(XL0-XL(1))*(XL(3)-XL0)*(XL(3)-XL(1))
      A3=(XL0-XL(1))*(XL0-XL(2))*(XL(2)-XL(1))
      YINT=(YL(1)*A1+YL(2)*A2+YL(3)*A3)/A0
      RETURN
      END