      program cubic spline
c
c	now reading N and oversampling factor
c	from file new_cubic.in !!!!
c
      implicit none
      integer i,n,m,idum,novers,nhead,iprod,NMAX,lr,icount
      PARAMETER (NMAX=749580)
      PARAMETER(lr=16)
      real*8 yp1,ypn,xint,yint,dx,start,step
      real*8 x(NMAX),y(NMAX),y2(NMAX)
      character*40 fin,fout1,fout2,fout3
c
c      open(unit=02,file='makemaster.par',status='old')
c
c      read(2,*)
c      read(2,*)
c      read(2,*)
c      read(2,*)
c      read(2,1) fin
c 1    format(A40)
c
c      write(*,1) fin
c
c     read(2,*)
c      read(2,*) nhead
c
       write(*,*) '*** CUBIC SPLINE ****'

c      read(2,*)
c      read(2,*) n

c	n = 19500
c	novers = 10

	OPEN(42,file='new_cubic.in',type='old')
	READ(42,*) n
	READ(42,*) novers
	CLOSE(42)

	WRITE(*,*)'Number of points to interpolate: ',n
      	write(*,*)'Oversampling factor: ',novers
	
      if (n.gt.NMAX) then
         write(*,*) "#############################################"
         write(*,*) "### Error: number of data points > 749580 ###"
         write(*,*) "### Program aborted ..................... ###"
         write(*,*) "#############################################"
         goto 9999
      endif
c      read(2,*)
c      read(2,*) start
c      read(2,*)
c      read(2,*) step
c      write(*,*) 'step= ',step
c
      open(unit=01,file='spline.in',status='old')

      do i = 1,n,1

	read(1,*,err=1111) x(i),y(i)

c         read(1,*,err=1111) idum,y(i)
c         x(i) = start+dfloat(i-1)*step
c         write(*,*) i,x(i),y(i)

      enddo
	CLOSE(01)
c
c      if (iprod.eq.2) then
c         open(unit=07,file=fout1,status='new')
c         open(unit=08,file=fout2,status='new')
c         do i = 1,n,1
c            write(7,11) x(i),y(i)
c         enddo
c         close(07)
c      endif
   11 format(1x,f15.9,1x,f8.6) 
   12 format(1x,f15.9,1x,f11.9)
c
 1111 continue
c
      CALL SPLINE(x,y,n,1.d30,1.d30,y2)
c
c      open(unit=03,file=fout3,status='new',access='direct',
c     *     form='unformatted',recl=lr)

	OPEN(02,file='spline.out')
c
      dx = (x(n)-x(1))/dfloat(n-1)
      m = novers*(n-1)+1
      icount = 0
      do i = 1,m,1
         xint = x(1)+dfloat(i-1)/dfloat(novers)*dx
         CALL SPLINT(x,y,y2,n,xint,yint)
         if (yint.lt.1.0d-8) then
            icount = icount+1
            yint=1.0d-8
         endif

c         write(*,*) i,xint,yint
	WRITE(02,333) i,xint,yint

333	FORMAT(1X,I6,1X,F7.1,1X,F20.10)

c         write(unit=3,rec=i) xint,yint
c         if (iprod.eq.2) then
c            write(8,12) xint,yint
c         endif
      enddo
      close(02)
c      close(08)
c
c      write(*,*) '##########################################'
      write(*,*) '### Check for spline values < 1.0d-8 : ###'
      write(*,*) '### Values found = ',icount,'          ###'
      write(*,*) '### These values were set to 1.0d-8 .  ###'
c      write(*,*) '##########################################'
c
 9999 continue
c
      end
C==============================================================================
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
      INTEGER n,NMAX
      REAL*8 yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=749580)
      INTEGER i,k
      REAL*8 p,qn,sig,un,u(NMAX)
      if (yp1.gt.0.99d30) then
        y2(1)=0.d0
        u(1)=0.d0
      else
        y2(1)=-0.5d0
        u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.d0
        y2(i)=(sig-1.)/p
        u(i)=(6.d0*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
11    continue
      if (ypn.gt.0.99d30) then
        qn=0.d0
        un=0.d0
      else
        qn=0.5d0
        un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 0=M,173+9.
C==============================================================================
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      INTEGER n
      REAL*8 x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      REAL*8 a,b,h
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.) pause 'bad xa input in splint'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *2)/6.d0
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 0=M,173+9.
