dbvalu Function

    real(wp) function dbvalu(t,a,n,k,ideriv,x,inbv,work,iflag)

    implicit none

    integer,intent(in)               :: n       !! number of b-spline coefficients.
                                                !! (sum of knot multiplicities-k)
    real(wp),dimension(:),intent(in) :: t       !! knot vector of length n+k
    real(wp),dimension(n),intent(in) :: a       !! b-spline coefficient vector of length n
    integer,intent(in)               :: k       !! order of the b-spline, k >= 1
    integer,intent(in)               :: ideriv  !! order of the derivative, 0 <= ideriv <= k-1.
                                                !! ideriv = 0 returns the b-spline value
    real(wp),intent(in)              :: x       !! argument, t(k) <= x <= t(n+1)
    integer,intent(inout)            :: inbv    !! an initialization parameter which must be set
                                                !! to 1 the first time dbvalu is called.
                                                !! inbv contains information for efficient process-
                                                !! ing after the initial call and inbv must not
                                                !! be changed by the user.  distinct splines require
                                                !! distinct inbv parameters.
    real(wp),dimension(:)            :: work    !! work vector of length 3*k
    integer,intent(out)              :: iflag   !!   0: no errors
                                                !! 401: k does not satisfy k>=1
                                                !! 402: n does not satisfy n>=k
                                                !! 403: ideriv does not satisfy 0<=ideriv<k
                                                !! 404: x is not greater than or equal to t(k)
                                                !! 405: x is not less than or equal to t(n+1)
                                                !! 406: a left limiting value cannot be obtained at t(k)

    integer :: i,iderp1,ihi,ihmkmj,ilo,imk,imkpj,ipj,&
               ip1,ip1mj,j,jj,j1,j2,kmider,kmj,km1,kpk,mflag
    real(wp) :: fkmj

    dbvalu = 0.0_wp

    if (k<1) then
        write(error_unit,'(A)') 'dbvalu - k does not satisfy k>=1'
        iflag = 401
        return
    endif

    if (n<k) then
        write(error_unit,'(A)') 'dbvalu - n does not satisfy n>=k'
        iflag = 402
        return
    endif

    if (ideriv<0 .or. ideriv>=k) then
        write(error_unit,'(A)') 'dbvalu - ideriv does not satisfy 0<=ideriv<k'
        iflag = 403
        return
    endif

    kmider = k - ideriv

    ! find *i* in (k,n) such that t(i) <= x < t(i+1)
    ! (or, <= t(i+1) if t(i) < t(i+1) = t(n+1)).

    km1 = k - 1
    call dintrv(t, n+1, x, inbv, i, mflag)
    if (x<t(k)) then
        write(error_unit,'(A)') 'dbvalu - x is not greater than or equal to t(k)'
        iflag = 404
        return
    endif

    if (mflag/=0) then

        if (x>t(i)) then
            write(error_unit,'(A)') 'dbvalu - x is not less than or equal to t(n+1)'
            iflag = 405
            return
        endif

        do
            if (i==k) then
                write(error_unit,'(A)') 'dbvalu - a left limiting value cannot be obtained at t(k)'
                iflag = 406
                return
            endif
            i = i - 1
            if (x/=t(i)) exit
        enddo

    endif

    ! difference the coefficients *ideriv* times
    ! work(i) = aj(i), work(k+i) = dp(i), work(k+k+i) = dm(i), i=1.k

    imk = i - k
    do j=1,k
        imkpj = imk + j
        work(j) = a(imkpj)
    enddo

    if (ideriv/=0) then
        do j=1,ideriv
            kmj = k - j
            fkmj = real(kmj,wp)
            do jj=1,kmj
                ihi = i + jj
                ihmkmj = ihi - kmj
                work(jj) = (work(jj+1)-work(jj))/(t(ihi)-t(ihmkmj))*fkmj
            enddo
        enddo
    endif

    ! compute value at *x* in (t(i),(t(i+1)) of ideriv-th derivative,
    ! given its relevant b-spline coeff. in aj(1),...,aj(k-ideriv).

    if (ideriv/=km1) then
        ip1 = i + 1
        kpk = k + k
        j1 = k + 1
        j2 = kpk + 1
        do j=1,kmider
            ipj = i + j
            work(j1) = t(ipj) - x
            ip1mj = ip1 - j
            work(j2) = x - t(ip1mj)
            j1 = j1 + 1
            j2 = j2 + 1
        enddo
        iderp1 = ideriv + 1
        do j=iderp1,km1
            kmj = k - j
            ilo = kmj
            do jj=1,kmj
                work(jj) = (work(jj+1)*work(kpk+ilo)+work(jj)*&
                            work(k+jj))/(work(kpk+ilo)+work(k+jj))
                ilo = ilo - 1
            enddo
        enddo
    endif

    iflag = 0
    dbvalu = work(1)

    endfunction dbvalu