anchor.diff:  
           !
           ! Line version
  +        ! (the 'eps' in the line definition helps the test to pass in both single and double precision,
  +        !  otherwise round-off can shift the results).
           !
  -        lx = [-1.0_dp ,  5.5_dp,  5.5_dp]
  -        ly = [3.5_dp  ,  3.5_dp, 11.5_dp]
  +        lx = [-1.0_dp ,  5.5_dp,  5.5_dp] + eps
  +        ly = [3.5_dp  ,  3.5_dp, 11.5_dp] - eps
           lz = lx + ly
           grid = 0.0_dp
           call burn_line_into_grid(lx, ly, lz, grid, lower_left, upper_right)
  -        expected_grid = 0.0
  +        expected_grid = 0.0_dp
           expected_grid(1:6,4) = [4.12500000_dp, 4.93750000_dp,  5.75000000_dp, 7.37500000_dp, 8.18750000_dp, 9.00000000_dp]
           expected_grid(6,5:10) = [9.88888931_dp, 10.7777777_dp, 11.6666670_dp, 13.4444447_dp, 14.3333330_dp, 15.2222223_dp]
           if(all(abs(grid - expected_grid) < 1.0e-05_dp)) then
  
positive.code:                 zip1 = z(i+1)
  
                  ! Interpolate np+1 points along the line segment between points i and (i+1)
                  np = ceiling(max((xip1-xi)/dx(1), (yip1-yi)/dx(2))) + 1
                  do j = 0, np
                      xj = xi + (j*(xip1-xi))/np
                      yj = yi + (j*(yip1-yi))/np
                      zj = zi + (j*(zip1-zi))/np
                      call burn_xyz_into_grid(xj, yj, zj, grid, lower_left, upper_right, burnt)
                  end do
              else
  
negative.code:     !! provides routines to burn these into grids. For example, this can be used to ensure that 
      !! an elevation grid includes a continuous representation of a breakwall or riverwall.
      !!
      use global_mod, only: ip, dp, charlen, force_double
      use logging_mod
      use stop_mod
      use file_io_mod, only: read_csv_into_array
  