 
 ///////////////////////////////////////////////////////////////////////////////
 
+__global__ 
+void sum_pressure (const Real* pres_red, const Real* pres_black, 
+					Real* pres_sum) 
+{
+	int row = (blockIdx.x * blockDim.x) + threadIdx.x + 1;
+	int col = (blockIdx.y * blockDim.y) + threadIdx.y + 1;
+
+	// shared memory for block's sum
+	__shared__ Real sum_cache[BLOCK_SIZE];
+			
+	int NUM_2 = NUM >> 1;
+
+	Real pres_r = pres_red(col, row);
+	Real pres_b = pres_black(col, row);
+
+	// add squared pressure
+	sum_cache[threadIdx.x] = (pres_r * pres_r) + (pres_b * pres_b);
+
+	// synchronize threads in block to ensure all thread values stored
+	__syncthreads();
+
+	// add up values for block
+	int i = BLOCK_SIZE >> 1;
+	while (i != 0) {
+		if (threadIdx.x < i) {
+			sum_cache[threadIdx.x] += sum_cache[threadIdx.x + i];
+		}
+		__syncthreads();
+		i >>= 1;
+	}
+
+	// store block's summed values
+	if (threadIdx.x == 0) {
+		pres_sum[blockIdx.y + (gridDim.y * blockIdx.x)] = sum_cache[0];
+	}
+
+} // end sum_pressure
+
+///////////////////////////////////////////////////////////////////////////////
+__global__ 
+void set_horz_pres_BCs (Real* pres_red, Real* pres_black) 
+{
+	int col = (blockIdx.x * blockDim.x) + threadIdx.x + 1;
+	col = (col * 2) - 1;
+
+	int NUM_2 = NUM >> 1;
+
+	// p_i,0 = p_i,1
+	pres_black(col, 0) = pres_red(col, 1);
+	pres_red(col + 1, 0) = pres_black(col + 1, 1);
+
+	// p_i,jmax+1 = p_i,jmax
+	pres_red(col, NUM_2 + 1) = pres_black(col, NUM_2);
+	pres_black(col + 1, NUM_2 + 1) = pres_red(col + 1, NUM_2);
+
+} // end set_horz_pres_BCs
+
+//////////////////////////////////////////////////////////////////////////////
+
+__global__
+void set_vert_pres_BCs (Real* pres_red, Real* pres_black) 
+{
+	int row = (blockIdx.x * blockDim.x) + threadIdx.x + 1;
+
+	int NUM_2 = NUM >> 1;
+
+	// p_0,j = p_1,j
+	pres_black(0, row) = pres_red(1, row);
+	pres_red(0, row) = pres_black(1, row);
+
+	// p_imax+1,j = p_imax,j
+	pres_black(NUM + 1, row) = pres_red(NUM, row);
+	pres_red(NUM + 1, row) = pres_black(NUM, row);
+
+} // end set_pressure_BCs
+
+///////////////////////////////////////////////////////////////////////////////
+
 /** Function to update pressure for red cells
  * 
- * \param[in]			dt					time-step size
- * \param[in]			F						array of discretized x-momentum eqn terms
- * \param[in]			G						array of discretized y-momentum eqn terms
- * \param[in]			pres_black	pressure values of black cells
- * \param[inout]	pres_red		pressure values of red cells
- * \param[inout]	norm_L2			variable holding summed residuals
+ * \param[in]		dt			time-step size
+ * \param[in]		F			array of discretized x-momentum eqn terms
+ * \param[in]		G			array of discretized y-momentum eqn terms
+ * \param[in]		pres_black	pressure values of black cells
+ * \param[inout]	pres_red	pressure values of red cells
+ */
+__global__
+void red_kernel (const Real dt, const Real* F, 
+				 const Real* G, const Real* pres_black,
+				 Real* pres_red) 
+{
+	int row = (blockIdx.x * blockDim.x) + threadIdx.x + 1;
+	int col = (blockIdx.y * blockDim.y) + threadIdx.y + 1;
+	
+	int NUM_2 = NUM >> 1;			
+	
+	Real p_ij = pres_red(col, row);
+	
+	Real p_im1j = pres_black(col - 1, row);
+	Real p_ip1j = pres_black(col + 1, row);
+	Real p_ijm1 = pres_black(col, row - (col & 1));
+	Real p_ijp1 = pres_black(col, row + ((col + 1) & 1));
+	
+	// right-hand side
+	Real rhs = (((F(col, (2 * row) - (col & 1))
+			    - F(col - 1, (2 * row) - (col & 1))) / dx)
+			  + ((G(col, (2 * row) - (col & 1))
+			    - G(col, (2 * row) - (col & 1) - 1)) / dy)) / dt;
+	
+	pres_red(col, row) = p_ij * (ONE - omega) + omega * 
+		(((p_ip1j + p_im1j) / (dx * dx)) + ((p_ijp1 + p_ijm1) / (dy * dy)) - 
+		rhs) / ((TWO / (dx * dx)) + (TWO / (dy * dy)));
+	
+} // end red_kernel
+
+///////////////////////////////////////////////////////////////////////////////
+
+/** Function to update pressure for black cells
+ * 
+ * \param[in]		dt			time-step size
+ * \param[in]		F			array of discretized x-momentum eqn terms
+ * \param[in]		G			array of discretized y-momentum eqn terms
+ * \param[in]		pres_red	pressure values of red cells
+ * \param[inout]	pres_black	pressure values of black cells
  */
 __global__ 
