c	spline interpolation
c
c       input: f(i) is the functional value and z(i) is the elevation
c       at location x(i); i=1,n
c
        parameter (MAXLEN = 2001)
        dimension f(MAXLEN),h(MAXLEN),x(MAXLEN),p(MAXLEN),z(MAXLEN)
        dimension t(MAXLEN)
	character*50 ifile,ofile
        print 16,'Enter input ascii x,z,f file name:'
16	format(1x,a,$)
	read (5,17) ifile
17	format(a50)
        open (7,file=ifile,status='old',form='formatted')
        print 16,'Enter output ascii x,z,f file name:'
	read (5,17) ofile
        open (11,file=ofile,status='unknown',form='formatted')
	h(1)=0
        read(7,*,end=20) x(1),z(1),f(1)
        n=1
10      read(7,*,end=20) x(n+1),z(n+1),f(n+1)
        if (x(n+1).ne.x(n)) n=n+1
        if(n.eq.MAXLEN) stop 'Too many points in profile'
	go to 10
20	continue
        write(11,*) x(1),z(1),f(1)
        dx=(x(n)-x(1))/float(n-1)
        sign=1.0
        if(dx.lt.0.0) sign=-1.0
	do 30 i=2,n
30      h(i)=sign*(x(i)-x(i-1))
c
c	find the derivatives p(i) of f(i) at the locations x(i); i=1,n
c
	call splnon(f,h,p,n)
c
c	perform the interpolation and write the results
c
	do 50 i=2,n
	c0=f(i-1)
	c1=p(i-1)
	c2=(3.0*(f(i)-c0)-h(i)*(p(i)+2.0*c1))/h(i)**2
	c3=(2.0*(c0-f(i))+h(i)*(p(i)+c1))/h(i)**3
        j=0
40	j=j+1
        y=sign*(x(1)+float(j)*dx)
        if(y.le.sign*x(i-1)) go to 40
        if(y.gt.sign*x(i)) go to 50
        g=y-sign*x(i-1)
        t(j)=c0+g*(c1+g*(c2+g*c3))
c        print *,i,j,t(j)
	go to 40
50	continue
        t(n)=f(n)
        do 60 i=2,n
60      h(i)=sign*(x(i)-x(i-1))
c
c       find the derivatives p(i) of z(i) at the locations x(i); i=1,n
c
        call splnon(z,h,p,n)
c
c	perform the interpolation and write the results
c
        do 80 i=2,n
        c0=z(i-1)
	c1=p(i-1)
        c2=(3.0*(z(i)-c0)-h(i)*(p(i)+2.0*c1))/h(i)**2
        c3=(2.0*(c0-z(i))+h(i)*(p(i)+c1))/h(i)**3
	j=0
70      j=j+1
        y=sign*(x(1)+float(j)*dx)
        if(y.le.sign*x(i-1)) go to 70
        if(y.gt.sign*x(i)) go to 80
        g=y-sign*x(i-1)
        q=c0+g*(c1+g*(c2+g*c3))
        write(11,*) sign*y,q,t(j)
        go to 70
80      continue
        stop
	end
      subroutine splnon(f,h,p,n)
c
c  finds the horizontal derivative for irregularly-spaced data.
c
      parameter (MAXLEN = 2000)
      dimension f(n),h(n),p(n),c(MAXLEN),d(MAXLEN)
      c(1)=.5
      d(1)=1.5*(f(2)-f(1))/h(2)
      n1=n-1
      do 10 i=2,n1
      dnm=2.*(h(i)+h(i+1))-h(i+1)*c(i-1)
      c(i)=h(i)/dnm
10    d(i)=(3.*(h(i)*(f(i+1)-f(i))/h(i+1)+h(i+1)*(f(i)-f(i-1))/h(i))-
     1     h(i+1)*d(i-1))/dnm
      p(n)=(3.*(f(n)-f(n1))/h(n)-d(n1))/(2.-c(n1))
      do 9 i=1,n1
      k=n-i
9     p(k)=d(k)-c(k)*p(k+1)
      return
      end
