subroutine grind(b1,b2)

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

      common g(650),x(650),rho(650),yliq(4,650),ray(4,650)
      common/misc/l,lhat,lindex,nsurf
      common/dmisc/period,grav,pi,pi4,p43,eps,verg,eig,eigt,y3i,y3t,
     1  amass
      common/ray/iray
      common/rs/r(650)
      common/splinq/vq(6)
      real*8 y(8),work(51),l,lhat,lindex
      integer iwork(5)
      external rkfliq

      yliq(1,1)=1.0 
      yliq(2,1)=eigt/(l*grav*p43*rho(1))
      y(1)=yliq(1,1)
      y(2)=yliq(2,1)
      yliq(3,1)=y3t 
      yliq(4,1)=l*yliq(3,1)
      y(3)=yliq(3,1)
      y(4)=yliq(4,1)
      iflag=1
      neqn=4
      if(iray.eq.0)go to 4
        do 6 j=1,4
          y(j+4)=0.
          ray(j,1)=0.
6       continue
        neqn=8
4     nsurf1=nsurf-1
      do 1 j=1,nsurf1
        k=j+1
        xstart=x(j)
        xfin=x(k)
        relerr=1.d-8
        abserr=dabs(y(1))
        do 2 i=2,neqn 
          t=dabs(y(i))
          abserr=dmin1(abserr,t)
2       continue
        abserr=1.d-8*abserr
        abserr=dmax1(abserr,1.d-20)
      call rkf(rkfliq,neqn,y,xstart,xfin,relerr,abserr,iflag,work,iwork)
        iflag=1
        do 3 i=1,4
          yliq(i,k)=y(i)
3       continue
        if (iray.eq.0) go to 1
        do 5 i=1,4
          ray(i,k)=y(i+4)
5       continue
1     continue
      rt=x(nsurf)
      call splintl(rt,vq)
      temp=1./vq(5)-1.
      b1=yliq(1,nsurf)*vq(4)+yliq(3,nsurf)*(l+1.)+yliq(4,nsurf)
*
* these are saio and cox bc's, changed
* from osaki and hansen's
*
      b2=yliq(1,nsurf)*((4.+r(nsurf)*eigt/g(nsurf))/temp-1.)
     1 +(1.-lhat*g(nsurf)/r(nsurf)/eigt/temp)*yliq(2,nsurf) 
     2 + ((l+1.)/temp-1.)*yliq(3,nsurf) 

      return
      end

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