subroutine splintl(xint,yout)
*
*     interpolates for liquid 
*
      implicit double precision(a-h,o-z)

      common g(650),x(650),rho(650),yliq(4,650),ray(4,650)
      common/savingl/k,kq
      common/setup/isetup
      common/mfac/m
      common/yfacs/y1(650),y2(650),y3(650),y4(650),y5(650)
      common/rs/r(650)
c       common/savel/m,c1(4,650),y1(650),c2(4,650),y2(650),
c     1  c3(4,650),y3(650),c4(4,650),y4(650),c5(4,650),
c     2  y5(650),c6(4,650)
      real*8 yout(6),y6(650)
      real*8 c1(4,650),c2(4,650),c3(4,650),c4(4,650),c5(4,650)
      real*8 c6(4,650)
*
* compute spline coefficients for a model if we haven't done so
*
      if (isetup.eq.1)then
        do 180 i=1,m
          y6(i)=r(i)
180     continue
        call consl(x,y1,m,c1)
        call consl(x,y2,m,c2)
        call consl(x,y3,m,c3)
        call consl(x,y4,m,c4)
        call consl(x,y5,m,c5)
        call consl(x,y6,m,c6)
        isetup=0
      endif
      mm=m-1
      if (xint.ge.x(1).and.xint.lt.x(m)) go to 40 
      if (xint .lt. (1.+1.e-8)*x(1)) go to 110
      if (xint .ge. x(m)) go to 20
      k=1 
      go to 70
20    if (xint .gt. (1.+1.e-6)*x(m)) go to 110
      k=mm
      go to 70
40    il=1
      ir=m
50    k=il+((ir-il)/2)
      if (xint.ge.x(k)) go to 60
      ir=k
      go to 50
60    if (xint.lt.x(k+1)) go to 70
      il=k
      go to 50
70    continue
      x1=x(k+1)-xint
      xx=xint-x(k)
      x12=x1*x1
      xx2=xx*xx
      yout(1)=x1*(c1(1,k)*x12 +c1(3,k))+xx*(c1(2,k)*xx2 +c1(4,k))
      yout(2)=x1*(c2(1,k)*x12 +c2(3,k))+xx*(c2(2,k)*xx2 +c2(4,k))
      yout(3)=x1*(c3(1,k)*x12 +c3(3,k))+xx*(c3(2,k)*xx2 +c3(4,k))
      yout(4)=x1*(c4(1,k)*x12 +c4(3,k))+xx*(c4(2,k)*xx2 +c4(4,k))
      yout(5)=x1*(c5(1,k)*x12 +c5(3,k))+xx*(c5(2,k)*xx2 +c5(4,k))
      yout(6)=x1*(c6(1,k)*x12 +c6(3,k))+xx*(c6(2,k)*xx2 +c6(4,k))
      return
110   continue
      return
700   continue

      stop 
      end

************************************************************************