anchor.diff:  
           end do
   
  +        ! Test nearest_index_sorted
  +        call random_number(r_x)
  +        call sort(r_x, size(r_x))
  +        call random_number(r_y)
  +        do j = 1, 2 ! Run twice -- the second time we have pre-specified the correct values of i_y(i)
  +            test_failed = .false.
  +            do i = 1, size(r_y)
  +                call nearest_index_sorted(size(r_x, kind=ip), r_x, r_y(i), i_y(i))
  +                ! Check that it really is the nearest
  +                if(any(abs(r_y(i) - r_x) < abs(r_y(i) - r_x(i_y(i))) )) then
  +                    test_failed = .true.
  +                    !print*, r_x(i_y(i)), r_y(i), abs(r_y(i) - r_x(i_y(i))), minval(abs(r_y(i) - r_x))
  +                end if
  +            end do
  +            print*, merge('FAIL', 'PASS', test_failed)
  +        end do
  +
       end subroutine
   
   end module
  
positive.code:     !!
      use global_mod, only: dp, ip
      use stop_mod, only: generic_stop
      implicit none
  
  
  
negative.code:     ! 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(inout) :: output
  
          integer :: i, upper, lower
  
          if(y <= x(1)) then
              ! Quick exit
              output = 1
          else if (y >= x(n)) then
              ! Quick exit
              output = n
          else
              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
  