subroutine rkfliq(rt,y,yp)

      implicit double precision(a-h,o-z)

      common/misc/l,lhat,lindex,nsurf
      common/dmisc/period,grav,pi,pi4,p43,eps,verg,eig,eigt,y3i,y3t,
     1  amass
      common/splinq/vq(6)
      common/ray/iray
      real*8 l,lhat,lindex,y(8),yp(8),z(4)
*
* order x,g/r,g/r/sl**2,n**i*r/g,u,fac,r
*
      call splintl(rt,vq)
      yp(1)=y(1)*(vq(2)-3.+lindex)+y(2)*(lhat*vq(1)/eigt-vq(2))+
     1  y(3)*vq(2)
      yp(2)=y(1)*(eigt/vq(1)-vq(3))+y(2)*(1.-vq(4)+vq(3)+lindex)-y(3)*
     1 vq(3)
      yp(3)=y(3)*(1.-vq(4)+lindex)+y(4) 
      yp(4)=y(1)*vq(4)*vq(3)+y(2)*vq(4)*vq(2)+y(3)*(lhat-vq(4)*vq(2)) 
     1  +y(4)*(lindex-vq(4))
      do 1 i=1,4
        yp(i)=vq(5)*yp(i)
1     continue
      if (iray.eq.0) return
      rrt=1./vq(6)**lindex
      do 3 i=1,4
        z(i)=y(i)*rrt 
3     continue
      rhot=vq(4)*vq(1)/pi4/grav
      yp(5)=rhot*vq(6)**5*(z(1)**2+lhat*z(2)**2*(vq(1)/eigt)**2)
      rrt=rhot*vq(1)*vq(6)**5 
      yp(6)=rrt*vq(2)*(z(2)-z(3))**2
      yp(7)=rrt*vq(3)*z(1)**2 
      yp(8)=-rrt*(z(4)+(l+1.)*z(3))**2/vq(4)
      do 2 i=5,8
        yp(i)=yp(i)*vq(5)
2     continue

      return
      end

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