anchor.diff:          class(domain_type), target, intent(inout):: domain
   
           integer(ip):: i, j, k
  -        real(dp), allocatable:: x(:,:), y(:,:)
  -
  -        character(len=charlen):: input_elevation_file, polygon_filename
  +        real(dp), allocatable:: x(:), y(:)
  +
  +        ! For reading background elevation
  +        character(len=charlen):: input_elevation_file
  +        type(multi_raster_type) :: elevation_data
   
           ! Add this elevation to all points in houses/ csv polygons
           real(dp), parameter:: house_height = 3.0_dp
  
positive.code:         domain%U(domain%nx(1), :, ELV) = domain%U(domain%nx(1)-1, :, ELV)
          domain%U(:, domain%nx(2), ELV) = domain%U(:, domain%nx(2)-1, ELV) 
  
  
          !
          ! Set friction
          !
  
negative.code:         call read_csv_into_array(polygon_coords, polygon_filename)
          inside_point_counter = 0
          do j = 1, domain%nx(2)
              y = domain%y(j)
              call points_in_poly(polygon_coords(1,:), polygon_coords(2,:), x, y, is_inside_poly)
              inside_point_counter = inside_point_counter + count(is_inside_poly)
  
              do i = 1, domain%nx(1) 
  