 	{
 		ret = Linear(dist.x, av, bv);
 	}
-	
+
 	else if (fabs(dist.x) < kEpsilon && cv != kFloatMissing && av != kFloatMissing)
 	{
-		ret = Linear(dist.y, cv, dv);
+		ret = Linear(dist.y, cv, av);
 	}
 
 	else if (fabs(dist.x - 1) < kEpsilon && av != kFloatMissing && bv != kFloatMissing)

}

__device__
double BiLinearInterpolation(const double* __restrict__ d_source, himan::info_simple& sourceInfo, const point& gp)
{
	// 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));

	// Neighbor values

	double av = d_source[Index(a,sourceInfo.size_x)];
	double bv = d_source[Index(b,sourceInfo.size_x)];
	double cv = d_source[Index(c,sourceInfo.size_x)];
	double dv = d_source[Index(d,sourceInfo.size_x)];	

	// Distance of interpolated point to neighboring points

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

	double ret = kFloatMissing;

	// All values present, regular bilinear interpolation
	
	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);

		printf("More than one point missing for gp x: %f y:%f --> a:%f b:%f c:%f d:%f | dx:%f dy:%f\n", gp.x, gp.y, av, bv, cv, dv, dist.x, dist.y);
	}

	if (ret == kFloatMissing) { 
	
		printf("gpx:%f gpy:%f [%ld %ld] |  dist x:%f y:%f\n", gp.x, gp.y, size_x, size_y, dist.x, dist.y);
		printf("av:%f bv:%f cv:%f dv:%f | interp:%f\n", av, bv, cv, dv, ret);
		printf("ax:%f ay:%f bx:%f by:%f cx:%f cy:%f dx:%f dy:%f\n", a.x, a.y, b.x, b.y, c.x, c.y, d.x, d.y);
		printf("is inside grid: a:%d b:%d c:%d d:%d\n", IsInsideGrid(a, size_x, size_y), IsInsideGrid(b, size_x, size_y), IsInsideGrid(c, size_x, size_y), IsInsideGrid(d, size_x, size_y));
	
	}
	
#endif

	return ret;
}

__device__
double NearestPointValueInterpolation(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));
