dbspvn Subroutine

    subroutine dbspvn(t,jhigh,k,index,x,ileft,vnikx,work,iwork,iflag)

    implicit none

    real(wp),intent(in)  :: t(*)     !! knot vector of length n+k, where
                                     !! n = number of b-spline basis functions
                                     !! n = sum of knot multiplicities-k
                                     !! dimension t(ileft+jhigh)
    integer,intent(in)   :: jhigh    !! order of b-spline, 1 <= jhigh <= k
    integer,intent(in)   :: k        !! highest possible order
    integer,intent(in)   :: index    !! index = 1 gives basis functions of order jhigh
                                     !!       = 2 denotes previous entry with work, iwork
                                     !!         values saved for subsequent calls to
                                     !!         dbspvn.
    real(wp),intent(in)  :: x        !! argument of basis functions, t(k) <= x <= t(n+1)
    integer,intent(in)   :: ileft    !! largest integer such that t(ileft) <= x < t(ileft+1)
    real(wp),intent(out) :: vnikx(k) !! vector of length k for spline values.
    real(wp),intent(out) :: work(*)  !! a work vector of length 2*k
    integer,intent(out)  :: iwork    !! a work parameter.  both work and iwork contain
                                     !! information necessary to continue for index = 2.
                                     !! when index = 1 exclusively, these are scratch
                                     !! variables and can be used for other purposes.
    integer,intent(out) :: iflag     !!   0: no errors
                                     !! 201: k does not satisfy k>=1
                                     !! 202: jhigh does not satisfy 1<=jhigh<=k
                                     !! 203: index is not 1 or 2
                                     !! 204: x does not satisfy t(ileft)<=x<=t(ileft+1)

    integer :: imjp1, ipj, jp1, jp1ml, l
    real(wp) :: vm, vmprev

    ! content of j, deltam, deltap is expected unchanged between calls.
    ! work(i) = deltap(i),
    ! work(k+i) = deltam(i), i = 1,k

    if (k<1) then
        write(error_unit,'(A)') 'dbspvn - k does not satisfy k>=1'
        iflag = 201
        return
    endif
    if (jhigh>k .or. jhigh<1) then
        write(error_unit,'(A)') 'dbspvn - jhigh does not satisfy 1<=jhigh<=k'
        iflag = 202
        return
    endif
    if (index<1 .or. index>2) then
        write(error_unit,'(A)') 'dbspvn - index is not 1 or 2'
        iflag = 203
        return
    endif
    if (x<t(ileft) .or. x>t(ileft+1)) then
        write(error_unit,'(A)') 'dbspvn - x does not satisfy t(ileft)<=x<=t(ileft+1)'
        iflag = 204
        return
    endif

    iflag = 0

    if (index==1) then
        iwork = 1
        vnikx(1) = 1.0_wp
        if (iwork>=jhigh) return
    endif

    do
        ipj = ileft + iwork
        work(iwork) = t(ipj) - x
        imjp1 = ileft - iwork + 1
        work(k+iwork) = x - t(imjp1)
        vmprev = 0.0_wp
        jp1 = iwork + 1
        do l=1,iwork
            jp1ml = jp1 - l
            vm = vnikx(l)/(work(l)+work(k+jp1ml))
            vnikx(l) = vm*work(l) + vmprev
            vmprev = vm*work(k+jp1ml)
        enddo
        vnikx(jp1) = vmprev
        iwork = jp1
        if (iwork>=jhigh) exit
    enddo

    endsubroutine dbspvn