      function velocitn(x1,tau1)
      implicit NONE
      real*4 x1,tau1
      real*4 velocitn
c--------------------------------------------------------------      
c     Purpose:
c              Find velocity at (x1,tau1)
c
c--------------------------------------------------------------      
c     Input :  
c              x1  x-coordinate
c              tau1  tau-coordinate
c-------------------------------------------------------------      

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

      real*4 yf4(4),yp1f4(4),yp2f4(4),yppf4(4),c(4,4),u,t,vel
      real*4 x1c,tau1c,eps/1.0e-6/

      integer ncol,nlin,num1,num2,num3,num4,i
      
c     correct coordinates relative to velocity model.

      x1c = x1 - xminco
      tau1c = tau1 - tauminco

C	write(0,*)'taulength=',taulength
C	write(0,*)'xlength=',xlength

C	write(0,*)'x1c=',x1c
C	write(0,*)'tau1c=',tau1c


c     check that point is within the model, if not: make them:

      x1c = amax1(amin1(x1c,xlength*(1.0-eps)),0.0)
      tau1c = amax1(amin1(tau1c,taulength*(1.0-eps)),0.0)

C	write(0,*)'after: x1c=',x1c
C	write(0,*)'after: tau1c=',tau1c

c     find nodes around the point.

      ncol = int(x1c*rdelx) + 1
      nlin = int(tau1c*rdeltau) + 1

c     check that point is within the model, if not: make it:

c      ncol = max0(min0(ncol,numx-1),1)
c      nlin = max0(min0(nlin,numtau-1),1)

      num1 = (ncol-1)*numtau + nlin
      num4 = num1 + 1
      num2 = ncol*numtau + 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) = dvdtau(num1) 
      yp2f4(4) = dvdtau(num4) 
      yp2f4(2) = dvdtau(num2) 
      yp2f4(3) = dvdtau(num3) 

      yppf4(1) = d2vdxtau(num1) 
      yppf4(4) = d2vdxtau(num4) 
      yppf4(2) = d2vdxtau(num2) 
      yppf4(3) = d2vdxtau(num3) 

c     compute coefficients 

      call bcucof(yf4,yp1f4,yp2f4,yppf4,delx,deltau,c)

      t =(x1c - real(ncol-1)*delx)*rdelx
      u =(tau1c - real(nlin-1)*deltau)*rdeltau 

      vel = 0.d0
      do i = 4,1,-1
        vel = t*vel+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
      enddo

      velocitn = vel
      return
      end








