binary_search_real Function

    function binary_search_real(arr, val, thresh) &
        result(pos)

        real(dp), intent(in) :: arr(:)
        real(dp), intent(in) :: val
        real(dp), intent(in) :: thresh
        integer :: pos

        integer :: hi, lo

        ! The search range
        lo = lbound(arr, 1)
        hi = ubound(arr, 1)

        ! Account for poor usage (i.e. array len == 0)
        if(hi < lo) then
            pos = -lo
            return
        endif

        ! Narrow the search range down in steps.
        do while(hi /= lo)
            pos = int(real(hi + lo, dp) / 2_dp)

            if(abs(arr(pos) - val) < thresh) then
                exit
            else if(val > arr(pos)) then
                ! val is "greater" than arr(:len,pos).
                ! The lowest position val can take is hence pos + 1 (i.e. if
                ! val is greater than pos by smaller than pos + 1).
                lo = pos + 1
            else
                ! arr(:,pos) is "greater" than val.
                ! The highest position val can take is hence pos (i.e. if val is
                ! smaller than pos but greater than pos - 1).  This is why
                ! we differ slightly from a standard binary search (where lo
                ! is set to be pos+1 and hi to be pos-1 accordingly), as
                ! a standard binary search assumes that the element you are
                ! searching for actually appears in the array being
                ! searched...
                hi = pos
            endif
        enddo

        ! If we have narrowed down to one position, and it is not the item,
        ! then return -pos to indicate that the item is not present, but that
        ! this is the location it should be in.
        if(hi == lo) then
            if(abs(arr(hi) - val) < thresh) then
                pos = hi
            else if(val > arr(hi)) then
                pos = -hi - 1
            else
                pos = -hi
            endif
        endif

    end function binary_search_real