binarysearch Function

  function binarysearch(length, array, value, delta)
    ! Given an array and a value, returns the index of the element that
    ! is closest to, but less than, the given value.
    ! Uses a binary search algorithm.
    ! "delta" is the tolerance used to determine if two values are equal
    ! if ( abs(x1 - x2) <= delta) then
    !    assume x1 = x2
    ! endif
    ! Should never return an index < 1 or > length!

    implicit none
    integer, intent(in) :: length
    real(dp), dimension(length), intent(in) :: array
    real(dp), intent(in) :: value
    real(dp), intent(in), optional :: delta
    integer :: binarysearch
    integer :: left, middle, right
    real(dp) :: d

    if (present(delta) .eqv. .true.) then
        d = delta
    else
        d = 1e-9
    endif

    left = 1
    right = length
    do
        if (left > right) then
            exit
        endif
        middle = nint((real(left+right)) / 2.0)
        if ( abs(array(middle) - value) <= d) then
            binarySearch = middle
            return
        else if (array(middle) > value) then
            right = middle - 1
        else
            left = middle + 1
        end if
    end do
    binarysearch = min(max(right,1),length)

  end function binarysearch