c  SPLI2ROW.FOR
c  
c  Read an x,y,f(x) file, interpolates equispaced values, and writes
c  as a 1-row grid.  Note the the y value is a dummy, and the xyz file
c  is assumed to be 1D.
c
c       input: f(i) is the functional value at location x(i); i=1,n
c
	dimension f(1000),h(1000),x(1000),p(1000),id(14),pgm(2)
	character*50 ifile,ofile
      print *,' Enter ascii x,y,f(x) input file name:'
      print 100
100   format(' *',$)
      read (5,101) ifile
101   format(a50)
      open (7,file=ifile,status='old',form='formatted')
      open (11,form='formatted',status='scratch')
      print *,' Enter output 1-row grid file name:'
      print 100
      read (5,101) ofile
      print*,'Enter id:'
      print 100
      read(5,16)(id(i),i=1,14)
16    format(14a4)
      open (12,file=ofile,status='new',form='unformatted')
      print*,' Enter new dx:'
      print 100
      read*,dx
      n=1
10    read(7,*,end=20) x(n),dum,f(n)
      n=n+1
      go to 10
20    n=n-1
      x0=x(1)
      y0=dum
      dy=1
      nrow=1
      nz=1
      write(11,*) x(1),f(1)
      h(1)=0
      do 30 i=2,n
30      h(i)=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
	ncol=1
	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=x(i-1)+(j*dx)
	if(y.gt.x(i)) go to 50
	ncol=ncol+1
	g=y-x(i-1)
	q=c0+g*(c1+g*(c2+g*c3))
      write(11,*) y,q
	go to 40
50      continue
      pgm(1)='spli'
      pgm(2)='ne  '
      write(12)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      rewind 11
      do 51 i=1,ncol
51    read(11,*)dum,f(i)
      write(12)x0,(f(i),i=1,ncol)
	stop
	end
      subroutine splnon(f,h,p,n)
c
c  finds the horizontal derivative for irregularly-spaced data.
c
      dimension f(n),h(n),p(n),c(1000),d(1000)
      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
