 			
 			// TODO: Maybe add special cases if only one or two neighbors are missing?			
 
-#ifdef DEBUG
+#ifdef EXTRADEBUG
 			// Neighbor point indexes in linear format
 
-			int aidx = index(a,sourceInfo.size_x);
-			int bidx = index(b,sourceInfo.size_x);
-			int cidx = index(c,sourceInfo.size_x);
-			int didx = index(d,sourceInfo.size_x);
-
-			/*
-			if (interp < 100)
+			int aidx = Index(a,sourceInfo.size_x);
+			int bidx = Index(b,sourceInfo.size_x);
+			int cidx = Index(c,sourceInfo.size_x);
+			int didx = Index(d,sourceInfo.size_x);
+			
+			if (i == 0 && j == 0)
 			{
-				printf("x:%d y:%d gpx:%f gpy:%f\n", x, y, gp.x, gp.y);
+				printf("x:%d y:%d gpx:%f gpy:%f\n", i, j, gp.x, gp.y);
 				printf("a x:%d y:%d val:%f\n", int(a.x), int(a.y), av);
 				printf("b x:%d y:%d val:%f\n", int(b.x), int(b.y), bv);
 				printf("c x:%d y:%d val:%f\n", int(c.x), int(c.y), cv);
 				printf("d x:%d y:%d val:%f\n", int(d.x), int(d.y), dv);
 				printf("dist x:%f y:%f\n", dist.x, dist.y);
 				printf("interp:%f\n", interp);
-			} */
+			}
 #endif
 		}
 
