diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index f6ef66eb..21e58066 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -3,10 +3,10 @@ color lattice boltzmann model
  */
 #include "models/ColorModel.h"
 
-ScaLBL_ColorModel::ScaLBL_ColorModel():
-  Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),rhoA(0),rhoB(0),alpha(0),beta(0),
+ScaLBL_ColorModel::ScaLBL_ColorModel(int RANK, int NP, MPI_Comm COMM):
+  rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),rhoA(0),rhoB(0),alpha(0),beta(0),
   Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0),
-  Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0)
+  Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
 {
 
 }
@@ -35,7 +35,7 @@ void ScaLBL_ColorModel::ReadParams(string filename){
     Restart = color_db->getScalar<bool>( "Restart" );
     din = color_db->getScalar<double>( "din" );
     dout = color_db->getScalar<double>( "dout" );
-    flux = color_db->getScalar<double>( "flux" );;
+    flux = color_db->getScalar<double>( "flux" );
     inletA=1.f;
     inletB=0.f;
     outletA=0.f;
@@ -65,6 +65,8 @@ void ScaLBL_ColorModel::ReadParams(string filename){
     N = Nx*Ny*Nz;
     for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1;               // initialize this way
     Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
+    
+    id = new char[N];
 
     MPI_Barrier(comm);
     Dm->CommInit(comm);
@@ -81,7 +83,6 @@ void ScaLBL_ColorModel::ReadInput(){
     sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
     sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
     // .......... READ THE INPUT FILE .......................................
-    id = new char[N];
     //...........................................................................
     if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
     //.......................................................................
@@ -201,7 +202,7 @@ void ScaLBL_ColorModel::Create(){
         
         int Npad=(Np/16 + 2)*16;
         if (rank==0)    printf ("Set up memory efficient layout \n");
-    	Map.resize(Nx,Ny,Nz);       Map.fill(0);
+    	Map.resize(Nx,Ny,Nz);       Map.fill(-2);
         auto neighborList= new int[18*Npad];
         Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
         MPI_Barrier(comm);
@@ -510,3 +511,14 @@ void ScaLBL_ColorModel::Run(){
 
     // ************************************************************************
 }
+
+void ScaLBL_ColorModel::WriteDebug(){
+	// Copy back final phase indicator field and convert to regular layout
+	DoubleArray PhaseField(Nx,Ny,Nz);
+	ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField);
+	FILE *OUTFILE;
+	sprintf(LocalRankFilename,"Phase.%05i.raw",rank);
+	OUTFILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,OUTFILE);
+	fclose(OUTFILE);
+}
\ No newline at end of file
diff --git a/models/ColorModel.h b/models/ColorModel.h
index f99b683e..17067185 100644
--- a/models/ColorModel.h
+++ b/models/ColorModel.h
@@ -18,7 +18,7 @@ Implementation of color lattice boltzmann model
 
 class ScaLBL_ColorModel{
 public:
-	ScaLBL_ColorModel();
+	ScaLBL_ColorModel(int RANK, int NP, MPI_Comm COMM);
 	~ScaLBL_ColorModel();	
 	
 	// functions in they should be run
@@ -28,6 +28,7 @@ public:
 	void Create();
 	void Initialize();
 	void Run();
+	void WriteDebug();
 	
 	bool Restart,pBC;
 	int timestep,timestepMax;
@@ -37,7 +38,7 @@ public:
 	double din,dout,inletA,inletB,outletA,outletB;
 	
 	int Nx,Ny,Nz,N,Np;
-	int nprocx,nprocy,nprocz;
+	int rank,nprocx,nprocy,nprocz,nprocs;
 	double Lx,Ly,Lz;
 
 	std::shared_ptr<Domain> Dm;   // this domain is for analysis
diff --git a/tests/TestColorBubble.cpp b/tests/TestColorBubble.cpp
index 550ba8db..96313ac7 100644
--- a/tests/TestColorBubble.cpp
+++ b/tests/TestColorBubble.cpp
@@ -12,55 +12,58 @@
 
 using namespace std;
 inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius){
-		// initialize a bubble
-  int i,j,k,n;
-  int Np=0;
-  double sum=0.f;
-  int rank = ColorModel.Mask->rank();
-		int nprocx = ColorModel.Mask->nprocx();
-		int nprocy = ColorModel.Mask->nprocy();
-		int nprocz = ColorModel.Mask->nprocz();
-		int Nx = ColorModel.Mask->Nx;
-		int Ny = ColorModel.Mask->Ny;
-		int Nz = ColorModel.Mask->Nz;
-		if (rank == 0) cout << "Setting up bubble..." << endl;
-		sum=0; Np=0;
-		for (k=0;k<Nz;k++){
-		  for (j=0;j<Ny;j++){
-		    for (i=0;i<Nx;i++){
-		      n = k*Nx*Ny + j*Nz + i;
-		      ColorModel.Averages->SDs(i,j,k) = 100.f;
-		      // Initialize phase positions field
-		      if (ColorModel.Averages->SDs(i,j,k) < 0.0){
-			ColorModel.Mask->id[n] = 0;
-		      }
-		      else {
-			sum++;
-			Np++;
-		      }
-		    }
-		  }
+	// initialize a bubble
+	int i,j,k,n;
+	int Np=0;
+	double sum=0.f;
+	int rank = ColorModel.Mask->rank();
+	int nprocx = ColorModel.Mask->nprocx();
+	int nprocy = ColorModel.Mask->nprocy();
+	int nprocz = ColorModel.Mask->nprocz();
+	int Nx = ColorModel.Mask->Nx;
+	int Ny = ColorModel.Mask->Ny;
+	int Nz = ColorModel.Mask->Nz;
+	if (rank == 0) cout << "Setting up bubble..." << endl;
+	sum=0; Np=0;
+	for (k=0;k<Nz;k++){
+		for (j=0;j<Ny;j++){
+			for (i=0;i<Nx;i++){
+				n = k*Nx*Ny + j*Nz + i;
+				ColorModel.Averages->SDs(i,j,k) = 100.f;
+				// Initialize phase positions field
+				if (ColorModel.Averages->SDs(i,j,k) < 0.0){
+					ColorModel.id[n] = 0;
+					ColorModel.id[n] = 0;
+				}
+				else {
+					sum++;
+					Np++;
+				}
+			}
 		}
-		// Initialize the bubble
-		for (k=0;k<Nz;k++){
-		  for (j=0;j<Ny;j++){
-		    for (i=0;i<Nx;i++){
-		      n = k*Nx*Ny + j*Nz + i;
-		      int iglobal= i+(Nx-2)*ColorModel.Dm->iproc();
-		      int jglobal= j+(Ny-2)*ColorModel.Dm->jproc();
-		      int kglobal= k+(Nz-2)*ColorModel.Dm->kproc();
-		      // Initialize phase position field for parallel bubble test
-		      if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
-			  +(jglobal-0.5*(Ny-2)*nprocy)*(jglobal-0.5*(Ny-2)*nprocy)
-			  +(kglobal-0.5*(Nz-2)*nprocz)*(kglobal-0.5*(Nz-2)*nprocz) < BubbleRadius*BubbleRadius){
-                        ColorModel.Mask->id[n] = 2;
-		      }
-		      else{
-                        ColorModel.Mask->id[n]=1;
-		      }
-		    }
-		  }
+	}
+	// Initialize the bubble
+	for (k=0;k<Nz;k++){
+		for (j=0;j<Ny;j++){
+			for (i=0;i<Nx;i++){
+				n = k*Nx*Ny + j*Nz + i;
+				int iglobal= i+(Nx-2)*ColorModel.Dm->iproc();
+				int jglobal= j+(Ny-2)*ColorModel.Dm->jproc();
+				int kglobal= k+(Nz-2)*ColorModel.Dm->kproc();
+				// Initialize phase position field for parallel bubble test
+				if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
+						+(jglobal-0.5*(Ny-2)*nprocy)*(jglobal-0.5*(Ny-2)*nprocy)
+						+(kglobal-0.5*(Nz-2)*nprocz)*(kglobal-0.5*(Nz-2)*nprocz) < BubbleRadius*BubbleRadius){
+					ColorModel.id[n] = 2;
+					ColorModel.id[n] = 2;
+				}
+				else{
+					ColorModel.id[n]=1;
+					ColorModel.id[n]=1;
+				}
+			}
 		}
+	}
 
 }
 //***************************************************************************************
@@ -75,7 +78,7 @@ int main(int argc, char **argv)
 	MPI_Comm comm = MPI_COMM_WORLD;
 	MPI_Comm_rank(comm,&rank);
 	MPI_Comm_size(comm,&nprocs);
-	int check;
+	int check=0;
 	{
 		if (rank == 0){
 			printf("********************************************************\n");
@@ -87,7 +90,7 @@ int main(int argc, char **argv)
 			}
 		}
 		auto filename = argv[1];
-		ScaLBL_ColorModel ColorModel;
+		ScaLBL_ColorModel ColorModel(rank,nprocs,comm);
 		ColorModel.ReadParams(filename);
 		//ColorModel.ReadInput();
 		double radius=15.5;
@@ -95,6 +98,7 @@ int main(int argc, char **argv)
 		ColorModel.Create();
 		ColorModel.Initialize();
 		ColorModel.Run();
+		ColorModel.WriteDebug();
 	}
 	// ****************************************************
 	MPI_Barrier(comm);