anchor.diff:          ! note that the final weights are a mixture of the x and y weights
           bilin_weights=(/x1*f2,x0*f2,x3*f1,x2*f1/)
       end function bilin_weights
  +
  +    !>------------------------------------------------------------
  +    !! Calculate the weights to use for inverse distance interpolation between surrounding points x, y
  +    !!
  +    !! Takes a point (yi,xi) and a set of 4 surrounding points (y[4],x[4]) as input
  +    !!
  +    !! @param   real    yi  input y location to interpolate to.
  +    !! @param   real    y   Array of surrounding y-locations to interpolate from.
  +    !! @param   real    xi  input x location to interpolate to.
  +    !! @param   real    x   Array of surrounding x-locations to interpolate from.
  +    !! @retval  real    weights
  +    !!
  +    !!------------------------------------------------------------
  +    function idw_weights(yi,y,xi,x)
  +        implicit none
  +        real,intent(in)::yi,y(:),xi,x(:)
  +        real, allocatable, dimension(:) :: idw_weights
  +
  +        real, allocatable, dimension(:) :: distances
  +        integer :: n
  +        real    :: total_inverse_distance
  +
  +        n = size(x)
  +        allocate(distances(n))
  +        allocate(idw_weights(n))
  +
  +        distances = 1.0/sqrt((y-yi)**2 + (x-xi)**2)
  +
  +        total_inverse_distance = sum(distances)
  +        ! print*, distances
  +        ! print*, distances / total_inverse_distance
  +
  +        idw_weights = distances / total_inverse_distance
  +
  +    end function idw_weights
   
       !>------------------------------------------------------------
       !! Find the minimum x
  
positive.code:                         lon(k) = lo%lon(xy%x(k), xy%y(k))
                      enddo
                      ! and calculate the weights to apply to each gridcell
                      lo%geolut%w(:,i,j) = bilin_weights(hi%lat(i,j), lat, hi%lon(i,j), lon)
                  endif
              enddo
          enddo
  
      end subroutine geo_LUT
  
      !>------------------------------------------------------------
  
negative.code: 
              do i=1, nx
                  curpos = find_location(lo, hi%lat(i,j), hi%lon(i,j), lastpos)
  
                  if (curpos%x<1) then
                      ! something broke, try assuming that the grids are the same possibly wrapping (for ideal)
                      write(*,*) "Error in Geographic interpolation, check input lat / lon grids", i,j
  