SUBROUTINE BSPLNK(X,K,B,L,T)
*     Returned are B(1)...B(K), the B-spline values for given X and order K 
*     K = order of spline (1=const, 2=linear...) up to limit K=KMAX=10 
*     The array T(KNK) contains the knot positions.
*     In T(K) =<  x  < T(KNK+1-K) the K B-spline values add up to 1. 
*     Multiple knots allowed: (K - mult)-th derivative is discontinuous.
*     Multiple knots allow to define a B-spline sum with certain
*     discontinuities.
* 
*     The B-spline values  B(1) ... B(K)      correspond to the B-splines
*                          B(L) ... B(L-1+K)
*     In total there are  NKN-K  B-splines. The index L can be determined
*     by the call
*        L=MAX(K,MIN(LEFT(X,T,NKN),NKN-K))
*     before the BSPLNK call.
*
*     Example: K = 4 = cubic B-splines
*     X-limits A = T(4) and B = T(KNK-3)
*     Knots T(1) ... T(3) are <= T(4), often defined with equidistant 
*     knot distances; Knots T(KNK-2),T(KNK-1),T(KNK) are >= T(KNK-3),
*     often with equidistant knot distances.
*     The knot multiplicity at knot T(I) is 1, if T(I-1) < T(I) < T(I+1);
*     derivatives 1, 2 and 3 (= K-1) agree at a knot of multiplicity 1.
      REAL B(K),T(*)
      PARAMETER (KMAX=10)
      DOUBLE PRECISION DL(KMAX-1),DR(KMAX-1),SV,TR
*     ...
      IF(K.GT.KMAX) STOP ' BSPLNK: K too large '
      B(1)=1.0
      DO J=1,K-1
       DR(J)=T(L+J)-X
       DL(J)=X-T(L-J+1)
       SV=0.0D0
       DO I=1,J
        TR  =B(I)/(DR(I)+DL(J-I+1))
        B(I)=SV+TR*DR(I)
        SV  =   TR*DL(J-I+1)
       END DO
       B(J+1)=SV
       IF(B(  1)+1.0.EQ.1.0) B(  1)=0.0 ! avoid underflows
       IF(B(J+1)+1.0.EQ.1.0) B(J+1)=0.0
      END DO
      END