      SUBROUTINE QINTRP(XM,X0,XP,FM,F0,FP,X,F)
      implicit real(selected_real_kind(15)) (a-h,o-z)
      D1=(FP-F0)/(XP-X0)
      D2=(FM-F0)/(XM-X0)
      B=(D1-D2)/(XP-XM)
      A=D1-B*(XP-X0)
      DEL=(X-X0)
      F=F0+DEL*(A+DEL*B)
      RETURN
      END
