diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index 06afeb4d..a8dd8dd8 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -62,7 +62,7 @@ extern "C" void ScaLBL_D3Q19_Pressure_BC_Z(double *disteven, double *distodd, do
 extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID,  double *disteven, double *distodd, double flux,
 								  int Nx, int Ny, int Nz);
 
-extern "C" double ScaLBL_D3Q19_Flux_BC_Z(double *disteven, double *distodd, double flux,
+extern "C" double ScaLBL_D3Q19_Flux_BC_Z(char *ID, double *disteven, double *distodd, double flux,
 								   int Nx, int Ny, int Nz, int outlet);
 
 extern "C" void ScaLBL_Color_Init(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, int Nz);
diff --git a/cpu/D3Q19.cpp b/cpu/D3Q19.cpp
index 52473181..9428913f 100644
--- a/cpu/D3Q19.cpp
+++ b/cpu/D3Q19.cpp
@@ -293,7 +293,7 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID,  double *disteven, double *di
 	double din = 0.f;
 	N = Nx*Ny*Nz;
 
-	double A = 1.f*double(Nx*Ny);
+	double A = 1.f*double((Nx-2)*(Ny-2));
 	double sum = 0.f;
 	char id;
 	for (n=Nx*Ny; n<2*Nx*Ny; n++){
@@ -303,26 +303,26 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID,  double *disteven, double *di
 			// Read distributions from "opposite" memory convention
 			//........................................................................
 			//........................................................................
-			f1 = distodd[n];
-			f3 = distodd[N+n];
-			f5 = distodd[2*N+n];
-			f7 = distodd[3*N+n];
-			f9 = distodd[4*N+n];
-			f11 = distodd[5*N+n];
-			f13 = distodd[6*N+n];
-			f15 = distodd[7*N+n];
-			f17 = distodd[8*N+n];
+			f2 = distodd[n];
+			f4 = distodd[N+n];
+			f6 = distodd[2*N+n];
+			f8 = distodd[3*N+n];
+			f10 = distodd[4*N+n];
+			f12 = distodd[5*N+n];
+			//f14 = distodd[6*N+n];
+			f16 = distodd[7*N+n];
+			//f18 = distodd[8*N+n];
 			//........................................................................
 			f0 = disteven[n];
-			f2 = disteven[N+n];
-			f4 = disteven[2*N+n];
-			f6 = disteven[3*N+n];
-			f8 = disteven[4*N+n];
-			f10 = disteven[5*N+n];
-			f12 = disteven[6*N+n];
-			f14 = disteven[7*N+n];
-			f16 = disteven[8*N+n];
-			f18 = disteven[9*N+n];
+			f1 = disteven[N+n];
+			f3 = disteven[2*N+n];
+			//f5 = disteven[3*N+n];
+			f7 = disteven[4*N+n];
+			f9 = disteven[5*N+n];
+			//f11 = disteven[6*N+n];
+			f13 = disteven[7*N+n];
+			//f15 = disteven[8*N+n];
+			f17 = disteven[9*N+n];
 			//...................................................
 
 			// Determine the outlet flow velocity
@@ -335,7 +335,7 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_z(char *ID,  double *disteven, double *di
 	return din;
 }
 
-extern "C" double ScaLBL_D3Q19_Flux_BC_Z(double *disteven, double *distodd, double flux,
+extern "C" double ScaLBL_D3Q19_Flux_BC_Z(char *ID, double *disteven, double *distodd, double flux,
 		int Nx, int Ny, int Nz, int outlet){
 	// Note that this routine assumes the distributions are stored "opposite"
 	// odd distributions in disteven and even distributions in distodd.
@@ -344,40 +344,41 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_Z(double *disteven, double *distodd, doub
 	double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
 	double f10,f11,f12,f13,f14,f15,f16,f17,f18;
 	double dout = 0.f;
-
 	N = Nx*Ny*Nz;
 
 	// Loop over the boundary - threadblocks delineated by start...finish
-	double A = 1.f*double(Nx*Ny);
+	double A = 1.f*double((Nx-2)*(Ny-2));
 	double sum = 0.f;
+    char id;
 	for (n=outlet; n<N-Nx*Ny; n++){
+        id = ID[n];
+        if (id>0){
+            //........................................................................
+            // Read distributions from "opposite" memory convention
+            //........................................................................
+            f2 = distodd[n];
+            f4 = distodd[N+n];
+            //f6 = distodd[2*N+n];
+            f8 = distodd[3*N+n];
+            f10 = distodd[4*N+n];
+            //f12 = distodd[5*N+n];
+            f14 = distodd[6*N+n];
+            //f16 = distodd[7*N+n];
+            f18 = distodd[8*N+n];
+            //........................................................................
+            f0 = disteven[n];
+            f1 = disteven[N+n];
+            f3 = disteven[2*N+n];
+            f5 = disteven[3*N+n];
+            f7 = disteven[4*N+n];
+            f9 = disteven[5*N+n];
+            f11 = disteven[6*N+n];
+            //f13 = disteven[7*N+n];
+            f15 = disteven[8*N+n];
+            //f17 = disteven[9*N+n];
 
-		//........................................................................
-		// Read distributions from "opposite" memory convention
-		//........................................................................
-		f1 = distodd[n];
-		f3 = distodd[N+n];
-		f5 = distodd[2*N+n];
-		f7 = distodd[3*N+n];
-		f9 = distodd[4*N+n];
-		f11 = distodd[5*N+n];
-		f13 = distodd[6*N+n];
-		f15 = distodd[7*N+n];
-		f17 = distodd[8*N+n];
-		//........................................................................
-		f0 = disteven[n];
-		f2 = disteven[N+n];
-		f4 = disteven[2*N+n];
-		f6 = disteven[3*N+n];
-		f8 = disteven[4*N+n];
-		f10 = disteven[5*N+n];
-		f12 = disteven[6*N+n];
-		f14 = disteven[7*N+n];
-		f16 = disteven[8*N+n];
-		f18 = disteven[9*N+n];
-
-		sum += (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f5+f11+f14+f15+f18));
-
+            sum += (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f5+f11+f14+f15+f18));
+        }
 	}
 	dout = sum/(A*(1.0+flux));
 	return dout;
@@ -402,26 +403,26 @@ extern "C" void ScaLBL_D3Q19_Pressure_BC_z(double *disteven, double *distodd, do
 		// Read distributions from "opposite" memory convention
 		//........................................................................
 		//........................................................................
-		f1 = distodd[n];
-		f3 = distodd[N+n];
-		f5 = distodd[2*N+n];
-		f7 = distodd[3*N+n];
-		f9 = distodd[4*N+n];
-		f11 = distodd[5*N+n];
-		f13 = distodd[6*N+n];
-		f15 = distodd[7*N+n];
-		f17 = distodd[8*N+n];
+		f2 = distodd[n];
+		f4 = distodd[N+n];
+		f6 = distodd[2*N+n];
+		f8 = distodd[3*N+n];
+		f10 = distodd[4*N+n];
+		f12 = distodd[5*N+n];
+		f14 = distodd[6*N+n];
+		f16 = distodd[7*N+n];
+		f18 = distodd[8*N+n];
 		//........................................................................
 		f0 = disteven[n];
-		f2 = disteven[N+n];
-		f4 = disteven[2*N+n];
-		f6 = disteven[3*N+n];
-		f8 = disteven[4*N+n];
-		f10 = disteven[5*N+n];
-		f12 = disteven[6*N+n];
-		f14 = disteven[7*N+n];
-		f16 = disteven[8*N+n];
-		f18 = disteven[9*N+n];
+		f1 = disteven[N+n];
+		f3 = disteven[2*N+n];
+		f5 = disteven[3*N+n];
+		f7 = disteven[4*N+n];
+		f9 = disteven[5*N+n];
+		f11 = disteven[6*N+n];
+		f13 = disteven[7*N+n];
+		f15 = disteven[8*N+n];
+		f17 = disteven[9*N+n];
 		//...................................................
 		//........Determine the inlet flow velocity.........
 		//			uz = -1 + (f0+f3+f4+f1+f2+f7+f8+f10+f9
@@ -447,11 +448,19 @@ extern "C" void ScaLBL_D3Q19_Pressure_BC_z(double *disteven, double *distodd, do
 		f15 = f16 + 0.16666666666666678*(uy+uz)-Cyz;
 		f18 = f17 + 0.16666666666666678*(uz-uy)+Cyz;
 		//........Store in "opposite" memory location..........
-		distodd[2*N+n] = f5;
+		/*		distodd[2*N+n] = f5;
 		distodd[5*N+n] = f11;
 		disteven[7*N+n] = f14;
 		distodd[7*N+n] = f15;
 		disteven[9*N+n] = f18;
+
+		*/
+
+		disteven[3*N+n] = f5;
+		disteven[6*N+n] = f11;
+		distodd[6*N+n] = f14;
+		disteven[8*N+n] = f15;
+		distodd[8*N+n] = f18;
 		/*
 		printf("Site=%i\n",n);
 		printf("ux=%f, uy=%f, uz=%f\n",ux,uy,uz);
@@ -480,26 +489,26 @@ extern "C" void ScaLBL_D3Q19_Pressure_BC_Z(double *disteven, double *distodd, do
 		//........................................................................
 		// Read distributions from "opposite" memory convention
 		//........................................................................
-		f1 = distodd[n];
-		f3 = distodd[N+n];
-		f5 = distodd[2*N+n];
-		f7 = distodd[3*N+n];
-		f9 = distodd[4*N+n];
-		f11 = distodd[5*N+n];
-		f13 = distodd[6*N+n];
-		f15 = distodd[7*N+n];
-		f17 = distodd[8*N+n];
+		f2 = distodd[n];
+		f4 = distodd[N+n];
+		f6 = distodd[2*N+n];
+		f8 = distodd[3*N+n];
+		f10 = distodd[4*N+n];
+		f12 = distodd[5*N+n];
+		f14 = distodd[6*N+n];
+		f16 = distodd[7*N+n];
+		f18 = distodd[8*N+n];
 		//........................................................................
 		f0 = disteven[n];
-		f2 = disteven[N+n];
-		f4 = disteven[2*N+n];
-		f6 = disteven[3*N+n];
-		f8 = disteven[4*N+n];
-		f10 = disteven[5*N+n];
-		f12 = disteven[6*N+n];
-		f14 = disteven[7*N+n];
-		f16 = disteven[8*N+n];
-		f18 = disteven[9*N+n];
+		f1 = disteven[N+n];
+		f3 = disteven[2*N+n];
+		f5 = disteven[3*N+n];
+		f7 = disteven[4*N+n];
+		f9 = disteven[5*N+n];
+		f11 = disteven[6*N+n];
+		f13 = disteven[7*N+n];
+		f15 = disteven[8*N+n];
+		f17 = disteven[9*N+n];
 		//........Determine the outlet flow velocity.........
 		//			uz = 1 - (f0+f3+f4+f1+f2+f7+f8+f10+f9+
 		//					  2*(f6+f16+f17+f12+f13))/dout;
@@ -531,11 +540,11 @@ extern "C" void ScaLBL_D3Q19_Pressure_BC_Z(double *disteven, double *distodd, do
 		f17 = f18 - 0.16666666666666678*(uz-uy)-Cyz;
 
 		//........Store in "opposite" memory location..........
-		disteven[3*N+n] = f6;
-		disteven[6*N+n] = f12;
-		distodd[6*N+n] = f13;
-		disteven[8*N+n] = f16;
-		distodd[8*N+n] = f17;
+		distodd[2*N+n] = f6;
+		distodd[5*N+n] = f12;
+		disteven[7*N+n] = f13;
+		distodd[7*N+n] = f16;
+		disteven[9*N+n] = f17;
 		//...................................................
 
 	}
@@ -705,6 +714,11 @@ extern "C" void ScaLBL_D3Q19_Velocity(char *ID, double *disteven, double *distod
 				distodd[q*N+n] = -1.0;
 			}
 			disteven[9*N+n] = -1.0;
+            
+            //For ID[n]<0 - solid nodes
+			vel[n] = 0.0;
+			vel[N+n] = 0.0;
+			vel[2*N+n] = 0.0;
 		}
 	}
 }
diff --git a/gpu/D3Q19.cu b/gpu/D3Q19.cu
index b125a2ca..b137eeef 100644
--- a/gpu/D3Q19.cu
+++ b/gpu/D3Q19.cu
@@ -557,7 +557,7 @@ __global__ void dvc_D3Q19_Flux_BC_z(double *disteven, double *distodd, double fl
 	double f10,f12,f13,f16,f17;
 
 	//double A = 1.f*double(Nx*Ny);
-	double factor = 1.f/(double(Nx*Ny)*(1.0-flux));
+	double factor = 1.f/(double((Nx-2)*(Ny-2))*(1.0-flux));
 
 	double sum = 0.f;
 
@@ -567,26 +567,26 @@ __global__ void dvc_D3Q19_Flux_BC_z(double *disteven, double *distodd, double fl
 	if (n < 2*Nx*Ny){
 
 		//........................................................................
-		f1 = distodd[n];
-		f3 = distodd[N+n];
-//		f5 = distodd[2*N+n];
-		f7 = distodd[3*N+n];
-		f9 = distodd[4*N+n];
-//		f11 = distodd[5*N+n];
-		f13 = distodd[6*N+n];
-//		f15 = distodd[7*N+n];
-		f17 = distodd[8*N+n];
+		f2 = distodd[n];
+		f4 = distodd[N+n];
+		f6 = distodd[2*N+n];
+		f8 = distodd[3*N+n];
+		f10 = distodd[4*N+n];
+		f12 = distodd[5*N+n];
+		//f14 = distodd[6*N+n];
+		f16 = distodd[7*N+n];
+		//f18 = distodd[8*N+n];
 		//........................................................................
 		f0 = disteven[n];
-		f2 = disteven[N+n];
-		f4 = disteven[2*N+n];
-		f6 = disteven[3*N+n];
-		f8 = disteven[4*N+n];
-		f10 = disteven[5*N+n];
-		f12 = disteven[6*N+n];
-//		f14 = disteven[7*N+n];
-		f16 = disteven[8*N+n];
-//		f18 = disteven[9*N+n];
+		f1 = disteven[N+n];
+		f3 = disteven[2*N+n];
+		//f5 = disteven[3*N+n];
+		f7 = disteven[4*N+n];
+		f9 = disteven[5*N+n];
+		//f11 = disteven[6*N+n];
+		f13 = disteven[7*N+n];
+		//f15 = disteven[8*N+n];
+		f17 = disteven[9*N+n];
 		//...................................................
 		// compute local sum to determine the density value to set pressure
 		//sum = (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f6+f12+f13+f16+f17))/(A*(1.0-flux));
@@ -617,7 +617,7 @@ __global__ void dvc_D3Q19_Flux_BC_Z(double *disteven, double *distodd, double fl
 	N = Nx*Ny*Nz;
 	n = outlet +  blockIdx.x*blockDim.x + threadIdx.x;
 
-	double factor = 1.f/(double(Nx*Ny)*(1.0+flux));
+	double factor = 1.f/(double((Nx-2)*(Ny-2))*(1.0+flux));
 	double sum = 0.f;
 
 	// Loop over the boundary - threadblocks delineated by start...finish
@@ -625,26 +625,26 @@ __global__ void dvc_D3Q19_Flux_BC_Z(double *disteven, double *distodd, double fl
    		//........................................................................
 		// Read distributions from "opposite" memory convention
 		//........................................................................
-		f1 = distodd[n];
-		f3 = distodd[N+n];
-		f5 = distodd[2*N+n];
-		f7 = distodd[3*N+n];
-		f9 = distodd[4*N+n];
-		f11 = distodd[5*N+n];
-//		f13 = distodd[6*N+n];
-		f15 = distodd[7*N+n];
-//		f17 = distodd[8*N+n];
+		f2 = distodd[n];
+		f4 = distodd[N+n];
+		//f6 = distodd[2*N+n];
+		f8 = distodd[3*N+n];
+		f10 = distodd[4*N+n];
+		//f12 = distodd[5*N+n];
+		f14 = distodd[6*N+n];
+		//f16 = distodd[7*N+n];
+		f18 = distodd[8*N+n];
 		//........................................................................
 		f0 = disteven[n];
-		f2 = disteven[N+n];
-		f4 = disteven[2*N+n];
-//		f6 = disteven[3*N+n];
-		f8 = disteven[4*N+n];
-		f10 = disteven[5*N+n];
-//		f12 = disteven[6*N+n];
-		f14 = disteven[7*N+n];
-//		f16 = disteven[8*N+n];
-		f18 = disteven[9*N+n];
+		f1 = disteven[N+n];
+		f3 = disteven[2*N+n];
+		f5 = disteven[3*N+n];
+		f7 = disteven[4*N+n];
+		f9 = disteven[5*N+n];
+		f11 = disteven[6*N+n];
+		//f13 = disteven[7*N+n];
+		f15 = disteven[8*N+n];
+		//f17 = disteven[9*N+n];
 
 		// Local sum (based on the consistency condition)
 		//sum = (f0+f1+f2+f3+f4+f7+f8+f9+f10 + 2*(f5+f11+f14+f15+f18))/(A*(1.0+flux));
@@ -682,26 +682,26 @@ __global__  void dvc_ScaLBL_D3Q19_Pressure_BC_z(double *disteven, double *distod
 			// Read distributions from "opposite" memory convention
 			//........................................................................
 			//........................................................................
-			f1 = distodd[n];
-			f3 = distodd[N+n];
-			f5 = distodd[2*N+n];
-			f7 = distodd[3*N+n];
-			f9 = distodd[4*N+n];
-			f11 = distodd[5*N+n];
-			f13 = distodd[6*N+n];
-			f15 = distodd[7*N+n];
-			f17 = distodd[8*N+n];
-			//........................................................................
-			f0 = disteven[n];
-			f2 = disteven[N+n];
-			f4 = disteven[2*N+n];
-			f6 = disteven[3*N+n];
-			f8 = disteven[4*N+n];
-			f10 = disteven[5*N+n];
-			f12 = disteven[6*N+n];
-			f14 = disteven[7*N+n];
-			f16 = disteven[8*N+n];
-			f18 = disteven[9*N+n];
+            f2 = distodd[n];
+            f4 = distodd[N+n];
+            f6 = distodd[2*N+n];
+            f8 = distodd[3*N+n];
+            f10 = distodd[4*N+n];
+            f12 = distodd[5*N+n];
+            f14 = distodd[6*N+n];
+            f16 = distodd[7*N+n];
+            f18 = distodd[8*N+n];
+            //........................................................................
+            f0 = disteven[n];
+            f1 = disteven[N+n];
+            f3 = disteven[2*N+n];
+            f5 = disteven[3*N+n];
+            f7 = disteven[4*N+n];
+            f9 = disteven[5*N+n];
+            f11 = disteven[6*N+n];
+            f13 = disteven[7*N+n];
+            f15 = disteven[8*N+n];
+            f17 = disteven[9*N+n];
 			//...................................................
 			//........Determine the inlet flow velocity.........
 			//			uz = -1 + (f0+f3+f4+f1+f2+f7+f8+f10+f9
@@ -727,11 +727,11 @@ __global__  void dvc_ScaLBL_D3Q19_Pressure_BC_z(double *disteven, double *distod
 			f15 = f16 + 0.16666666666666678*(uy+uz)-Cyz;
 			f18 = f17 + 0.16666666666666678*(uz-uy)+Cyz;
 			//........Store in "opposite" memory location..........
-			distodd[2*N+n] = f5;
-			distodd[5*N+n] = f11;
-			disteven[7*N+n] = f14;
-			distodd[7*N+n] = f15;
-			disteven[9*N+n] = f18;
+            disteven[3*N+n] = f5;
+            disteven[6*N+n] = f11;
+            distodd[6*N+n] = f14;
+            disteven[8*N+n] = f15;
+            distodd[8*N+n] = f18;
 		}
 	}
 
@@ -753,26 +753,26 @@ __global__  void dvc_ScaLBL_D3Q19_Pressure_BC_Z(double *disteven, double *distod
 		//........................................................................
 		// Read distributions from "opposite" memory convention
 		//........................................................................
-		f1 = distodd[n];
-		f3 = distodd[N+n];
-		f5 = distodd[2*N+n];
-		f7 = distodd[3*N+n];
-		f9 = distodd[4*N+n];
-		f11 = distodd[5*N+n];
-		f13 = distodd[6*N+n];
-		f15 = distodd[7*N+n];
-		f17 = distodd[8*N+n];
+		f2 = distodd[n];
+		f4 = distodd[N+n];
+		f6 = distodd[2*N+n];
+		f8 = distodd[3*N+n];
+		f10 = distodd[4*N+n];
+		f12 = distodd[5*N+n];
+		f14 = distodd[6*N+n];
+		f16 = distodd[7*N+n];
+		f18 = distodd[8*N+n];
 		//........................................................................
 		f0 = disteven[n];
-		f2 = disteven[N+n];
-		f4 = disteven[2*N+n];
-		f6 = disteven[3*N+n];
-		f8 = disteven[4*N+n];
-		f10 = disteven[5*N+n];
-		f12 = disteven[6*N+n];
-		f14 = disteven[7*N+n];
-		f16 = disteven[8*N+n];
-		f18 = disteven[9*N+n];
+		f1 = disteven[N+n];
+		f3 = disteven[2*N+n];
+		f5 = disteven[3*N+n];
+		f7 = disteven[4*N+n];
+		f9 = disteven[5*N+n];
+		f11 = disteven[6*N+n];
+		f13 = disteven[7*N+n];
+		f15 = disteven[8*N+n];
+		f17 = disteven[9*N+n];
 		//........Determine the outlet flow velocity.........
 		//			uz = 1 - (f0+f3+f4+f1+f2+f7+f8+f10+f9+
 		//					  2*(f6+f16+f17+f12+f13))/dout;
@@ -804,11 +804,11 @@ __global__  void dvc_ScaLBL_D3Q19_Pressure_BC_Z(double *disteven, double *distod
 		f17 = f18 - 0.16666666666666678*(uz-uy)-Cyz;
 
 		//........Store in "opposite" memory location..........
-		disteven[3*N+n] = f6;
-		disteven[6*N+n] = f12;
-		distodd[6*N+n] = f13;
-		disteven[8*N+n] = f16;
-		distodd[8*N+n] = f17;
+		distodd[2*N+n] = f6;
+		distodd[5*N+n] = f12;
+		disteven[7*N+n] = f13;
+		distodd[7*N+n] = f16;
+		disteven[9*N+n] = f17;
 		//...................................................
 	}
 }
diff --git a/tests/TestFluxBC.cpp b/tests/TestFluxBC.cpp
index ac5e49d7..73a75804 100644
--- a/tests/TestFluxBC.cpp
+++ b/tests/TestFluxBC.cpp
@@ -70,6 +70,15 @@ int main (int argc, char **argv)
 		if (rank==0)	printf ("Copying phase ID to device \n");
 		char *ID;
 		ScaLBL_AllocateDeviceMemory((void **) &ID, N);						// Allocate device memory
+        // Don't compute in the halo
+        for (k=0;k<Nz;k++){
+            for (j=0;j<Ny;j++){
+                for (i=0;i<Nx;i++){
+                    int n = k*Nx*Ny+j*Nx+i;
+                    if (i==0 || i==Nx-1 || j==0 || j==Ny-1 || k==0 || k==Nz-1)	id[n] = 0;
+                }
+            }
+        }
 		// Copy to the device
 		ScaLBL_CopyToDevice(ID, id, N);
 		//...........................................................................
@@ -100,13 +109,13 @@ int main (int argc, char **argv)
 		//printf("kproc=%i \n",Dm.kproc);
 		if (pBC && Dm.kproc == 0){
 		  din = ScaLBL_D3Q19_Flux_BC_z(ID, f_even,f_odd,flux,Nx,Ny,Nz);
-			printf("Computed inlet pressure: %f \n", din);
-			ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
+		  printf("Computed inlet pressure: %.10g \n", din);
+		  ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 		}
 
 		if (pBC && Dm.kproc == nprocz-1){
-			dout = ScaLBL_D3Q19_Flux_BC_Z(f_even,f_odd,flux,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-			printf("Computed outlet pressure: %f \n", dout);
+			dout = ScaLBL_D3Q19_Flux_BC_Z(ID,f_even,f_odd,flux,Nx,Ny,Nz,Nx*Ny*(Nz-2));
+		    printf("Computed outlet pressure: %.10g \n", dout);
 			ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 		}
 
@@ -126,13 +135,17 @@ int main (int argc, char **argv)
 		for (j=0; j<Ny; j++){
 			for (i=0; i<Nx; i++){
 				n = Nx*Ny+j*Nx + i;
-				sum += vel[offset+n];
+				sum += -1.f*vel[offset+n];//extract the z-component of the velocity
+                                          //NOTE: After applying the pressure BC, the 
+                                          //distributions are stored in the "opposite" memory
+                                          //thus the pre-factor '-1.0' for getting the 
+                                          //velocity in the correct directions
 			}
 		}
 		double err;
 
 		double value;
-		value = sum/(Nx*Ny)/din;
+		value = sum/((Nx-2)*(Ny-2))/din;
 		printf("Inlet Flux: input=%f, output=%f \n",flux,value);
 		err = fabs(flux - value);
 		if (err > 1e-14){
@@ -145,10 +158,14 @@ int main (int argc, char **argv)
 		for (j=0; j<Ny; j++){
 			for (i=0; i<Nx; i++){
 				n = (Nz-2)*Nx*Ny+j*Nx + i;
-				sum += vel[offset+n];
+				sum += -1.f*vel[offset+n];//extract the z-component of the velocity
+                                          //NOTE: After applying the pressure BC, the 
+                                          //distributions are stored in the "opposite" memory
+                                          //thus the pre-factor '-1.0' for getting the 
+                                          //velocity in the correct directions
 			}
 		}
-		value = sum/(Nx*Ny)/dout;
+		value = sum/((Nx-2)*(Ny-2))/dout;
 		err = fabs(flux - value);
 		printf("Outlet Flux: input=%f, output=%f \n",flux,value);
 		err = fabs(flux - value);
diff --git a/tests/lbpm_color_macro_simulator.cpp b/tests/lbpm_color_macro_simulator.cpp
index 03ba424b..5a1ef5ae 100644
--- a/tests/lbpm_color_macro_simulator.cpp
+++ b/tests/lbpm_color_macro_simulator.cpp
@@ -648,13 +648,13 @@ int main(int argc, char **argv)
 	if (BoundaryCondition==1 && Mask.kproc == 0)	{
 		ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 		ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,0);
+		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 	}
 		
 	if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
 		ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 		ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,Nz-1);
+		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 	}
 
 	if (rank==0 && BoundaryCondition==2){
@@ -664,13 +664,13 @@ int main(int argc, char **argv)
 	if (BoundaryCondition==2 && Mask.kproc == 0)	{
 		ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 		ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,0);
+		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 	}
 
 	if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
 		ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 		ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,Nz-1);
+		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 	}
 
 	// Set dynamic pressure boundary conditions
@@ -687,13 +687,13 @@ int main(int argc, char **argv)
 		if (Mask.kproc == 0)	{
 			ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 			ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,0);
+			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 
 		}
 		if (Mask.kproc == nprocz-1){
 			ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 			ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,Nz-1);
+			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 		}
 	}
 
@@ -721,7 +721,7 @@ int main(int argc, char **argv)
 	    if (rank==0) printf("Flux = %.3e, Computed inlet pressure: %f \n",flux,din);
 	    ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 		ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,0);
+		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 
 	  }
 
@@ -729,7 +729,7 @@ int main(int argc, char **argv)
 	   // if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout);
 	    ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 		ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,Nz-1);
+		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 	  }
 	}
 
@@ -900,22 +900,24 @@ int main(int argc, char **argv)
 		if (BoundaryCondition==1 && Mask.kproc == 0)	{
 			ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 			ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 		}
 		if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
 			ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 			ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 		}
 
 		// Velocity boundary conditions
 		if (BoundaryCondition==2 && Mask.kproc == 0)	{
 			ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 			ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,0);
+			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 		}
 		if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
 			ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 			ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,Nz-1);
+			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 		}
 
 		if (BoundaryCondition==3){
@@ -927,10 +929,12 @@ int main(int argc, char **argv)
 			if (Mask.kproc == 0)	{
 				ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 				ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			    ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 			}
 			if (Mask.kproc == nprocz-1){
 				ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 				ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			    ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 			}
 		}
 
@@ -950,14 +954,14 @@ int main(int argc, char **argv)
 		  if (pBC && Dm.kproc == 0){
 		    ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 			ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,0);
+			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 		  }
 
 		  if (pBC && Dm.kproc == nprocz-1){
 		   // if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout);
 		    ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 			ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,Nz-1);
+			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 		  }
 		}
 
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index c808e668..73181948 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -638,11 +638,13 @@ int main(int argc, char **argv)
 	if (BoundaryCondition==1 && Mask.kproc == 0)	{
 		ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 		ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 	}
 		
 	if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
 		ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 		ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 	}
 
 	if (rank==0 && BoundaryCondition==2){
@@ -651,16 +653,14 @@ int main(int argc, char **argv)
 	}
 	if (BoundaryCondition==2 && Mask.kproc == 0)	{
 		ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-		//ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 		ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
-		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,0);
+		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 	}
 
 	if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
 		ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-		//ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 		ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);	
-		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,Nz-1);
+		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 	}
 	// Set dynamic pressure boundary conditions
 	double dp, slope;
@@ -676,10 +676,12 @@ int main(int argc, char **argv)
 		if (Mask.kproc == 0)	{
 			ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 			ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 		}
 		if (Mask.kproc == nprocz-1){
 			ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 			ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 		}
 	}
 
@@ -848,10 +850,12 @@ int main(int argc, char **argv)
 		if (BoundaryCondition==1 && Mask.kproc == 0)	{
 			ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 			ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 		}
 		if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
 			ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 			ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 		}
 
 		// Velocity boundary conditions
@@ -875,10 +879,12 @@ int main(int argc, char **argv)
 			if (Mask.kproc == 0)	{
 				ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
 				ScaLBL_Color_BC_z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			    ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 			}
 			if (Mask.kproc == nprocz-1){
 				ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
 				ScaLBL_Color_BC_Z(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			    ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 			}
 		}