SUBROUTINE spline(xb,yb,n,yp1,ypn,y2)

USE Numeric_kinds_Module
USE CommonInfo_Module

IMPLICIT NONE

! Only need to call once

INTEGER,INTENT(IN):: n
INTEGER i,k
REAL(KIND=dp) yp1,ypn
REAL(KIND=dp) p,qn,sig,un
REAL(KIND=dp) xb(n),yb(n),y2(n),u(n)

IF (yp1.gt.0.99d30) THEN
  y2(1)=0.d0
  u(1)=0.d0
ELSE
!  y2(1)=-0.5d0 ! Why should this not be zero?
 y2(1)=0.d0
  u(1)=(3.d0/(xb(2)-xb(1)))*((yb(2)-yb(1))/(xb(2)-xb(1))-yp1)
ENDIF

DO i=2,n-1
  sig=(xb(i)-xb(i-1))/(xb(i+1)-xb(i-1))
  p=sig*y2(i-1)+2.d0
  y2(i)=(sig-1.d0)/p
  u(i)=(6.d0*((yb(i+1)-yb(i))/(xb(i+1)-xb(i))-(yb(i)-yb(i-1))/(xb(i)-xb(i-1)))/(xb(i+1)-xb(i-1))-sig*u(i-1))/p
ENDDO

IF (ypn.gt.0.99d30) THEN
  qn=0.d0
  un=0.d0
ELSE
!  qn=-0.5d0 ! How about this one too?
  qn=-0.5d0
  un=(3.d0/(xb(n)-xb(n-1)))*(ypn-(yb(n)-yb(n-1))/(xb(n)-xb(n-1)))
ENDIF

y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)

DO k=n-1,1,-1
  y2(k)=y2(k)*y2(k+1)+u(k)
ENDDO

END SUBROUTINE spline