diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index 1836c45b..fb7a471e 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -86,10 +86,10 @@ extern "C" void ScaLBL_D3Q7_ColorCollideMass(char *ID, double *A_even, double *A
 
 extern "C" void ScaLBL_ComputePhaseField(char *ID, double *Phi, double *Den, int N);
 
-extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, double *A_odd,
+extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *Vel, double *A_even, double *A_odd,
 								  double *B_even, double *B_odd, int Nx, int Ny, int Nz);
 
-extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, double *A_odd,
+extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *Vel, double *A_even, double *A_odd,
 								  double *B_even, double *B_odd, int Nx, int Ny, int Nz);
 
 extern "C" void ScaLBL_D3Q19_Velocity_BC_z(double *disteven, double *distodd, double uz,
diff --git a/cpu/Color.cpp b/cpu/Color.cpp
index 56a6ee74..f0d2723f 100644
--- a/cpu/Color.cpp
+++ b/cpu/Color.cpp
@@ -103,10 +103,12 @@ extern "C" void ScaLBL_Color_InitDistance(char *ID, double *Den, double *Phi, do
 
 
 //*************************************************************************
-extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, double *A_odd,
+extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *Velocity, double *A_even, double *A_odd,
 								  double *B_even, double *B_odd, int Nx, int Ny, int Nz)
 {
 	int i,j,k,n,N;
+	double na,nb,ux,uy,uz;
+
 	N = Nx*Ny*Nz;
 	// Fill the inlet with component a
 	for (k=0; k<1; k++){
@@ -126,6 +128,7 @@ extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, doub
 				Den[n] = 1.0;
 				Den[N+n] = 0.0;
 
+				/*
 				A_even[n] = 0.3333333333333333;
 				A_odd[n] = 0.1111111111111111;	
 				A_even[N+n] = 0.1111111111111111;
@@ -135,21 +138,80 @@ extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, doub
 				A_even[3*N+n] = 0.1111111111111111;
 				
 				B_even[n] = 0.0;
-				B_odd[n] = 0.0;	
+				B_odd[n] = 0.0;
 				B_even[N+n] = 0.0;
-				B_odd[N+n] = 0.0;	
+				B_odd[N+n] = 0.0;
 				B_even[2*N+n] = 0.0;
 				B_odd[2*N+n] = 0.0;
 				B_even[3*N+n] = 0.0;
+				*/
+
+				na = 1.0;
+				nb = 0.0;
+				//....Load the flow velocity...........
+				ux = Velocity[n];
+				uy = Velocity[N+n];
+				uz = Velocity[2*N+n];
+
+				A_even[n] = 0.3333333333333333*na;
+				B_even[n] = 0.3333333333333333*nb;
+				// Non-Stationary equilibrium distributions
+				//feq[0] = 0.1111111111111111*(1+4.5*ux);
+				//feq[1] = 0.1111111111111111*(1-4.5*ux);
+				//feq[2] = 0.1111111111111111*(1+4.5*uy);
+				//feq[3] = 0.1111111111111111*(1-4.5*uy);
+				//feq[4] = 0.1111111111111111*(1+4.5*uz);
+				//feq[5] = 0.1111111111111111*(1-4.5*uz);
+
+				//...............................................
+				// q = 0,2,4
+				// Cq = {1,0,0}, {0,1,0}, {0,0,1}
+				a1 = na*(0.1111111111111111*(1+4.5*ux));
+				b1 = nb*(0.1111111111111111*(1+4.5*ux));
+				a2 = na*(0.1111111111111111*(1-4.5*ux));
+				b2 = nb*(0.1111111111111111*(1-4.5*ux));
+
+				A_odd[n] 	= a1;
+				A_even[N+n] = a2;
+				B_odd[n] 	= b1;
+				B_even[N+n] = b2;
+				//...............................................
+				// q = 2
+				// Cq = {0,1,0}
+				a1 = na*(0.1111111111111111*(1+4.5*uy));
+				b1 = nb*(0.1111111111111111*(1+4.5*uy));
+				a2 = na*(0.1111111111111111*(1-4.5*uy));
+				b2 = nb*(0.1111111111111111*(1-4.5*uy));
+
+				A_odd[N+n] 	= a1;
+				A_even[2*N+n] = a2;
+				B_odd[N+n] 	= b1;
+				B_even[2*N+n] = b2;
+				//...............................................
+				// q = 4
+				// Cq = {0,0,1}
+				a1 = na*(0.1111111111111111*(1+4.5*uz));
+				b1 = nb*(0.1111111111111111*(1+4.5*uz));
+				a2 = na*(0.1111111111111111*(1-4.5*uz));
+				b2 = nb*(0.1111111111111111*(1-4.5*uz));
+
+				A_odd[2*N+n] = a1;
+				A_even[3*N+n] = a2;
+				B_odd[2*N+n] = b1;
+				B_even[3*N+n] = b2;
+
+
 			}					
 		}
 	}
 }
 //*************************************************************************
-extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, double *A_odd,
+extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *Velocity, double *A_even, double *A_odd,
 								  double *B_even, double *B_odd, int Nx, int Ny, int Nz)
 {
 	int i,j,k,n,N;
+	double na,nb,ux,uy,uz;
+
 	N = Nx*Ny*Nz;
 	// Fill the outlet with component b
 	for (k=Nz-3; k<Nz-1; k++){
@@ -161,7 +223,7 @@ extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, doub
 				Den[n] = 0.0;
 				Den[N+n] = 1.0;
 				
-				A_even[n] = 0.0;
+/*				A_even[n] = 0.0;
 				A_odd[n] = 0.0;	
 				A_even[N+n] = 0.0;
 				A_odd[N+n] = 0.0;	
@@ -176,6 +238,61 @@ extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, doub
 				B_even[2*N+n] = 0.1111111111111111;
 				B_odd[2*N+n] = 0.1111111111111111;
 				B_even[3*N+n] = 0.1111111111111111;
+				*/
+				na = 0.0;
+				nb = 1.0;
+				//....Load the flow velocity...........
+				ux = Velocity[n];
+				uy = Velocity[N+n];
+				uz = Velocity[2*N+n];
+
+				A_even[n] = 0.3333333333333333*na;
+				B_even[n] = 0.3333333333333333*nb;
+				// Non-Stationary equilibrium distributions
+				//feq[0] = 0.1111111111111111*(1+4.5*ux);
+				//feq[1] = 0.1111111111111111*(1-4.5*ux);
+				//feq[2] = 0.1111111111111111*(1+4.5*uy);
+				//feq[3] = 0.1111111111111111*(1-4.5*uy);
+				//feq[4] = 0.1111111111111111*(1+4.5*uz);
+				//feq[5] = 0.1111111111111111*(1-4.5*uz);
+
+				//...............................................
+				// q = 0,2,4
+				// Cq = {1,0,0}, {0,1,0}, {0,0,1}
+				a1 = na*(0.1111111111111111*(1+4.5*ux));
+				b1 = nb*(0.1111111111111111*(1+4.5*ux));
+				a2 = na*(0.1111111111111111*(1-4.5*ux));
+				b2 = nb*(0.1111111111111111*(1-4.5*ux));
+
+				A_odd[n] 	= a1;
+				A_even[N+n] = a2;
+				B_odd[n] 	= b1;
+				B_even[N+n] = b2;
+				//...............................................
+				// q = 2
+				// Cq = {0,1,0}
+				a1 = na*(0.1111111111111111*(1+4.5*uy));
+				b1 = nb*(0.1111111111111111*(1+4.5*uy));
+				a2 = na*(0.1111111111111111*(1-4.5*uy));
+				b2 = nb*(0.1111111111111111*(1-4.5*uy));
+
+				A_odd[N+n] 	= a1;
+				A_even[2*N+n] = a2;
+				B_odd[N+n] 	= b1;
+				B_even[2*N+n] = b2;
+				//...............................................
+				// q = 4
+				// Cq = {0,0,1}
+				a1 = na*(0.1111111111111111*(1+4.5*uz));
+				b1 = nb*(0.1111111111111111*(1+4.5*uz));
+				a2 = na*(0.1111111111111111*(1-4.5*uz));
+				b2 = nb*(0.1111111111111111*(1-4.5*uz));
+
+				A_odd[2*N+n] = a1;
+				A_even[3*N+n] = a2;
+				B_odd[2*N+n] = b1;
+				B_even[3*N+n] = b2;
+
 				
 			}					
 		}
diff --git a/tests/lbpm_color_macro_simulator.cpp b/tests/lbpm_color_macro_simulator.cpp
index bfefbf31..1a89d0a6 100644
--- a/tests/lbpm_color_macro_simulator.cpp
+++ b/tests/lbpm_color_macro_simulator.cpp
@@ -655,13 +655,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_Color_BC_z(Phi,Den,Velocity,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_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 	}
 
@@ -671,13 +671,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_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 		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,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 	}
 
@@ -694,13 +694,13 @@ int main(int argc, char **argv)
 		// set the initial boundary conditions
 		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_Color_BC_z(Phi,Den,Velocity,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_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 		}
 	}
@@ -728,7 +728,7 @@ int main(int argc, char **argv)
 	  if (pBC && Dm.kproc == 0){
 	    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_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 		ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 
 	  }
@@ -736,7 +736,7 @@ int main(int argc, char **argv)
 	  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_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 	  }
 	}
@@ -907,24 +907,24 @@ int main(int argc, char **argv)
 		// Pressure boundary conditions
 		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_Color_BC_z(Phi,Den,Velocity,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_Color_BC_Z(Phi,Den,Velocity,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_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 			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,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 		}
 
@@ -936,12 +936,12 @@ int main(int argc, char **argv)
 			// set the initial boundary conditions
 			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_Color_BC_z(Phi,Den,Velocity,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_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 			    ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 			}
 		}
@@ -961,14 +961,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_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 			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_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 			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 73181948..53ebb827 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -637,13 +637,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_Color_BC_z(Phi,Den,Velocity,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_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 	}
 
@@ -653,13 +653,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_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 		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,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 		ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 	}
 	// Set dynamic pressure boundary conditions
@@ -675,12 +675,12 @@ int main(int argc, char **argv)
 		// set the initial boundary conditions
 		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_Color_BC_z(Phi,Den,Velocity,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_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 		}
 	}
@@ -849,24 +849,24 @@ int main(int argc, char **argv)
 		// Pressure boundary conditions
 		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_Color_BC_z(Phi,Den,Velocity,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_Color_BC_Z(Phi,Den,Velocity,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_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 			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,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,Nz-1);
 		}
 
@@ -878,12 +878,12 @@ int main(int argc, char **argv)
 			// set the initial boundary conditions
 			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_Color_BC_z(Phi,Den,Velocity,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_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
 			    ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
 			}
 		}