anchor.diff:          !     f1(unique_Mw + dMw/2) * (unique_Mw + dMw/2 <= max_mw_max) )
           rate_curve_approx_Mws(1:N2) = unique_Mw(1:N2) - dMw(1)*0.5_dp
           rate_curve_approx_Mws(N2+1) = unique_Mw(N2) + dMw(1)*0.5_dp
  -        call linear_interpolation(N1, rate_curve_Mw, rate_curve_exrate, N2+1, rate_curve_approx_Mws, rate_curve_approx)
  +        call linear_interpolation_input_x_increasing(N1, rate_curve_Mw, rate_curve_exrate, N2+1, &
  +            rate_curve_approx_Mws, rate_curve_approx)
           where(rate_curve_approx_Mws > max_mw_max) rate_curve_approx = 0.0_dp
           unique_rates(1:N2) = rate_curve_approx(1:N2) - rate_curve_approx(2:(N2+1))
   
  
positive.code: 
          do i = 1, n_output
              ! set j = nearest index to output_x(i) in input_x
              call nearest_index_sorted(n_input, input_x, output_x(i), j)
              ! interpolate
              if (input_x(j) > output_x(i)) then
                  if(j > 1) then
  
negative.code:     ! 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_x_increasing(n, x, y, output)
          integer(ip), intent(in) :: n
          real(dp), intent(in) :: x(n)
          real(dp), intent(in) :: y      
  