 		// with i and j we can get the grid point coordinates in the source grid
 		
 		point gp = d_grid[Index(i,j,targetInfo.size_x)];
-		
+
 		if (sourceInfo.wraps_globally && (gp.x < 0 || gp.x > sourceInfo.size_x - 1))
 		{
 			// wrap x if necessary

		
		double interp = kFloatMissing;
		
		if (
			// if interpolated grid points are negative, it means that we are outside
			// of the source area
				
			// sometime first grid point is -0, so we substrace a small value from first
			// grid point accept that value as well
				
			gp.x >= (0 - kEpsilon) && gp.y >= (0 - kEpsilon) &&
				
			// if interpolated grid points are larger than source grid in x or y
			// direction, it means again that we are outside of the area
			
			__double2uint_ru(gp.x) <= sourceInfo.size_x && __double2uint_ru(gp.y) <= sourceInfo.size_y
		)
		{
			
			//targetInfo.interpolation = himan::kNearestPointValue;

}

__device__
bool IsInsideGrid(point& gp, size_t size_x, size_t size_y)
{
	// if interpolated grid points are negative, it means that we are outside the grid
	
	// sometime first grid point is -0, so we subtract a small value from first
	// grid point accept that value as well
	
	if (gp.x >= (0 - kEpsilon) && gp.y >= (0 - kEpsilon) &&
				
		// if interpolated grid points are larger than source grid in x or y
		// direction, it means again that we are outside of the area
			
		((fabs(gp.x - (size_x - 1)) < kEpsilon || __double2uint_ru(gp.x) < size_x) && (fabs(gp.y - (size_y - 1)) < kEpsilon || __double2uint_ru(gp.y) < size_y)))
	{
		return true;
	}

#ifdef EXTRADEBUG
	bool lc = gp.x >= (0 - kEpsilon) && gp.y >= (0 - kEpsilon);
	bool uc = (fabs(gp.x - (size_x - 1)) < kEpsilon || __double2uint_ru(gp.x) < size_x) && (fabs(gp.y - (size_y - 1)) < kEpsilon || __double2uint_ru(gp.y) < size_y);
	
	printf("gp x:%f y:%f discarded [%ld,%ld]: lower cond --> x:%d y:%d upper cond x:%d y:%d\n", 
		gp.x, gp.y, size_x, size_y, gp.x >= (0 - kEpsilon), gp.y >= (0 - kEpsilon), lc, uc);
	
#endif

	return false;
}

__device__
double NearestPointInterpolation(const double* __restrict__ d_source, himan::info_simple& sourceInfo, const point& gp)
{
	int rx = rint(gp.x);
	int ry = rint(gp.y);

	assert(rx >= 0 && rx <= sourceInfo.size_x);
	assert(ry >= 0 && ry <= sourceInfo.size_y);

	return d_source[Index(rx,ry,sourceInfo.size_x)];

}

__device__
double BiLinearInterpolation(const double* __restrict__ d_source, himan::info_simple& sourceInfo, const point& gp)
{
	double ret = kFloatMissing;

	// Find all four neighboring points

	point a(floor(gp.x), ceil(gp.y));
	point b(ceil(gp.x), ceil(gp.y));
	point c(floor(gp.x), floor(gp.y));
	point d(ceil(gp.x), floor(gp.y));
	
	// Assure neighboring points are inside grid and get values
	
	size_t size_x = sourceInfo.size_x;
	size_t size_y = sourceInfo.size_y;
	
	double av = kFloatMissing, bv = kFloatMissing, cv = kFloatMissing, dv = kFloatMissing;
	
	if (IsInsideGrid(a, size_x, size_y))
	{
		av = d_source[Index(a,size_x)];
	}
	if (IsInsideGrid(b, size_x, size_y))
	{
		bv = d_source[Index(b,size_x)];
	}
	if (IsInsideGrid(c, size_x, size_y))
	{
		cv = d_source[Index(c,size_x)];
	}
	if (IsInsideGrid(d, size_x, size_y))
	{
		dv = d_source[Index(d,size_x)];
	}

	// Distance of interpolated point to neighboring points

	point dist(gp.x - c.x, gp.y - c.y);

	assert(dist.x >= 0 && dist.x <= 1);
	assert(dist.y >= 0 && dist.y <= 1);
	
	// If interpolated point is very close to source grid point, pick 
	// the point value directly
	
	// This is preferred since nearest point is faster than bilinear, and
	// if wanted grid point =~ source grid point, the bilinear interpolation
	// value will be very close to nearest point value
	
	if ( 
			(dist.x < kEpsilon || fabs(dist.x-1) < kEpsilon) && 
			(dist.y < kEpsilon || fabs(dist.y-1) < kEpsilon) )
	{
		ret = NearestPointInterpolation(d_source, sourceInfo, gp);
	}
	
	// All values present, regular bilinear interpolation
	
	else if (av != kFloatMissing && bv != kFloatMissing && cv != kFloatMissing && dv != kFloatMissing)
	{
		ret = BiLinear(dist.x, dist.y, av, bv, cv, dv);
	}
	
	// x or y is at grid edge

	else if (fabs(dist.y) < kEpsilon && cv != kFloatMissing && dv != kFloatMissing)
	{
		ret = Linear(dist.x, cv, dv);
