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
	else
	    Y2(1) = -0.5
	    U(1)  = (3.0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
	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))
C	    if(abs(denom1).lt.1.0e-30) denom1 = 1.0e-30
	    if(denom1.lt.1.0e-30.and.denom1.gt.0.0) denom1 = 1.0e-30
	    if(denom1.gt.-1.0e-30.and.denom1.le.0.0) denom1 = -1.0e-30

	    denom2 = X(I)-X(I-1)
C	    if(abs(denom2).lt.1.0e-30) denom2 = 1.0e-30
	    if(denom2.lt.1.0e-30.and.denom2.gt.0.0) denom2 = 1.0e-30
	    if(denom2.gt.-1.0e-30.and.denom2.le.0.0) denom2 = -1.0e-30

	    denom3 = (X(I+1)-X(I-1))-SIG*U(I-1)
C	    if(abs(denom3).lt.1.0e-30) denom3 = 1.0e-30
	    if(denom3.lt.1.0e-30.and.denom3.gt.0.0) denom3 = 1.0e-30
	    if(denom3.gt.-1.0e-30.and.denom3.le.0.0) denom3 = -1.0e-30

C	    if(abs(P).lt.1.0e-30) P = 1.0e-30
	    if(P.lt.1.0e-30.and.P.gt.0.0) P = 1.0e-30
	    if(P.gt.-1.0e-30.and.P.le.0.0) P = -1.0e-30

	    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
	end do

	if (YPN .gt. 0.99e30) then
	    QN = 0.0
	    UN = 0.0
	else
	    QN = 0.5
	    UN = (3.0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
	end if

	denom4 = (QN*Y2(N-1)+1.0)
C	if(abs(denom4).lt.1.0e-30) denom4 = 1.0e-30
C	if(denom4.lt.1.0e-5.and.denom4.gt.0.0) denom4 = 1.0e-30
C	if(denom4.gt.-1.0e-5.and.denom4.le.0.0) denom4 = -1.0e-30

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

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

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

	return
	end