!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
SUBROUTINE splint(xa,ya,y2a,n,x,y)
!***********************************************************************
!                                                                      *
!        Given tabulated FUNCTION ya(xa) defined at n points, and its  *
!        second order derivatives y2a as derived from spline, and      *
!        given a value of x, this routine RETURNs a cubic-spline       *
!        interpolated value y.                                         *
!                                                                      *
!        From:      "Numerical recepies", p. 89                         *
!***********************************************************************
implicit none
integer::n,kl0,khi,k
real::x,y,h,a,b
REAL, DIMENSION (n)::  xa, ya, y2a
!-----------------------------------------------------------------------
      kl0=1.
      khi=n

!------------- Obsolet-------------------------------------------------
!
! 1      IF(khi-kl0.gt.1) THEN
!          k=(khi+kl0)/2
!          IF(xa(k).gt.x) THEN
!            khi=k
!          ELSE
!            kl0=k
!          END IF
!        goto 1
!        END IF
!-----------------------------------------------------------------------
 loop:  DO 
          IF(khi-kl0.le.1) EXIT loop
            k=(khi+kl0)/2
            IF(xa(k).gt.x) THEN
              khi=k
            ELSE
              kl0=k
            END IF
        END DO loop
!-----------------------------------------------------------------------

      h=xa(khi)-xa(kl0)
      IF(h.eq.0.) THEN
         print *,xa(khi),xa(kl0),khi,kl0
         PRINT*, 'bad xa input'
         stop
      END IF

      a=(xa(khi)-x)/h
      b=(x-xa(kl0))/h

      y=a*ya(kl0)+b*ya(khi)   + &
        ((a**3-a)*y2a(kl0)+(b**3-b)*y2a(khi))*(h**2)/6
        
!-----------------------------------------------------------------------
END SUBROUTINE splint