subroutine inter(x,y,z,f,irow,ndim,eps)
c interpolation subroutine. x=input x. y=output y. z,f are table of valu
c  (z(x),f(y)).  irow=no. of pairs in table.  ndim=no. of pts. in interp
c  eps=tolerance in y (y-eps<y<y+eps).

      real z(irow),f(irow),arg(9000),val(9000),work(9000)
      if(irow)9,9,1
    9 write(6,8)
      go to 19
    1 n=ndim
      if(n-irow) 3,3,2
    2 n=irow
      ndim=irow
    3 b=0.
      do 5 i=1,irow
      delta=abs(z(i)-x)
      if(delta-b) 5,5,4
    4 b=delta
    5 work(i)=delta
      b=b+1.0
      do 10 j=1,n
      delta=b
      do 7 i=1,irow
      if(work(i)-delta) 6,7,7
    6 ii=i
      delta=work(i)
    7 continue
      arg(j)=z(ii)
      val(j)=f(ii)
   10 work(ii)=b
      if(ndim-1) 19,17,11
   11 delt2=0.
      do 16 j=2,ndim
      delt1=delt2
      iend=j-1
      do 12 i=1,iend
      h=arg(i)-arg(j)
      if(h) 12,23,12
   12 val(j)=(val(i)*(x-arg(j))-val(j)*(x-arg(i)))/h
      delt2=abs(val(j)-val(iend))
      if(j-2) 16,16,13
   13 if(delt2-eps) 20,20,14
   14 if(j-5) 16,15,15
   15 if(delt2-delt1) 16,21,21
   16 continue
   17 j=ndim
   18 y=val(j)
   19 return
   20 continue
      go to 18
   21 write (6,26)
   22 j=iend
      go to 18
   23 write (6,27)
      go to 22
    8 format(23h int error irow invalid)
   26 format(33h rounding errors cause inaccuracy)
   27 format (30h two identical argument values)
      end