      subroutine findvel(y,velr,dvdxr,dvdzr,d2vdq2)
      implicit NONE
      real*8 y(*),velr,dvdxr,dvdzr,d2vdq2
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             dvdzr - 1st derivative with respect to global z.
c             d2vdq2 - 2nd derivative with respect to ray centred 
c                      coordinate.
c--------------------------------------------------------------            

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

      real*4 x1,z1,t,u
      real*4 yf4(4),yp1f4(4),yp2f4(4),yppf4(4),c(4,4)
      real*4 vel2,d2vdx2,d2vdz2,dvdxz

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

      x1 = y(1) - xminco
      z1 = y(2) - zminco
      
c     find nodes around the point.
      
      ncol = int(x1*rdelx) + 1
      nlin = int(z1*rdelz) + 1


c     make the  point is within the model.

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

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

         num1 = (ncol-1)*numz + nlin
         num4 = num1 + 1
         num2 = ncol*numz + nlin
         num3 = num2 + 1
         
         yf4(1) = velp(num1) 
         yf4(4) = velp(num4) 
         yf4(2) = velp(num2) 
         yf4(3) = velp(num3) 
         
         yp1f4(1) = dvdx(num1) 
         yp1f4(4) = dvdx(num4) 
         yp1f4(2) = dvdx(num2) 
         yp1f4(3) = dvdx(num3) 
         
         
         yp2f4(1) = dvdz(num1) 
         yp2f4(4) = dvdz(num4) 
         yp2f4(2) = dvdz(num2) 
         yp2f4(3) = dvdz(num3) 
         
         
         yppf4(1) = d2vdxz(num1) 
         yppf4(4) = d2vdxz(num4) 
         yppf4(2) = d2vdxz(num2) 
         yppf4(3) = d2vdxz(num3) 
         
         
c     compute coefficients c
         
         call bcucof(yf4,yp1f4,yp2f4,yppf4,delx,delz,c)
         
         t = (x1 - real(ncol-1)*delx)*rdelx
         u = (z1 - real(nlin-1)*delz)*rdelz
         
         velr=0.
         dvdxr=0.
         dvdzr=0.
         d2vdx2 = 0.0
         d2vdz2 = 0.0
         do i = 4,1,-1
            velr = t*velr + ((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
            dvdzr = t*dvdzr + (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)
            d2vdz2 = t*d2vdz2 + 6.*c(i,4)*u + 2.*c(i,3)
         enddo
         dvdxz = 0.0
         do i=4,2,-1
            dvdxz = t*dvdxz + real(i-1)*
     *                          (3.*c(i,4)*u + 2.*c(i,3))*u + c(i,2)
         enddo

         vel2   = velr*velr
         dvdxr  = dvdxr*rdelx
         dvdzr  = dvdzr*rdelz
         d2vdx2 = d2vdx2*rdelx*rdelx
         d2vdz2 = d2vdz2*rdelz*rdelz
         dvdxz  = dvdxz*rdelx*rdelz
         
         d2vdq2 = vel2*y(4)*y(4)*d2vdx2 - 2.*vel2*y(4)*y(3)*dvdxz
     *        + vel2*y(3)*y(3)*d2vdz2
C      endif
      return
      end



