 	PrepareInfo(sourceInfo, d_source, stream);
 	PrepareInfo(targetInfo);
 
-	const int blocksize = 16;
-	const dim3 blocks(blocksize,blocksize);
-	
-	dim3 grids;
-	
-	grids.x = ceil((targetInfo->size_x+blocks.x-1)/blocks.x);
-	grids.y = ceil((targetInfo->size_y+blocks.y-1)/blocks.y);
+	if (!sourceInfo->j_scans_positive)
+	{
+		// Force +x-y --> +x+y
+
+		size_t N = sourceInfo->size_x * sourceInfo->size_y * 0.5 ;
+
+		int bs = 256;
+		int gs = N/bs + (N % bs == 0?0:1);
+		Flip <<<gs,bs,0,stream>>>(d_source, sourceInfo->size_x, sourceInfo->size_y);
+
+		sourceInfo->j_scans_positive = true;
+	}
+	
+	const int bs = 256;
+	const int gs = N/bs + (N % bs == 0?0:1);
 
 	// Do bilinear transform on CUDA device
-	InterpolateCudaKernel <<<grids,blocks,0,stream>>>(d_source, d_target, *sourceInfo, *targetInfo, offset);
+	InterpolateCudaKernel <<<gs,bs,0,stream>>>(d_source, d_target, d_grid, *sourceInfo, *targetInfo);
  
 	CUDA_CHECK(cudaStreamSynchronize(stream));
-
+	
 	himan::ReleaseInfo(sourceInfo);
 	himan::ReleaseInfo(targetInfo, d_target, stream);
 	
-	CUDA_CHECK(cudaStreamSynchronize(stream));
-
 	CUDA_CHECK(cudaFree(d_source));
 	CUDA_CHECK(cudaFree(d_target));
-
+	CUDA_CHECK(cudaFree(d_grid));
+	
 	CUDA_CHECK(cudaStreamDestroy(stream));
-
+	
 	return true;
 }
			point d(ceil(gp.x), floor(gp.y));

			// Neighbor values
			
			double av = source[index(a,sourceInfo.size_x)];
			double bv = source[index(b,sourceInfo.size_x)];
			double cv = source[index(c,sourceInfo.size_x)];
			double dv = source[index(d,sourceInfo.size_x)];

			// Distance of interpolated point to neighboring points
			
			point dist(gp.x - c.x, gp.y - c.y);

			if (av != kFloatMissing && bv != kFloatMissing && cv != kFloatMissing && dv != kFloatMissing)

		NFmiPoint sp(info->south_pole_lon, info->south_pole_lat);
		area = new NFmiRotatedLatLonArea(bl, tr, sp, NFmiPoint(0,0), NFmiPoint(1,1), true);
	}
	else if (info->projection == himan::kStereographicProjection)
	{
		area = new NFmiStereographicArea(bl, (info->size_x - 1) * info->di, (info->size_y - 1) * info->dj, info->orientation);
	}
	else
	{
		throw std::runtime_error("Invalid projection for cuda interpolation");
	}

	assert(area);
	return area;
}

point* CreateGrid(himan::info_simple* sourceInfo, himan::info_simple* targetInfo)
{

	NFmiArea* sourceArea = CreateArea(sourceInfo);
	NFmiArea* targetArea = CreateArea(targetInfo);
	
	NFmiGrid sourceGrid(sourceArea, sourceInfo->size_x, sourceInfo->size_y, kBottomLeft);
	NFmiGrid targetGrid(targetArea, targetInfo->size_x, targetInfo->size_y, kBottomLeft);
/*
	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())
				<< "Source J scans positive: " << sourceInfo->j_scans_positive << std::endl
				<< "Target area BL: " << targetArea->BottomLeftLatLon()
				<< "Target area TR: " << targetArea->TopRightLatLon()
				<< "Target grid BL (relative): " << sourceGrid.LatLonToGrid(targetArea->BottomLeftLatLon())
				<< "Target grid TR (relative): " << sourceGrid.LatLonToGrid(targetArea->TopRightLatLon())
				<< "Target J scans positive: " << targetInfo->j_scans_positive << std::endl
				;
*/	
	delete (sourceArea);
	delete (targetArea);
	
	point* ret = new point[targetGrid.XNumber() * targetGrid.YNumber()];
	
	targetGrid.Reset();
	
	int i = 0;
	
	while(targetGrid.Next())
	{
		NFmiPoint gp = sourceGrid.LatLonToGrid(targetGrid.LatLon());
		
		ret[i].x = gp.X();
		ret[i].y = gp.Y();
		
		i++;
	}

	return ret;
}

__device__
