      dimension x(2048),p(2048),d(2048),h(2048),pp(2048),f(2048)
      character id*56, pgm*8, ifile*50
      data ddval/1.e20/
      print *,'enter ifile:'
      read(5,100) ifile
  100 format(a)
      open(10,file=ifile,form='unformatted',status='old')
      read(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
      if(nc.gt.2048) stop 'grid dimension exceeds 2048'
      print *,'enter ofile:'
      read(5,100) ifile
      open(11,file=ifile,form='unformatted',status='unknown')
      print *,'enter id:'
      read(5,100) id
      pgm='sp_x    '
      write(11) id,pgm,nc,nr,nz,xo,dx,yo,dy
      do 90 j=1,nr
      read(10) dlx,(x(i),i=1,nc)
      do 10 i=1,nc
      if(x(i).ge.ddval) go to 20
   10 continue
      call spline(x,p,nc,dx)
      go to 90
   20 continue
      nn=0
      do 30 i=1,nc
      p(i)=x(i)
      if(x(i).ge.ddval) go to 30
      nn=nn+1
      f(nn)=x(i)
      d(nn)=float(i)
   30 continue
      if(nn.eq.0) go to 90
      h(1)=0.
      do 40 i=2,nn
      h(i)=d(i)-d(i-1)
   40 continue
c     find the derivatives pp(i) of f(i) at the locations d(i); i=1,nn
      call splnon(f,h,pp,nn)
      nn=0
      do 50 i=1,nc
      if(x(i).ge.ddval) then
        p(i)=x(i)
      else
        nn=nn+1
        p(i)=pp(nn)
      endif
   50 continue
   90 write(11) dlx,(p(i),i=1,nc)
      stop
      end
      subroutine spline(f,p,n,dx)
      dimension f(n),p(n),c(2048),d(2048)
      c(1)=.5
      d(1)=1.5*(f(2)-f(1))/dx
      n1=n-1
      do 10 i=2,n1
      dnm=(4.0-c(i-1))*dx
c      print *,'dnm = ',dnm,' d(i-1) = ',d(i-1)
      c(i)=dx/dnm
   10 d(i)=(3.0*(f(i+1)-f(i-1)) - dx*d(i-1))/dnm
c      print *,'c(n1) = ',c(n1)
      p(n)=(3.*(f(n)-f(n1))/dx - 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
      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
