db2val Subroutine

    subroutine db2val(xval,yval,idx,idy,tx,ty,nx,ny,kx,ky,bcoef,f,iflag,inbvx,inbvy,iloy)

    implicit none

    integer,intent(in)                   :: idx      !! x derivative of piecewise polynomial to evaluate.
    integer,intent(in)                   :: idy      !! y derivative of piecewise polynomial to evaluate.
    integer,intent(in)                   :: nx       !! the number of interpolation points in x. (same as in last call to [[db2ink]])
    integer,intent(in)                   :: ny       !! the number of interpolation points in y. (same as in last call to [[db2ink]])
    integer,intent(in)                   :: kx       !! order of polynomial pieces in x. (same as in last call to [[db2ink]])
    integer,intent(in)                   :: ky       !! order of polynomial pieces in y. (same as in last call to [[db2ink]])
    real(wp),intent(in)                  :: xval     !! x coordinate of evaluation point.
    real(wp),intent(in)                  :: yval     !! y coordinate of evaluation point.
    real(wp),dimension(nx+kx),intent(in) :: tx       !! sequence of knots defining the piecewise polynomial in the x direction. (same as in last call to [[db2ink]])
    real(wp),dimension(ny+ky),intent(in) :: ty       !! sequence of knots defining the piecewise polynomial in the y direction. (same as in last call to [[db2ink]])
    real(wp),dimension(nx,ny),intent(in) :: bcoef    !! the b-spline coefficients computed by [[db2ink]].
    real(wp),intent(out)                 :: f        !! interpolated value
    integer,intent(out)                  :: iflag    !! status flag: 0 : no errors, /=0 : error
    integer,intent(inout)                :: inbvx    !! initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
    integer,intent(inout)                :: inbvy    !! initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
    integer,intent(inout)                :: iloy     !! initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.

    integer :: k, lefty, mflag, kcol
    real(wp),dimension(ky) :: temp
    real(wp),dimension(3*max(kx,ky)) :: work

    f = 0.0_wp

    if (xval<tx(1) .or. xval>tx(nx+kx)) then
        write(error_unit,'(A)') 'db2val - x value out of bounds'
        iflag = 1
        return
    endif
    if (yval<ty(1) .or. yval>ty(ny+ky)) then
        write(error_unit,'(A)') 'db2val - y value out of bounds'
        iflag = 2
        return
    endif

    iflag = -1
    call dintrv(ty,ny+ky,yval,iloy,lefty,mflag); if (mflag /= 0) return

    kcol = lefty - ky
    do k=1,ky
        kcol = kcol + 1
        temp(k) = dbvalu(tx,bcoef(:,kcol),nx,kx,idx,xval,inbvx,work,iflag)
        if (iflag/=0) return !error
    enddo

    kcol = lefty - ky + 1
    f = dbvalu(ty(kcol:),temp,ky,ky,idy,yval,inbvy,work,iflag)

    endsubroutine db2val