!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
SUBROUTINE spline(x,y,n,yp1,ypn,y2)
!***********************************************************************
!                                                                      *
!        Falulate second derivatives of tabulated FUNCTION y(x)        *
!        defined at n points, given first derivatives at the           *
!        boundaries (yp1 and ypn).                                     *
!                                                                      *
!        From:      "Numerical Recipes", p. 88                         *
!***********************************************************************
implicit none
INTEGER, PARAMETER :: nmax=100
integer::n,i,k
real::yp1,ypn,qn,un,sig,p
REAL, DIMENSION (n) :: x, y, y2
REAL, DIMENSION (nmax) :: u
!-----------------------------------------------------------------------
      IF(yp1.gt..99e30) THEN
        y2(1)=0.
        u(1)=0.
      ELSE
        y2(1)=-0.5
        u(1)=(3./(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))
        p=sig*y2(i-1)+2.
        y2(i)=(sig-1.)/p
        u(i)=(6.*((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
      END DO

      IF(ypn.gt..99e30) THEN
        qn=0.
        un=0.
      ELSE
        qn=0.5
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      END IF
      
        y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      DO k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
      END DO
!-----------------------------------------------------------------------
END SUBROUTINE spline