 	{
 		// Force +x-y --> +x+y
 
+		// This is needed because latlon coordinates are created from newbase area
+		// and they are in +x+y. This also means that we have to flip the data
+		// back after interpolation.
+
+		// TODO: Do not flip the data twice, but create the grid in the correct scanning mode!
+		
 		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);
+		Swap <<<gs,bs,0,stream>>>(d_source, sourceInfo->size_x, sourceInfo->size_y);
 
 		sourceInfo->j_scans_positive = true;
 	}

	targetGrid.Reset();
	
	int i = 0;
	
	while(targetGrid.Next())
	{
		NFmiPoint latlon = targetGrid.LatLon();
		NFmiPoint gp = sourceGrid.LatLonToGrid(targetGrid.LatLon());

#ifdef EXTRADEBUG
		if (!sourceArea->IsInside(latlon))
		{
			std::cout << "Latlon " << latlon << " is outside source area!" << std::endl;


	// Do bilinear transform on CUDA device
	InterpolateCudaKernel <<<gs,bs,0,stream>>>(d_source, d_target, d_grid, *sourceInfo, *targetInfo);

	if (!targetInfo->j_scans_positive)
	{
		// Flip data back

		size_t N = targetInfo->size_x * targetInfo->size_y * 0.5 ;

		int bs = 256;
		int gs = N/bs + (N % bs == 0?0:1);
		Swap <<<gs,bs,0,stream>>>(d_target, targetInfo->size_x, targetInfo->size_y);
	}

	CUDA_CHECK(cudaStreamSynchronize(stream));
	
	himan::ReleaseInfo(sourceInfo);
