diff --git a/common/ScaLBL.h b/common/ScaLBL.h
index f7722da8..f52c2760 100644
--- a/common/ScaLBL.h
+++ b/common/ScaLBL.h
@@ -91,6 +91,8 @@ extern "C" void ScaLBL_D3Q19_Velocity_BC_z(double *disteven, double *distodd, do
 extern "C" void ScaLBL_D3Q19_Velocity_BC_Z(double *disteven, double *distodd, double uz,
 								   int Nx, int Ny, int Nz, int outlet);
 
+extern "C" void SetPhiSlice_z(double *Phi, double value, int Nx, int Ny, int Nz, int Slice);
+
 class ScaLBL_Communicator{
 public:
 	//......................................................................................
diff --git a/cpu/Color.cpp b/cpu/Color.cpp
index d4eab961..65c39d8f 100644
--- a/cpu/Color.cpp
+++ b/cpu/Color.cpp
@@ -1425,6 +1425,13 @@ extern "C" void ComputePhi(char *ID, double *Phi, double *Den, int N)
 	//...................................................................
 }
 
+extern "C" void SetPhiSlice_z(double *Phi, double value, int Nx, int Ny, int Nz, int Slice){
+	int n;
+	for (n=Slice*Nx*Ny; n<(Slice+1)*Nx*Ny; n++){
+		Phi[n] = value;
+	}
+}
+
 /*
 //*************************************************************************
 extern "C" void InitDenColor( int nblocks, int nthreads, int S,
diff --git a/gpu/Color.cu b/gpu/Color.cu
index 24a881f0..a250876c 100644
--- a/gpu/Color.cu
+++ b/gpu/Color.cu
@@ -1387,17 +1387,34 @@ __global__  void dvc_ComputePhi(char *ID, double *Phi, double *Den, int N)
 	}
 	//...................................................................
 }
+
+__global__  void dvc_SetPhiSlice_z(double *Phi, double value, int Nx, int Ny, int Nz, int Slice)
+{
+	int n = Slice*Nx*Ny +  blockIdx.x*blockDim.x + threadIdx.x;
+	if (n < (Slice+1)*Nx*Ny){
+		Phi[n] = value;
+	}
+}
+
+extern "C" void SetPhiSlice_z(double *Phi, double value, int Nx, int Ny, int Nz, int Slice){
+	int GRID = Nx*Ny / 512 + 1;
+	dvc_SetPhiSlice_z<<<GRID,512>>>(Phi,value,Nx,Ny,Nz,Slice);
+}
+
 extern "C" void InitDenColor(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, int Nz){
 	dvc_InitDenColor<<<NBLOCKS,NTHREADS >>>(ID, Den, Phi, das, dbs, Nx, Ny, Nz);
 }
+
 extern "C" void InitDenColorDistance(char *ID, double *Den, double *Phi, double *Distance,
 								double das, double dbs, double beta, double xp, int Nx, int Ny, int Nz){
 
 	dvc_InitDenColorDistance<<<NBLOCKS,NTHREADS >>>(ID, Den, Phi, Distance, das, dbs, beta, xp, Nx, Ny, Nz);
 }
+
 extern "C" void ComputeColorGradient(char *ID, double *phi, double *ColorGrad, int Nx, int Ny, int Nz){
 	dvc_ComputeColorGradient<<<NBLOCKS,NTHREADS >>>(ID, phi, ColorGrad, Nx, Ny, Nz);
 }
+
 extern "C" void ColorCollide( char *ID, double *disteven, double *distodd, double *ColorGrad,
 								double *Velocity, int Nx, int Ny, int Nz,double rlx_setA, double rlx_setB,
 								double alpha, double beta, double Fx, double Fy, double Fz, bool pBC){
@@ -1405,6 +1422,7 @@ extern "C" void ColorCollide( char *ID, double *disteven, double *distodd, doubl
 									alpha, beta, Fx, Fy, Fz, pBC);
 
 }
+
 extern "C" void ColorCollideOpt( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad,
 								double *Velocity, int Nx, int Ny, int Nz,double rlx_setA, double rlx_setB,
 								double alpha, double beta, double Fx, double Fy, double Fz){
@@ -1412,6 +1430,7 @@ extern "C" void ColorCollideOpt( char *ID, double *disteven, double *distodd, do
 									alpha, beta, Fx, Fy, Fz);
 
 }
+
 extern "C" void DensityStreamD3Q7(char *ID, double *Den, double *Copy, double *Phi, double *ColorGrad, double *Velocity,
 		double beta, int Nx, int Ny, int Nz, bool pBC){
 
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index 138c5f8b..d55f7464 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -670,7 +670,6 @@ int main(int argc, char **argv)
 		//..................................................................................
 		ComputeDensityD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
 		ComputeDensityD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
-		
 		//*************************************************************************
 		// 		Compute the phase indicator field 
 		//*************************************************************************
@@ -699,11 +698,14 @@ int main(int argc, char **argv)
 		// Velocity boundary conditions
 		if (BoundaryCondition==2 && kproc == 0)	{
 			ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
-			ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			//ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			SetPhiSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 		}
 		if (BoundaryCondition==2 && kproc == nprocz-1){
 			ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
-			ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			//ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			SetPhiSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
+
 		}
 		//...................................................................................