      subroutine findvel(y,velr,dvdxr,dvdtaur,d2vdq2,
     *                   sigr,dsigdxr,dsigdtaur)
      implicit NONE
      real*8 y(*),velr,dvdxr,dvdtaur,d2vdq2,sigr,dsigdxr,dsigdtaur
c--------------------------------------------------------------      
c     Purpose:
c              Find velocity and some of the derivatives from the 
c              model.
c--------------------------------------------------------------      
c     Input :  
c             y - last updated solution vector.
c--------------------------------------------------------------      
c     Output : 
c             velr - velocity function
c             dvdxr - 1st derivative with respect to global x
c             dvdtaur - 1st derivative with respect to global tau.
c             d2vdq2 - 2nd derivative with respect to ray centred 
c                      coordinate.
c--------------------------------------------------------------            

      include '../Src/raypar.inc'
      include '../Src/model.inc'

      real*4 x1,tau1,t,u
      real*4 yf4(4),yp1f4(4),yp2f4(4),yppf4(4),c(4,4)
      real*4 ys4(4),yp1s4(4),yp2s4(4),ypps4(4),cs(4,4)
      real*4 vel2,d2vdx2,d2vdtau2,dvdxtau

      integer ncol,nlin,num1,num2,num3,num4,i

      x1 = y(1) - xminco
      tau1 = y(2) - tauminco
      
c     find nodes around the point.
      
      ncol = int(x1*rdelx) + 1
      nlin = int(tau1*rdeltau) + 1

        ncol=min(numx-1,max(2,ncol))
        nlin=min(numtau-1,max(2,nlin))

c     check that point is within the model.
C
C      if (ncol.lt.1.or.ncol.ge.numx.or.nlin.lt.1.or.nlin.ge.numtau) then
C         velr  = 0.d0
C         dvdxr = 0.d0
C         dvdtaur = 0.d0
C         d2vdq2= 0.d0
C
C         sigr = 0.d0
C         dsigdxr = 0.d0
C         dsigdtaur = 0.d0
C      else

         num1 = (ncol-1)*numtau + nlin
         num4 = num1 + 1
         num2 = ncol*numtau + nlin
         num3 = num2 + 1
         
         yf4(1) = velp(num1) 
         yf4(2) = velp(num2) 
         yf4(3) = velp(num3) 
         yf4(4) = velp(num4) 
         
         yp1f4(1) = dvdx(num1) 
         yp1f4(2) = dvdx(num2) 
         yp1f4(3) = dvdx(num3) 
         yp1f4(4) = dvdx(num4) 
         
         yp2f4(1) = dvdtau(num1) 
         yp2f4(2) = dvdtau(num2) 
         yp2f4(3) = dvdtau(num3) 
         yp2f4(4) = dvdtau(num4) 
         
         
         yppf4(1) = d2vdxtau(num1) 
         yppf4(2) = d2vdxtau(num2) 
         yppf4(3) = d2vdxtau(num3) 
         yppf4(4) = d2vdxtau(num4) 

         ys4(1) = sigmap(num1) 
         ys4(2) = sigmap(num2) 
         ys4(3) = sigmap(num3) 
         ys4(4) = sigmap(num4) 
         
         yp1s4(1) = dsigdx(num1) 
         yp1s4(2) = dsigdx(num2) 
         yp1s4(3) = dsigdx(num3) 
         yp1s4(4) = dsigdx(num4) 
         
         yp2s4(1) = dsigdtau(num1) 
         yp2s4(2) = dsigdtau(num2) 
         yp2s4(3) = dsigdtau(num3) 
         yp2s4(4) = dsigdtau(num4) 

         ypps4(1) = d2sigdxtau(num1) 
         ypps4(2) = d2sigdxtau(num2) 
         ypps4(3) = d2sigdxtau(num3) 
         ypps4(4) = d2sigdxtau(num4) 

         
c     compute coefficients c
         
         call bcucof(yf4,yp1f4,yp2f4,yppf4,delx,deltau,c)
         call bcucof(ys4,yp1s4,yp2s4,ypps4,delx,deltau,cs)
         
         t = (x1 - real(ncol-1)*delx)*rdelx
         u = (tau1 - real(nlin-1)*deltau)*rdeltau
         
         velr=0.d0
         dvdxr=0.d0
         dvdtaur=0.d0
         d2vdx2 = 0.d0
         d2vdtau2 = 0.d0

         sigr=0.d0
         dsigdxr=0.d0
         dsigdtaur=0.d0

         do i = 4,1,-1
            velr= t*velr + ((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
            dvdtaur= t*dvdtaur + (3.*c(i,4)*u+2.*c(i,3))*u + c(i,2)
            dvdxr= u*dvdxr + (3.*c(4,i)*t+2.*c(3,i))*t + c(2,i)
            d2vdx2= u*d2vdx2 + 6.*c(4,i)*t + 2.*c(3,i)
            d2vdtau2= t*d2vdtau2 + 6.*c(i,4)*u + 2.*c(i,3)

            sigr= t*sigr + ((cs(i,4)*u+cs(i,3))*u+cs(i,2))*u+cs(i,1)
            dsigdtaur= t*dsigdtaur +(3.*cs(i,4)*u+2.*cs(i,3))*u+ cs(i,2)
            dsigdxr= u*dsigdxr + (3.*cs(4,i)*t+2.*cs(3,i))*t + cs(2,i)

         enddo
         dvdxtau = 0.0
         do i=4,2,-1
            dvdxtau = t*dvdxtau + real(i-1)*
     *                          (3.*c(i,4)*u + 2.*c(i,3))*u + c(i,2)
         enddo

         vel2   = velr*velr

         dvdxr  = dvdxr*rdelx
         dvdtaur  = dvdtaur*rdeltau
         dsigdxr  = dsigdxr*rdelx
         dsigdtaur  = dsigdtaur*rdeltau

         d2vdx2 = d2vdx2*rdelx*rdelx

         d2vdtau2 = d2vdtau2*rdeltau*rdeltau
         dvdxtau  = dvdxtau*rdelx*rdeltau
         
         d2vdq2 = vel2*y(4)*y(4)*d2vdx2 - 2.*vel2*y(4)*y(3)*dvdxtau
     *        + vel2*y(3)*y(3)*d2vdtau2
C      endif

      return
      end


