anchor.diff:      ! Suppose x is a SORTED vector of length n with x(i) <= x(i+1).
       ! We want to find the index corresponding to the value in 'x' that is nearest y.
       ! This is a useful operation e.g. for interpolation of sorted input data
  -    !
       pure subroutine nearest_index_sorted(n, x, y, output)
           integer(ip), intent(in) :: n
           real(dp), intent(in) :: x(n)
           real(dp), intent(in) :: y
  -        integer(ip), intent(out) :: output
  +        integer(ip), intent(inout) :: output
   
           integer :: i, upper, lower
   
  -        if(y < x(1)) then
  +        if(y <= x(1)) then
  +            ! Quick exit
               output = 1
  +        else if (y >= x(n)) then
  +            ! Quick exit
  +            output = n
           else
  -            if (y > x(n)) then
  -                output = n
  -            else
  -                lower = 1
  -                upper = n
  -                ! Deliberate integer division
  -                i = (lower + upper)/2 !floor(0.5_dp*(lower + upper))
  -                do while ((x(i) > y).or.(x(i+1) < y))
  -                   if(x(i) > y) then
  -                       upper = i
  -                   else
  -                       lower = i
  -                   end if
  -                   ! Deliberate integer division
  -                   i = (lower + upper)/2 !floor(0.5_dp*(lower + upper))
  -                end do
  -
  -                if ( y - x(i) > x(i+1) - y) then
  -                    output = i+1
  -                else
  -                    output = i
  -                end if
  -
  -            end if
  +            lower = 1
  +            upper = n
  +
  +            ! Guess based on the prior value of output (quickly brackets if
  +            ! the previous y-value was close)
  +            if(output >= 1) then
  +                if(output <= n .and. output >= 2) then
  +                    if(y >= x(output - 1)) lower = output - 1
  +                end if
  +                if(output <= n - 1) then
  +                    if(y <= x(output + 1)) upper = output + 1
  +                end if
  +            end if
  +
  +            ! Deliberate integer division
  +            i = (lower + upper)/2 
  +            do while ((x(i) > y).or.(x(i+1) < y))
  +               if(x(i) > y) then
  +                   upper = i
  +               else
  +                   lower = i
  +               end if
  +               ! Deliberate integer division
  +               i = (lower + upper)/2 !floor(0.5_dp*(lower + upper))
  +            end do
  +
  +            if ( y - x(i) > x(i+1) - y) then
  +                output = i+1
  +            else
  +                output = i
  +            end if
  +
           end if
   
       end subroutine
  
positive.code: 
          end do
  
      end subroutine
  
  end module
  
negative.code:     !!
      use global_mod, only: dp, ip
      use stop_mod, only: generic_stop
      use qsort_mod, only: sort ! Only for testing
      implicit none
  
  
  