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:         real(dp):: x(5), y(5), xout(9), yout_true(9), yout(9)
          type(linear_interpolator_type):: li
          integer(ip):: i, j
          logical:: copy_data
  
          ! Data to interpolate from
          x = [1.0_dp, 4.0_dp, 4.01_dp, 10.0_dp, 11.0_dp]
  
negative.code:     !!
      use global_mod, only: dp, ip
      use stop_mod, only: generic_stop
      use qsort_mod, only: sort ! Only for testing
      implicit none
  
  
  