 
 };
 
+__host__ __device__ 
+int Index(int x, int y, int sx, float dir = 1)
+{
+	return y * sx + dir * x;
+}
+
+__host__ __device__ 
+int Index(point p, int sx, float dir = 1)
+{
+	return Index(static_cast<int> (p.x),static_cast<int> (p.y), sx, dir);
+}
+
+__global__
+void Flip(double* __restrict__ arr, size_t ni, size_t nj)
+{
+	const int idx = blockIdx.x * blockDim.x + threadIdx.x;
+
+	// Flip with regards to x axis
+
+	if (idx < nj * ni * 0.5)
+	{
+		size_t half = static_cast<size_t> (floor(static_cast<double>(nj/2)));
+
+		const int i = fmod(static_cast<double> (idx), static_cast<double> (ni));
+		const int j = floor(static_cast<double> (idx / ni));
+
+		double upper = arr[idx];
+		double lower = arr[Index(i,nj-1-j,ni)];
+
+		arr[idx] = lower;
+		arr[Index(i,nj-1-j,ni)] = upper;
+		
+	}
+}
+
+__global__ void Print(double* __restrict__ arr, int i) { printf("%d %f\n", i, arr[i]); }
+
  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);
+	NFmiPoint bl, tr;
+	
+	bl.X(info->first_lon);
+	tr.X(bl.X() + (info->size_x - 1) * info->di);
+
+	if (info->j_scans_positive)
+	{	
+		bl.Y(info->first_lat);
+		tr.Y(bl.Y() + (info->size_y - 1) * info->dj);
+	}
+	else
+	{
+		tr.Y(info->first_lat);
+		bl.Y(tr.Y() - (info->size_y - 1) * info->dj);
+	}
 
 	NFmiArea* area = 0;
-	
+
 	if (info->projection == himan::kLatLonProjection)
 	{
 		area = new NFmiLatLonArea(bl, tr);

			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)

			
			// TODO: Maybe add special cases if only one or two neighbors are missing?			

#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 (i == 0 && j == 0)
			{
				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[idx] = interp ;

	}
}


bool InterpolateCuda(himan::info_simple* sourceInfo, himan::info_simple* targetInfo)
{
	cudaStream_t stream;
	CUDA_CHECK(cudaStreamCreate(&stream));

	/* 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, N * sizeof(double)));

#ifdef DEBUG
	CUDA_CHECK(cudaMemset(d_target, 0, targetInfo->size_x * targetInfo->size_y * 8));
#endif