-		target[index(x, y, targetInfo.size_x)] = interp ;
-
-	}
-}
+		target[idx] = interp ;
+
+	}
+}
+
 
 bool InterpolateCuda(himan::info_simple* sourceInfo, himan::info_simple* targetInfo)
 {
 	cudaStream_t stream;
 	CUDA_CHECK(cudaStreamCreate(&stream));
-	
-	NFmiArea* sourceArea = CreateArea(sourceInfo);
-	NFmiArea* targetArea = CreateArea(targetInfo);
-	
-	//NFmiPoint interpbl(targetInfo->first_lon, targetInfo->first_lat);
-	//NFmiPoint interptr(interpbl.X() + (targetInfo->size_x - 1) * targetInfo->di, interpbl.Y() + (targetInfo->size_y - 1) * targetInfo->dj);
-
-	NFmiPoint interpbl = targetArea->BottomLeftLatLon();
-	NFmiPoint interptr = targetArea->TopRightLatLon();
-		
-	NFmiGrid sourceGrid(sourceArea, sourceInfo->size_x, sourceInfo->size_y);
-
-	NFmiPoint gpbl = sourceGrid.LatLonToGrid(interpbl);
-
-#ifdef DEBUG
-	NFmiPoint gptr = sourceGrid.LatLonToGrid(interptr);
-#endif
-
-	if (targetInfo->projection == himan::kRotatedLatLonProjection)
-	{
-		interpbl = dynamic_cast<NFmiRotatedLatLonArea*> (targetArea)->ToRotLatLon(interpbl);
-		interptr = dynamic_cast<NFmiRotatedLatLonArea*> (targetArea)->ToRotLatLon(interptr);
-	}
-	
-	point offset(gpbl.X(), gpbl.Y());
-
-#ifdef DEBUG
-	std::cout	<< "Source area BL: " << sourceArea->BottomLeftLatLon()
-				<< "Source area TR: " << sourceArea->TopRightLatLon()
-				<< "Source grid BL: " << sourceGrid.LatLonToGrid(sourceArea->BottomLeftLatLon())
-				<< "Source grid TR: " << sourceGrid.LatLonToGrid(sourceArea->TopRightLatLon())
-				<< "Target area BL: " << interpbl
-				<< "Target area TR: " << interptr
-				<< "Target grid BL: " << gpbl
-				<< "Target grid TR: " << gptr
-				<< "Source grid dx/dy " << sourceInfo->di << "/" << sourceInfo->di << std::endl
-				<< "Target grid dx/dy " << targetInfo->di << "/" << targetInfo->dj << std::endl
-				<< "Offset x/y: " << offset.x << "/" << offset.y << std::endl
-				;
-#endif	
-	delete (sourceArea);
-	delete (targetArea);
+
+	/* Determine all grid point coordinates that need to be interpolated.
+	 * This is done with newbase by explicitly looping through the grid.
+	 * Initially I tried to implement it with just starting point and offset
+	 * but the code was awkward and would not work with stereographic projections
+	 * anyway.
+	 */
+
+	point* grid = CreateGrid(sourceInfo, targetInfo);
+	
+	const size_t N = targetInfo->size_x * targetInfo->size_y;
+	
+	point* d_grid = 0;
+	CUDA_CHECK(cudaMalloc((void**) &d_grid, sizeof(point) * N));
+	CUDA_CHECK(cudaMemcpyAsync(d_grid, grid, sizeof(point) * N, cudaMemcpyHostToDevice, stream));
 
 	double* d_source = 0;
 	double* d_target = 0;
 
 	CUDA_CHECK(cudaMalloc((void **) &d_source, sourceInfo->size_x * sourceInfo->size_y * sizeof(double)));
-	CUDA_CHECK(cudaMalloc((void **) &d_target, targetInfo->size_x * targetInfo->size_y * sizeof(double)));
+	CUDA_CHECK(cudaMalloc((void **) &d_target, N * sizeof(double)));
+
 #ifdef DEBUG
 	CUDA_CHECK(cudaMemset(d_target, 0, targetInfo->size_x * targetInfo->size_y * 8));
 #endif


};

 NFmiArea* CreateArea(himan::info_simple* info)
{
	NFmiPoint bl(info->first_lon, info->first_lat);
	NFmiPoint tr(bl.X() + (info->size_x - 1) * info->di, bl.Y() + (info->size_y - 1) * info->dj);

	NFmiArea* area = 0;
	
	if (info->projection == himan::kLatLonProjection)
	{
		area = new NFmiLatLonArea(bl, tr);

__device__ 
double BiLinear(double dx, double dy, double a, double b, double c, double d)
{
	// Method below is faster but gives visible interpolation artifacts

	//double ab = Linear(dx, a, b);
	//double cd = Linear(dx, c, d);
	//return Linear(dy, ab, cd);

	// This one gives smooth interpolation surfaces
	return (1 - dx) * (1 - dy) * c + dx * (1 - dy) * d + (1 - dx) * dy * a + dx * dy * b;
}

__global__ 
void InterpolateCudaKernel(const double* __restrict__ source, 
							double* __restrict__ target,
							const point* __restrict__ d_grid,
							himan::info_simple sourceInfo, 
							himan::info_simple targetInfo)
{

	// idx is our pointer to the TARGET data in linear format

	const int idx = blockIdx.x * blockDim.x + threadIdx.x;

	if (idx < targetInfo.size_x * targetInfo.size_y)
	{
		// next we need to get x and y of the 'idx' in the source grid coordinates
		// to do that we first determine the i and j of the target grid coordinates
		
		const int i = fmod(static_cast<double> (idx), static_cast<double> (targetInfo.size_x));
		const int j = floor(static_cast<double> (idx / targetInfo.size_x));

		// 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 (gp.x < 0 || gp.x > sourceInfo.size_x - 1)
		{
			// wrap x if necessary
			// this might happen f.ex. with EC where grid start at 0 meridian and 
			// we interpolate from say -10 to 40 longitude

			while (gp.x < 0) gp.x += sourceInfo.size_x;
			while (gp.x > sourceInfo.size_x-1) gp.x -= sourceInfo.size_x-1;
		}
		
		double interp = kFloatMissing;
		