-void red_kernel (const Real dt, const Real * F, const Real * G, const Real * pres_black,
-								 Real * pres_red, Real * norm_L2)
-{
-	int row = (blockIdx.y * blockDim.y) + threadIdx.y;
-	int col = (blockIdx.x * blockDim.x) + threadIdx.x;
-	
-	int ind_red = (col * (NUM >> 1)) + row;  					// local (red) index
-	int ind = (col * NUM) + (2 * row) + (col & 1);		// global index
-
-	Real p_ij = pres_red[ind_red];			
-
-	// left boundary
-	Real p_im1j;
-	Real F_im1j;
-	if (col == 0) {
-		p_im1j = p_ij;
-		F_im1j = ZERO;
+void black_kernel (const Real dt, const Real* F, 
+				   const Real* G, const Real* pres_red, 
+				   Real* pres_black) 
+{
+	int row = (blockIdx.x * blockDim.x) + threadIdx.x + 1;
+	int col = (blockIdx.y * blockDim.y) + threadIdx.y + 1;
+	
+	int NUM_2 = NUM >> 1;
+	
+	Real p_ij = pres_black(col, row);
+
+	Real p_im1j = pres_red(col - 1, row);
+	Real p_ip1j = pres_red(col + 1, row);
+	Real p_ijm1 = pres_red(col, row - ((col + 1) & 1));
+	Real p_ijp1 = pres_red(col, row + (col & 1));
+	
+	// right-hand side
+	Real rhs = (((F(col, (2 * row) - ((col + 1) & 1))
+			 	- F(col - 1, (2 * row) - ((col + 1) & 1))) / dx)
+			  + ((G(col, (2 * row) - ((col + 1) & 1))
+			    - G(col, (2 * row) - ((col + 1) & 1) - 1)) / dy)) / dt;
+	
+	pres_black(col, row) = p_ij * (ONE - omega) + omega * 
+		(((p_ip1j + p_im1j) / (dx * dx)) + ((p_ijp1 + p_ijm1) / (dy * dy)) - 
+		rhs) / ((TWO / (dx * dx)) + (TWO / (dy * dy)));
+	
+} // end black_kernel
+
+///////////////////////////////////////////////////////////////////////////////
+
+__global__
+void calc_residual (const Real dt, const Real* F, const Real* G, 
+					const Real* pres_red, const Real* pres_black,
+					Real* res_array)
+{
+	int row = (blockIdx.x * blockDim.x) + threadIdx.x + 1;
+	int col = (blockIdx.y * blockDim.y) + threadIdx.y + 1;
+	
+	int NUM_2 = NUM >> 1;
+
+	Real p_ij, p_im1j, p_ip1j, p_ijm1, p_ijp1, rhs, res, res2;
+
+	// red point
+	p_ij = pres_red(col, row);
+	
+	p_im1j = pres_black(col - 1, row);
+	p_ip1j = pres_black(col + 1, row);
+	p_ijm1 = pres_black(col, row - (col & 1));
+	p_ijp1 = pres_black(col, row + ((col + 1) & 1));
+
+	rhs = (((F(col, (2 * row) - (col & 1)) - F(col - 1, (2 * row) - (col & 1))) / dx)
+		+  ((G(col, (2 * row) - (col & 1)) - G(col, (2 * row) - (col & 1) - 1)) / dy)) / dt;
+
+	// calculate residual
+	res = ((p_ip1j - (TWO * p_ij) + p_im1j) / (dx * dx))
+		+ ((p_ijp1 - (TWO * p_ij) + p_ijm1) / (dy * dy)) - rhs;
+
+	// black point
+	p_ij = pres_black(col, row);
+
+	p_im1j = pres_red(col - 1, row);
+	p_ip1j = pres_red(col + 1, row);
+	p_ijm1 = pres_red(col, row - ((col + 1) & 1));
+	p_ijp1 = pres_red(col, row + (col & 1));
+	
+	// right-hand side
+	rhs = (((F(col, (2 * row) - ((col + 1) & 1)) - F(col - 1, (2 * row) - ((col + 1) & 1))) / dx)
+		+  ((G(col, (2 * row) - ((col + 1) & 1)) - G(col, (2 * row) - ((col + 1) & 1) - 1)) / dy)) / dt;
+
+	// calculate residual
+	res2 = ((p_ip1j - (TWO * p_ij) + p_im1j) / (dx * dx))
+		 + ((p_ijp1 - (TWO * p_ij) + p_ijm1) / (dy * dy)) - rhs;
+
+	// shared memory for block's sum
+	__shared__ Real sum_cache[BLOCK_SIZE];
+
+	sum_cache[threadIdx.x] = (res * res) + (res2 * res2);
+
+	// synchronize threads in block to ensure all residuals stored
+	__syncthreads();
+
+	// add up squared residuals for block
+	int i = BLOCK_SIZE >> 1;
+	while (i != 0) {
+		if (threadIdx.x < i) {
+			sum_cache[threadIdx.x] += sum_cache[threadIdx.x + i];
+		}
+		__syncthreads();
+		i >>= 1;
+	}
+
+	// store block's summed residuals
+	if (threadIdx.x == 0) {
+		res_array[blockIdx.y + (gridDim.y * blockIdx.x)] = sum_cache[0];
+	}
+} 
+
+///////////////////////////////////////////////////////////////////////////////
+
+__global__ 
+void calculate_u (const Real dt, const Real* F, 
+				  const Real* pres_red, const Real* pres_black, 
+				  Real* u, Real* max_u)
+{
+	int row = (blockIdx.x * blockDim.x) + threadIdx.x + 1;
+	int col = (blockIdx.y * blockDim.y) + threadIdx.y + 1;
+	
+	// allocate shared memory to store max velocities
+	__shared__ Real max_cache[BLOCK_SIZE];
+	max_cache[threadIdx.x] = ZERO;
+	
+	int NUM_2 = NUM >> 1;
+	Real new_u = ZERO;
+
+	if (col != NUM) {
+
+		Real p_ij, p_ip1j, new_u2;
+
+		// red point
+		p_ij = pres_red(col, row);
+		p_ip1j = pres_black(col + 1, row);
+
+		new_u = F(col, (2 * row) - (col & 1)) - (dt * (p_ip1j - p_ij) / dx);
+		u(col, (2 * row) - (col & 1)) = new_u;
+
+		// black point
+		p_ij = pres_black(col, row);
+		p_ip1j = pres_red(col + 1, row);
+
+		new_u2 = F(col, (2 * row) - ((col + 1) & 1)) - (dt * (p_ip1j - p_ij) / dx);
+		u(col, (2 * row) - ((col + 1) & 1)) = new_u2;
+
+		// check for max of these two
+		new_u = fmax(fabs(new_u), fabs(new_u2));
+
+		if ((2 * row) == NUM) {
+			// also test for max velocity at vertical boundary
+			new_u = fmax(new_u, fabs( u(col, NUM + 1) ));
+		}
 	} else {
-		p_im1j = pres_black[((col - 1) * (NUM >> 1)) + row];
-		F_im1j = F[((col - 1) * NUM) + (2 * row) + (col & 1)];
-	}
-
-	// right boundary
-	Real p_ip1j;
-	if (col == (NUM - 1)) {
-		p_ip1j = p_ij;
+		// check for maximum velocity in boundary cells also
+		new_u = fmax(fabs( u(NUM, (2 * row)) ), fabs( u(0, (2 * row)) ));
+		new_u = fmax(fabs( u(NUM, (2 * row) - 1) ), new_u);
+		new_u = fmax(fabs( u(0, (2 * row) - 1) ), new_u);
+
+		new_u = fmax(fabs( u(NUM + 1, (2 * row)) ), new_u);
+		new_u = fmax(fabs( u(NUM + 1, (2 * row) - 1) ), new_u);
+
+	} // end if
+
+	// store maximum u for block from each thread
+	max_cache[threadIdx.x] = new_u;
+
+	// synchronize threads in block to ensure all velocities stored
+	__syncthreads();
+
+	// calculate maximum for block
+	int i = BLOCK_SIZE >> 1;
+	while (i != 0) {
+		if (threadIdx.x < i) {
+			max_cache[threadIdx.x] = fmax(max_cache[threadIdx.x], max_cache[threadIdx.x + i]);
+		}
+		__syncthreads();
+		i >>= 1;
+	}
+
+	// store block's maximum
+	if (threadIdx.x == 0) {
+		max_u[blockIdx.y + (gridDim.y * blockIdx.x)] = max_cache[0];
+	}
+
+	
+} // end calculate_u
+
+///////////////////////////////////////////////////////////////////////////////
+
+__global__ 
+void calculate_v (const Real dt, const Real* G, 
+				  const Real* pres_red, const Real* pres_black, 
+				  Real* v, Real* max_v)
+{
+	int row = (blockIdx.x * blockDim.x) + threadIdx.x + 1;
+	int col = (blockIdx.y * blockDim.y) + threadIdx.y + 1;
+	
+	// allocate shared memory to store maximum velocities
+	__shared__ Real max_cache[BLOCK_SIZE];
+	max_cache[threadIdx.x] = ZERO;
+
+	int NUM_2 = NUM >> 1;
+	Real new_v = ZERO;
+	
+	if (row != NUM_2) {
+		Real p_ij, p_ijp1, new_v2;
+
+		// red pressure point
+		p_ij = pres_red(col, row);
+		p_ijp1 = pres_black(col, row + ((col + 1) & 1));
+	
+		new_v = G(col, (2 * row) - (col & 1)) - (dt * (p_ijp1 - p_ij) / dy);
+		v(col, (2 * row) - (col & 1)) = new_v;
+
+
+		// black pressure point
+		p_ij = pres_black(col, row);
+		p_ijp1 = pres_red(col, row + (col & 1));
+		
+		new_v2 = G(col, (2 * row) - ((col + 1) & 1)) - (dt * (p_ijp1 - p_ij) / dy);
+		v(col, (2 * row) - ((col + 1) & 1)) = new_v2;
+
+
+		// check for max of these two
+		new_v = fmax(fabs(new_v), fabs(new_v2));
+
+		if (col == NUM) {
+			// also test for max velocity at vertical boundary
+			new_v = fmax(new_v, fabs( v(NUM + 1, (2 * row)) ));
+		}
+
 	} else {
-		p_ip1j = pres_black[((col + 1) * (NUM >> 1)) + row];
-	}
-
-	// bottom boundary
-	Real p_ijm1;
-	Real G_ijm1;
-	if (((2 * row) + (col & 1)) == 0) {
-		p_ijm1 = p_ij;
-		G_ijm1 = ZERO;
-	} else {
-		p_ijm1 = pres_black[(col * (NUM >> 1)) + row - ((col + 1) & 1)];
-		G_ijm1 = G[(col * NUM) + (2 * row) + (col & 1) - 1];
-	}
-
-	// top boundary
-	Real p_ijp1;
-	if (((2 * row) + (col & 1)) == (NUM - 1)) {
-		p_ijp1 = p_ij;
-	} else {
-		p_ijp1 = pres_black[(col * (NUM >> 1)) + row + (col & 1)];
-	}
-
-	// right-hand side
-	Real rhs = (((F[ind] - F_im1j) / dx) + ((G[ind] - G_ijm1) / dy)) / dt;
-
-	pres_red[ind_red] = p_ij * (ONE - omega) + omega * (
-										  ((p_ip1j + p_im1j) / (dx * dx)) + ((p_ijp1 + p_ijm1) / (dy * dy)) - 
-										  rhs) / ((TWO / (dx * dx)) + (TWO / (dy * dy)));
-	
-	// calculate residual (reuse rhs variable)
-	rhs = ((p_ip1j - (TWO * p_ij) + p_im1j) / (dx * dx)) + ((p_ijp1 - (TWO * p_ij) + p_ijm1) / (dy * dy)) - rhs;
-	
-	#ifdef SHARED
-		// store residual for block
-		__shared__ Real res_cache[BLOCK_SIZE];
-	
-		// store squared residual from each thread
-		res_cache[threadIdx.y] = rhs * rhs;
-	
-		// synchronize threads in block
+
+		if ((col & 1) == 1) {
+			// black point is on boundary, only calculate red point below it
+			Real p_ij = pres_red(col, row);
+			Real p_ijp1 = pres_black(col, row + ((col + 1) & 1));
+		
+			new_v = G(col, (2 * row) - (col & 1)) - (dt * (p_ijp1 - p_ij) / dy);
+			v(col, (2 * row) - (col & 1)) = new_v;
+
+		} else {
+			// red point is on boundary, only calculate black point below it
+			Real p_ij = pres_black(col, row);
+			Real p_ijp1 = pres_red(col, row + (col & 1));
+		
+			new_v = G(col, (2 * row) - ((col + 1) & 1)) - (dt * (p_ijp1 - p_ij) / dy);
+			v(col, (2 * row) - ((col + 1) & 1)) = new_v;
+		}
+
+		// get maximum v velocity
+		new_v = fabs(new_v);
+
+		// check for maximum velocity in boundary cells also
+		new_v = fmax(fabs( v(col, NUM) ), new_v);
+		new_v = fmax(fabs( v(col, 0) ), new_v);
+
+		new_v = fmax(fabs( v(col, NUM + 1) ), new_v);
+
+	} // end if
+		
+	// store absolute value of velocity
+	max_cache[threadIdx.x] = new_v;
+	
+	// synchronize threads in block to ensure all velocities stored
+	__syncthreads();
+
+	// calculate maximum for block
+	int i = BLOCK_SIZE >> 1;
+	while (i != 0) {
+		if (threadIdx.x < i) {
+			max_cache[threadIdx.x] = fmax(max_cache[threadIdx.x], max_cache[threadIdx.x + i]);
+		}
 		__syncthreads();
-	
-		// add up squared residuals for block
-		int i = BLOCK_SIZE >> 1;
-		while (i != 0) {
-			if (threadIdx.y < i) {
-				res_cache[threadIdx.y] += res_cache[threadIdx.y + i];
-			}
-			__syncthreads();
-			i >>= 1;
-		}
-	
-		// store block's summed residuals
-		if (threadIdx.y == 0) {
-			#ifdef ATOMIC
-				atomicAdd (norm_L2, res_cache[0]);
-			#else
-				norm_L2[blockIdx.x + (gridDim.x * blockIdx.y)] = res_cache[0];
-			#endif
-		}
-	#else
-		norm_L2[ind_red] = rhs * rhs;
-	#endif
-	
-} // end red_kernel
-
-///////////////////////////////////////////////////////////////////////////////
-
-/** Function to update pressure for black cells
- * 
- * \param[in]			dt					time-step size
- * \param[in]			F						array of discretized x-momentum eqn terms
- * \param[in]			G						array of discretized y-momentum eqn terms
- * \param[in]			pres_red		pressure values of red cells
- * \param[inout]	pres_black	pressure values of black cells
- * \param[inout]	norm_L2			variable holding summed residuals
- */
-__global__ 
-void black_kernel (const Real dt, const Real * F, const Real * G, const Real * pres_red,
-									 Real * pres_black, Real * norm_L2)
-{
-	int row = (blockIdx.y * blockDim.y) + threadIdx.y;
-	int col = (blockIdx.x * blockDim.x) + threadIdx.x;
-	
-	int ind_black = (col * (NUM >> 1)) + row;  						// local (black) index
-	int ind = (col * NUM) + (2 * row) + ((col + 1) & 1);	// global index
-	
-	Real p_ij = pres_black[ind_black];
-	
-	// left boundary
-	Real p_im1j;
-	Real F_im1j;
-	if (col == 0) {
-		p_im1j = p_ij;
-		F_im1j = ZERO;
-	} else {
-		p_im1j = pres_red[((col - 1) * (NUM >> 1)) + row];
-		F_im1j = F[((col - 1) * NUM) + (2 * row) + ((col + 1) & 1)];
-	}
-	
-	// right boundary
-	Real p_ip1j;
-	if (col == (NUM - 1)) {
-		p_ip1j = p_ij;
-	} else {
-		p_ip1j = pres_red[((col + 1) * (NUM >> 1)) + row];
-	}
-	
-	// bottom boundary
-	Real p_ijm1;
-	Real G_ijm1;
-	if (((2 * row) + ((col + 1) & 1)) == 0) {
-		p_ijm1 = p_ij;
-		G_ijm1 = ZERO;
-	} else {
-		p_ijm1 = pres_red[(col * (NUM >> 1)) + row - (col & 1)];
-		G_ijm1 = G[(col * NUM) + (2 * row) + ((col + 1) & 1) - 1];
-	}
-	
-	// top boundary
-	Real p_ijp1;
-	if (((2 * row) + ((col + 1) & 1)) == (NUM - 1)) {
-		p_ijp1 = p_ij;
-	} else {
-		p_ijp1 = pres_red[(col * (NUM >> 1)) + row + ((col + 1) & 1)];
-	}
-	
-	// right-hand side
-	Real rhs = (((F[ind] - F_im1j) / dx) + ((G[ind] - G_ijm1) / dy)) / dt;
-	
-	pres_black[ind_black] = p_ij * (ONE - omega) + omega * (
-										  		((p_ip1j + p_im1j) / (dx * dx)) + ((p_ijp1 + p_ijm1) / (dy * dy)) - 
-										  		rhs) / ((TWO / (dx * dx)) + (TWO / (dy * dy)));
-	
-	// calculate residual (reuse rhs variable)
-	rhs = ((p_ip1j - (TWO * p_ij) + p_im1j) / (dx * dx)) + ((p_ijp1 - (TWO * p_ij) + p_ijm1) / (dy * dy)) - rhs;
-	
-	#ifdef SHARED
-		// store residual for block
-		__shared__ Real res_cache[BLOCK_SIZE];
-	
-		// store squared residual from each thread
-		res_cache[threadIdx.y] = rhs * rhs;
-	
-		// synchronize threads in block
-		__syncthreads();
-	
-		// add up squared residuals for block
-		int i = BLOCK_SIZE >> 1;
-		while (i != 0) {
-			if (threadIdx.y < i) {
-				res_cache[threadIdx.y] += res_cache[threadIdx.y + i];
-			}
-			__syncthreads();
-			i >>= 1;
-		}
-	
-		// store block's summed residuals
-		if (threadIdx.y == 0) {
-			#ifdef ATOMIC
-				atomicAdd (norm_L2, res_cache[0]);
-			#else
-				norm_L2[blockIdx.x + (gridDim.x * blockIdx.y)] = res_cache[0];
-			#endif
-		}
-	#else
-		norm_L2[ind_black] = rhs * rhs;
-	#endif
-	
-} // end black_kernel
-
-///////////////////////////////////////////////////////////////////////////////
-
-__global__ 
-void calculate_u (const Real dt, const Real * F, 
-									const Real * pres_red, const Real * pres_black, 
-									#ifdef ATOMIC
-									Real * u, unsigned int * max_u_d)
-									#else
-									Real * u, Real * max_u_d)
-									#endif
-{
-	int row = (blockIdx.y * blockDim.y) + threadIdx.y;
-	int col = (blockIdx.x * blockDim.x) + threadIdx.x;
-	
-	#ifdef SHARED
-		#ifdef ATOMIC
-			__shared__ unsigned int res_cache[BLOCK_SIZE];
-			res_cache[threadIdx.y] = 0;
-		#else
-			__shared__ Real res_cache[BLOCK_SIZE];
-			res_cache[threadIdx.y] = ZERO;
-		#endif
-	#endif
-	
-	if (col != (NUM - 1)) {
-		int ind = (col * NUM) + row;
-		
-		Real p_ij, p_ip1j;
-		if (((row + col) & 1) == 0) {
-			// red pressure cell
-			p_ij = pres_red[(col * (NUM >> 1)) + ((row - (col & 1)) >> 1)];
-			
-			// p_ip1j is black cell
-			p_ip1j = pres_black[((col + 1) * (NUM >> 1)) + ((row - (col & 1)) >> 1)];
-		} else {
-			// black pressure cell
-			p_ij = pres_black[(col * (NUM >> 1)) + ((row - ((col + 1) & 1)) >> 1)];
-			
-			// p_ip1j is red cell
-			p_ip1j = pres_red[((col + 1) * (NUM >> 1)) + ((row - ((col + 1) & 1)) >> 1)];
-		}
-		
-		//u[ind] = F[ind] - (dt * (p_ip1j - p_ij) / dx);
-		Real u_ij = F[ind] - (dt * (p_ip1j - p_ij) / dx);
-		
-		u[ind] = u_ij;
-		
-		#ifdef SHARED
-		// store maximum u for block from each thread
-			#ifdef ATOMIC
-				res_cache[threadIdx.y] = floatFlip (fabs(u_ij));
-			#else
-				res_cache[threadIdx.y] = fabs(u_ij);
-			#endif
-		
-		// synchronize threads in block
-		__syncthreads();
-
-		// add up squared residuals for block
-		int i = BLOCK_SIZE >> 1;
-		while (i != 0) {
-			if (threadIdx.y < i) {
-				res_cache[threadIdx.y] = fmax(res_cache[threadIdx.y], res_cache[threadIdx.y + i]);
-			}
-			__syncthreads();
-			i >>= 1;
-		}
-
-		// store block's summed residuals
-		if (threadIdx.y == 0) {
-			#ifdef ATOMIC
-				atomicMax (max_u_d, res_cache[0]);
-			#else
-				max_u_d[blockIdx.x + (gridDim.x * blockIdx.y)] = res_cache[0];
-			#endif
-		}
-		#endif
-	} // end if
-	
-} // end calculate_u
-
-///////////////////////////////////////////////////////////////////////////////
-
-__global__ 
-void calculate_v (const Real dt, const Real * G, 
-									const Real * pres_red, const Real * pres_black, 
-									#ifdef ATOMIC
-									Real * v, unsigned int * max_v_d)
-									#else
-									Real * v, Real * max_v_d)
-									#endif
-{
-	int row = (blockIdx.y * blockDim.y) + threadIdx.y;
-	int col = (blockIdx.x * blockDim.x) + threadIdx.x;
-	
-	#ifdef SHARED
-		#ifdef ATOMIC
-			__shared__ unsigned int res_cache[BLOCK_SIZE];
-			res_cache[threadIdx.y] = 0;
-		#else
-			__shared__ Real res_cache[BLOCK_SIZE];
-			res_cache[threadIdx.y] = ZERO;
-		#endif
-	#endif
-	
-	if (row != (NUM - 1)) {
-		int ind = (col * NUM) + row;
-		
-		Real p_ij, p_ijp1;
-		if (((row + col) & 1) == 0) {
-			// red pressure cell
-			p_ij = pres_red[(col * (NUM >> 1)) + ((row - (col & 1)) >> 1)];
-			
-			// p_ijp1 is black cell
-			p_ijp1 = pres_black[(col * (NUM >> 1)) + ((row + 1 - ((col + 1) & 1)) >> 1)];
-		} else {
-			// black pressure cell
-			p_ij = pres_black[(col * (NUM >> 1)) + ((row - ((col + 1) & 1)) >> 1)];
-			
-			// p_ijp1 is red cell
-			p_ijp1 = pres_red[(col * (NUM >> 1)) + ((row + 1 - (col & 1)) >> 1)];
-		}
-		
-		//v[ind] = G[ind] - (dt * (p_ijp1 - p_ij) / dy);
-		Real v_ij = G[ind] - (dt * (p_ijp1 - p_ij) / dy);
-		
-		v[ind] = v_ij;
-		
-		#ifdef SHARED
-		// store maximum v for block for each thread
-		
-			#ifdef ATOMIC
-				res_cache[threadIdx.y] = floatFlip (fabs(v_ij));
-			#else
-				res_cache[threadIdx.y] = fabs(v_ij);
-			#endif
-		
-		// synchronize threads in block
-		__syncthreads();
-
-		// add up squared residuals for block
-		int i = BLOCK_SIZE >> 1;
-		while (i != 0) {
-			if (threadIdx.y < i) {
-				res_cache[threadIdx.y] = fmax(res_cache[threadIdx.y], res_cache[threadIdx.y + i]);
-			}
-			__syncthreads();
-			i >>= 1;
-		}
-
-		// store block's summed residuals
-		if (threadIdx.y == 0) {
-			#ifdef ATOMIC
-				atomicMax (max_v_d, res_cache[0]);
-			#else
-				max_v_d[blockIdx.x + (gridDim.x * blockIdx.y)] = res_cache[0];
-			#endif
-		}
-		#endif
-	} // end if
+		i >>= 1;
+	}
+
+	// store block's summed residuals
+	if (threadIdx.x == 0) {
+		max_v[blockIdx.y + (gridDim.y * blockIdx.x)] = max_cache[0];
+	}
 	
 } // end calculate_v
 

	}
	
	fclose(pfile);
	
	free(pres_red);
	free(pres_black);

int main (void)
{
	// iterations for Red-Black Gauss-Seidel with SOR
	int iter = 0;
	const int it_max = 10000;
	
	// SOR iteration tolerance
	const Real tol = 0.001;
