C+
C NAME:
C	nrSpline
C PURPOSE:
C	Calculate second derivatives in X locations of interpolating spline
C	for Y fnc-values (the 2nd derivative array serves as input to
C	the interpolation subroutine nrSplInt)
C CATEGORY:
C	Math: interpolation and extrapolation
C CALLING SEQUENCE:
	subroutine nrSpline_old(N,X,Y,YP1,YPN,Y2)
C INPUTS:
C	N	integer		# points and fnc-values
C	X(N)	real		x-values; MUST be sorted
C	Y(N)	real		function values
C	YP1	real		derivative at point 1
C				if YP1>=1.e30 then 2nd derivative is zero in point 1
C	YPN	real		derivative at point N
C				if YP1>=1.e30 then 2nd derivative is zero in point N
C OUTPUTS:
C	Y2(N)	real		second derivatives at locations X
C CALLS:
C	Say
C RESTRICTIONS:
C	The X array must be sorted.
C	NO CHECK IS MADE TO ENSURE THAT THIS IS THE CASE
C PROCEDURE:
C	See Numerical Recipes, par. 3.4, p. 88
C MODIFICATION HISTORY:
C	JUN-1993, Paul Hick (UCSD)
C-
	    integer	N
	    real	X (N)
	    real	Y (N)
	    real	Y2(N)

	parameter      (NMAX = 1000)
	real		U(NMAX)

	if (N .gt. NMAX) call Say('nrSpline','E','Stop','Too small parameter NMAX')

	if (YP1 .gt. 0.99e30) then
	    Y2(1) = 0.0
	    U(1)  = 0.0

	print *, 'In nrSpline 1a Y2 U',Y2(1),U(1)

	else
	    Y2(1) = -0.5
	    U(1)  = (3.0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)

	    print *, 'In nrSpline 1b Y2 U',Y2(1),U(1)

	end if

	do I=2,N-1
	    SIG = (X(I)-X(I-1))/(X(I+1)-X(I-1))

	    P = SIG*Y2(I-1)+2.0
	    Y2(I) = (SIG-1.0)/P

	    denom1 = (X(I+1)-X(I))-(Y(I)-Y(I-1))
	    if(abs(denom1).lt.1.0e-5) denom1 = 1.0e-5

	    denom2 = X(I)-X(I-1)
	    if(abs(denom2).lt.1.0e-5) denom2 = 1.0e-5

	    denom3 = (X(I+1)-X(I-1))-SIG*U(I-1)
	    if(abs(denom3).lt.1.0e-5) denom3 = 1.0e-5

	    if(abs(P).lt.1.0e-5) P = 1.0e-5

	    U(I) = (6.0*((Y(I+1)-Y(I))/denom1
     &		/(denom2))/denom3)/P

C	    U(I) = (6.0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))
C     &		/(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P

	    print *, ' '
	    print *, 'In nrSpline 2 I SIG P      ',I,SIG,P

	    if(SIG.eq.1.0) then 
		Y2(I) = (1.00001-1.0)/P
		U(I)  = -1.0e5
	    end if



	    print *, 'X(I) X(I-1) X(I)-X(I-1)    ',X(I),X(I-1),X(I)-X(I-1)
	    print *, 'X(I+1) X(I-1) X(I+1)-X(I-1)',X(I+1),X(I-1),X(I+1)-X(I-1)
	    print *, 'Y(I+1) Y(I) Y(I+1)-Y(I)    ',Y(I+1),Y(I),Y(I+1)-Y(I)
	    Print *, '(X(I+1)-X(I))-(Y(I)-Y(I-1))',(X(I+1)-X(I))-(Y(I)-Y(I-1))
	    Print *, '6.0*                       ',6.0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))/(X(I)-X(I-1)))
	    Print *, '(X(I+1)-X(I-1))-SIG*U(I-1))',(X(I+1)-X(I-1))-SIG*U(I-1)
	    print *, 'Y2(I),U(I)                 ',Y2(I),U(I)


	end do

	if (YPN .gt. 0.99e30) then
	    QN = 0.0
	    UN = 0.0

	    print *, 'In nrSpline 3a QN UN',QN, UN

	else
	    QN = 0.5
	    UN = (3.0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))

	    print *, 'In nrSpline 3b QN UN',QN, UN

	end if

	print *, 'In nrSpline 4a N,UN,QN,U(N-1),Y2(N-1)',N,UN,QN,U(N-1),Y2(N-1)


	denom4 = (QN*Y2(N-1)+1.0)
	if(abs(denom4).lt.1.0e-5) denom4 = 1.0e-5

	Y2(N) = (UN-QN*U(N-1))/denom4
C	Y2(N) = (UN-QN*U(N-1))/(QN*Y2(N-1)+1.0)

	print *, 'In nrSpline 4b N Y2(N)',N,Y2(N)

	do K=N-1,1,-1
	    Y2(K) = Y2(K)*Y2(K+1)+U(K)

	    print *, 'In nrSpline 5 K Y2(I)',K,Y2(K)

	end do

	return
	end