diff --git a/IO/Writer.cpp b/IO/Writer.cpp
index 80ba812f..012cb682 100644
--- a/IO/Writer.cpp
+++ b/IO/Writer.cpp
@@ -71,7 +71,8 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
             ERROR("Unknown mesh");
         }
         fclose(fid);
-        std::unique(mesh_entry.variables.begin(),mesh_entry.variables.end());
+        std::sort( mesh_entry.variables.begin(), mesh_entry.variables.end() );
+        mesh_entry.variables.erase( std::unique( mesh_entry.variables.begin(), mesh_entry.variables.end() ), mesh_entry.variables.end() );
         meshes_written.push_back(mesh_entry);
     }
     return meshes_written;
diff --git a/analysis/analysis.cpp b/analysis/analysis.cpp
index b54ebab8..62414872 100644
--- a/analysis/analysis.cpp
+++ b/analysis/analysis.cpp
@@ -339,9 +339,8 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_
     neighbors.push_back( rank_info.rank[1][2][1] );
     neighbors.push_back( rank_info.rank[1][1][0] );
     neighbors.push_back( rank_info.rank[1][1][2] );
-    std::unique(neighbors.begin(),neighbors.end());
-    if ( std::find(neighbors.begin(),neighbors.end(),rank) != neighbors.end() )
-        neighbors.erase(std::find(neighbors.begin(),neighbors.end(),rank));
+    std::sort( neighbors.begin(), neighbors.end() );
+    neighbors.erase( std::unique( neighbors.begin(), neighbors.end() ), neighbors.end() );
     // Create a map of all local ids to the neighbor ids
     std::map<int64_t,global_id_info_struct> map;
     std::set<int64_t> local;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 34ec9f13..e9060da6 100755
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,6 +1,7 @@
 # Copy files for the tests
 ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_color_simulator )
+ADD_LBPM_EXECUTABLE( lbpm_color_simulator_basic )
 ADD_LBPM_EXECUTABLE( lbpm_sphere_pp )
 ADD_LBPM_EXECUTABLE( lbpm_random_pp )
 ADD_LBPM_EXECUTABLE( lbpm_morphdrain_pp )
diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index 23faa789..b22032e3 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -61,18 +61,19 @@ public:
         std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
         BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ ):
         timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
-        phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_) { }
+        phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_)
+        {
+            MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
+        }
+    ~BlobIdentificationWorkItem1() { MPI_Comm_free(&newcomm); }
     virtual void run() {
         ThreadPool::WorkItem::d_state = 1;  // Change state to in progress
         // Compute the global blob id and compare to the previous version
         PROFILE_START("Identify blobs",1);
-        MPI_Comm newcomm;
-        MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
         double vF = 0.0;
         double vS = -1.0; // one voxel buffer region around solid
         IntArray& ids = new_index->second;
         new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,newcomm);
-        MPI_Comm_free(&newcomm);
         PROFILE_STOP("Identify blobs",1);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     }
@@ -85,6 +86,7 @@ private:
     const DoubleArray& dist;
     BlobIDstruct last_id, new_index, new_id;
     BlobIDList new_list;
+    MPI_Comm newcomm;
 };
 class BlobIdentificationWorkItem2: public ThreadPool::WorkItem
 {
@@ -93,13 +95,15 @@ public:
         std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
         BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ ):
         timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
-        phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_) { }
+        phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_)
+        {
+            MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
+        }
+    ~BlobIdentificationWorkItem2() { MPI_Comm_free(&newcomm); }
     virtual void run() {
         ThreadPool::WorkItem::d_state = 1;  // Change state to in progress
         // Compute the global blob id and compare to the previous version
         PROFILE_START("Identify blobs maps",1);
-        MPI_Comm newcomm;
-        MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
         const IntArray& ids = new_index->second;
         static int max_id = -1;
         new_id->first = new_index->first;
@@ -118,7 +122,6 @@ public:
             getNewIDs(map,max_id,*new_list);
             writeIDMap(map,timestep,id_map_filename);
         }
-        MPI_Comm_free(&newcomm);
         PROFILE_STOP("Identify blobs maps",1);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     }
@@ -131,6 +134,7 @@ private:
     const DoubleArray& dist;
     BlobIDstruct last_id, new_index, new_id;
     BlobIDList new_list;
+    MPI_Comm newcomm;
 };
 
 
@@ -140,7 +144,11 @@ class WriteVisWorkItem: public ThreadPool::WorkItem
 public:
     WriteVisWorkItem( int timestep_, std::vector<IO::MeshDataStruct>& visData_,
         TwoPhase& Avgerages_, fillHalo<double>& fillData_ ):
-        timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_) {}
+        timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_)
+        {
+            MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
+        }
+    ~WriteVisWorkItem() { MPI_Comm_free(&newcomm); }
     virtual void run() {
         ThreadPool::WorkItem::d_state = 1;  // Change state to in progress
         PROFILE_START("Save Vis",1);
@@ -156,10 +164,7 @@ public:
         fillData.copy(Averages.Press,PressData);
         fillData.copy(Averages.SDs,SignData);
         fillData.copy(Averages.Label_NWP,BlobData);
-        MPI_Comm newcomm;
-        MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
         IO::writeData( timestep, visData, 2, newcomm );
-        MPI_Comm_free(&newcomm);
         PROFILE_STOP("Save Vis",1);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     };
@@ -169,8 +174,10 @@ private:
     std::vector<IO::MeshDataStruct>& visData;
     TwoPhase& Averages;
     fillHalo<double>& fillData;
+    MPI_Comm newcomm;
 };
 
+
 // Helper class to run the analysis from within a thread
 // Note: Averages will be modified after the constructor is called
 class AnalysisWorkItem: public ThreadPool::WorkItem
@@ -180,6 +187,7 @@ public:
         BlobIDstruct ids, BlobIDList id_list_, double beta_ ):
         type(type_), timestep(timestep_), Averages(Averages_), 
         blob_ids(ids), id_list(id_list_), beta(beta_) { }
+    ~AnalysisWorkItem() { }
     virtual void run() {
         ThreadPool::WorkItem::d_state = 1;  // Change state to in progress
         Averages.NumberComponents_NWP = blob_ids->first;
diff --git a/tests/lbpm_color_simulator_basic.cpp b/tests/lbpm_color_simulator_basic.cpp
new file mode 100644
index 00000000..35100c44
--- /dev/null
+++ b/tests/lbpm_color_simulator_basic.cpp
@@ -0,0 +1,994 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/stat.h>
+#include <iostream>
+#include <exception>
+#include <stdexcept>
+#include <fstream>
+
+#include "common/ScaLBL.h"
+#include "common/Communication.h"
+#include "common/TwoPhase.h"
+#include "common/MPI_Helpers.h"
+
+#include "ProfilerApp.h"
+//#include "threadpool/thread_pool.h"
+
+//#include "lbpm_color_simulator.h"
+
+/*
+ * Simulator for two-phase flow in porous media
+ * James E. McClure 2013-2016
+ */
+
+using namespace std;
+
+//*************************************************************************
+// Implementation of Two-Phase Immiscible LBM using CUDA
+//*************************************************************************
+inline void PackID(int *list, int count, char *sendbuf, char *ID){
+	// Fill in the phase ID values from neighboring processors
+	// This packs up the values that need to be sent from one processor to another
+	int idx,n;
+
+	for (idx=0; idx<count; idx++){
+		n = list[idx];
+		sendbuf[idx] = ID[n];
+	}
+}
+//***************************************************************************************
+
+inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
+	// Fill in the phase ID values from neighboring processors
+	// This unpacks the values once they have been recieved from neighbors
+	int idx,n;
+
+	for (idx=0; idx<count; idx++){
+		n = list[idx];
+		ID[n] = recvbuf[idx];
+	}
+}
+
+//***************************************************************************************
+
+inline void ZeroHalo(double *Data, int Nx, int Ny, int Nz)
+{
+	int i,j,k,n;
+	for (k=0;k<Nz;k++){
+		for (j=0;j<Ny;j++){
+			i=0;
+			n = k*Nx*Ny+j*Nx+i;
+			Data[2*n] = 0.0;
+			Data[2*n+1] = 0.0;
+			i=Nx-1;
+			n = k*Nx*Ny+j*Nx+i;
+			Data[2*n] = 0.0;
+			Data[2*n+1] = 0.0;
+		}
+	}
+
+	for (k=0;k<Nz;k++){
+		for (i=0;i<Nx;i++){
+			j=0;
+			n = k*Nx*Ny+j*Nx+i;
+			Data[2*n] = 0.0;
+			Data[2*n+1] = 0.0;
+			j=Ny-1;
+			n = k*Nx*Ny+j*Nx+i;
+			Data[2*n] = 0.0;
+			Data[2*n+1] = 0.0;
+		}
+	}
+
+	for (j=0;j<Ny;j++){
+		for (i=0;i<Nx;i++){
+			k=0;
+			n = k*Nx*Ny+j*Nx+i;
+			Data[2*n] = 0.0;
+			Data[2*n+1] = 0.0;
+			k=Nz-1;
+			n = k*Nx*Ny+j*Nx+i;
+			Data[2*n] = 0.0;
+			Data[2*n+1] = 0.0;
+		}
+	}
+}
+//***************************************************************************************
+
+
+int main(int argc, char **argv)
+{
+  // Initialize MPI
+  int provided_thread_support = -1;
+  MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
+  MPI_Comm comm;
+  MPI_Comm_dup(MPI_COMM_WORLD,&comm);
+  int rank = comm_rank(comm);
+  int nprocs = comm_size(comm);
+  if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE )
+    std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
+  { // Limit scope so variables that contain communicators will free before MPI_Finialize
+
+	// parallel domain size (# of sub-domains)
+	int nprocx,nprocy,nprocz;
+	int iproc,jproc,kproc;
+
+	MPI_Request req1[18],req2[18];
+	MPI_Status stat1[18],stat2[18];
+
+	if (rank == 0){
+		printf("********************************************************\n");
+		printf("Running Color LBM	\n");
+		printf("********************************************************\n");
+	}
+
+    PROFILE_ENABLE(1);
+    //PROFILE_ENABLE_TRACE();
+    //PROFILE_ENABLE_MEMORY();
+    PROFILE_SYNCHRONIZE();
+    PROFILE_START("Main");
+    Utilities::setErrorHandlers();
+
+	// Variables that specify the computational domain  
+	string FILENAME;
+	unsigned int nBlocks, nthreads;
+	int Nx,Ny,Nz;		// local sub-domain size
+	int nspheres;		// number of spheres in the packing
+	double Lx,Ly,Lz;	// Domain length
+	double D = 1.0;		// reference length for non-dimensionalization
+	// Color Model parameters
+	int timestepMax;
+	double tau,Fx,Fy,Fz,tol,err;
+	double alpha, beta;
+	double das, dbs, phi_s;
+	double din,dout;
+	double wp_saturation;
+	int BoundaryCondition;
+	int InitialCondition;
+//	bool pBC,Restart;
+	int i,j,k;
+
+	// pmmc threshold values
+	//double fluid_isovalue,solid_isovalue;
+	//fluid_isovalue = 0.0;
+	//solid_isovalue = 0.0;
+	
+	int RESTART_INTERVAL=20000;
+	//int ANALYSIS_)INTERVAL=1000;	
+	int BLOB_ANALYSIS_INTERVAL=1000;
+	int timestep = -1;
+
+	if (rank==0){
+		//.............................................................
+		//		READ SIMULATION PARMAETERS FROM INPUT FILE
+		//.............................................................
+		ifstream input("Color.in");
+		if (input.is_open()){
+			// Line 1: model parameters (tau, alpha, beta, das, dbs)
+			input >> tau;			// Viscosity parameter
+			input >> alpha;			// Surface Tension parameter
+			input >> beta;			// Width of the interface
+			input >> phi_s;			// value of phi at the solid surface
+			// Line 2: wetting phase saturation to initialize
+			input >> wp_saturation;
+			// Line 3: External force components (Fx,Fy, Fz)
+			input >> Fx;
+			input >> Fy;
+			input >> Fz;
+			// Line 4: Pressure Boundary conditions
+			input >> InitialCondition;
+			input >> BoundaryCondition;
+			input >> din;
+			input >> dout;
+			// Line 5: time-stepping criteria
+			input >> timestepMax;		// max no. of timesteps
+			input >> RESTART_INTERVAL;	// restart interval
+			input >> tol;		      	// error tolerance
+			// Line 6: Analysis options
+			input >> BLOB_ANALYSIS_INTERVAL; // interval to analyze blob states
+			//.............................................................
+		}
+		else{
+			// Set default values
+		    // Print warning
+			printf("WARNING: No input file provided (Color.in is missing)! Default parameters will be used. \n");
+			tau = 1.0;
+			alpha=0.005;
+			beta= 0.9;
+			Fx = Fy = Fz = 0.0;
+			InitialCondition=0;
+			BoundaryCondition=0;
+			din=dout=1.0;
+			timestepMax=0;
+		}
+
+		//.......................................................................
+		// Reading the domain information file
+		//.......................................................................
+		ifstream domain("Domain.in");
+		if (input.is_open()){
+
+			domain >> nprocx;
+			domain >> nprocy;
+			domain >> nprocz;
+			domain >> Nx;
+			domain >> Ny;
+			domain >> Nz;
+			domain >> nspheres;
+			domain >> Lx;
+			domain >> Ly;
+			domain >> Lz;
+			//.......................................................................
+		}
+		else{
+			// Set default values
+		    // Print warning
+			printf("WARNING: No input file provided (Domain.in is missing)! Default parameters will be used. \n");
+			nprocx=nprocy=nprocz=1;
+			Nx=Ny=Nz=10;
+			nspheres=0;
+			Lx=Ly=Lz=1.0;
+		}
+	}
+	// **************************************************************
+	// Broadcast simulation parameters from rank 0 to all other procs
+	MPI_Barrier(comm);
+	//.................................................
+	MPI_Bcast(&tau,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&alpha,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&beta,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&das,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&dbs,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&phi_s,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&wp_saturation,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&BoundaryCondition,1,MPI_INT,0,comm);
+	MPI_Bcast(&InitialCondition,1,MPI_INT,0,comm);
+	MPI_Bcast(&din,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&dout,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&Fx,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&Fy,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&Fz,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&timestepMax,1,MPI_INT,0,comm);
+	MPI_Bcast(&RESTART_INTERVAL,1,MPI_INT,0,comm);
+	MPI_Bcast(&tol,1,MPI_DOUBLE,0,comm);
+	// Computational domain
+	MPI_Bcast(&Nx,1,MPI_INT,0,comm);
+	MPI_Bcast(&Ny,1,MPI_INT,0,comm);
+	MPI_Bcast(&Nz,1,MPI_INT,0,comm);
+	MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
+	MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
+	MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
+	MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
+	MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm);
+	//.................................................
+
+    // Get the rank info
+    const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
+
+	MPI_Barrier(comm);
+	
+	// **************************************************************
+	// **************************************************************
+	double Ps = -(das-dbs)/(das+dbs);
+	double rlxA = 1.f/tau;
+	double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA);
+	//double xIntPos = log((1.0+phi_s)/(1.0-phi_s))/(2.0*beta); 	
+	
+	// Set the density values inside the solid based on the input value phi_s
+ 	das = (phi_s+1.0)*0.5;
+	dbs = 1.0 - das;
+	
+	if (nprocs != nprocx*nprocy*nprocz){
+		printf("nprocx =  %i \n",nprocx);
+		printf("nprocy =  %i \n",nprocy);
+		printf("nprocz =  %i \n",nprocz);
+		INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!");
+	}
+
+	if (rank==0){
+		printf("********************************************************\n");
+		printf("tau = %f \n", tau);
+		printf("alpha = %f \n", alpha);		
+		printf("beta = %f \n", beta);
+		printf("das = %f \n", das);
+		printf("dbs = %f \n", dbs);
+		printf("gamma_{wn} = %f \n", 5.796*alpha);
+		printf("Force(x) = %f \n", Fx);
+		printf("Force(y) = %f \n", Fy);
+		printf("Force(z) = %f \n", Fz);
+		printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
+		printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
+		if (BoundaryCondition==0) printf("Periodic boundary conditions will applied \n");
+		if (BoundaryCondition==1) printf("Pressure boundary conditions will be applied \n");
+		if (BoundaryCondition==2) printf("Velocity boundary conditions will be applied \n");
+		if (BoundaryCondition==3) printf("Dynamic pressure boundary conditions will be applied \n");
+		if (InitialCondition==0) printf("Initial conditions assigned from phase ID file \n");
+		if (InitialCondition==1) printf("Initial conditions assigned from restart file \n");
+		printf("********************************************************\n");
+	}
+
+	// Initialized domain and averaging framework for Two-Phase Flow
+	bool pBC,velBC;
+	if (BoundaryCondition==1 || BoundaryCondition==3)
+								pBC=true;
+	else						pBC=false;
+	if (BoundaryCondition==2)	velBC=true;
+	else						velBC=false;
+
+	bool Restart;
+	if (InitialCondition==1)    Restart=true;
+	else 						Restart=false;
+    NULL_USE(pBC);  NULL_USE(velBC);
+
+	// Full domain used for averaging (do not use mask for analysis)
+	Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
+	for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
+	//	std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
+        TwoPhase Averages(Dm);
+	Dm.CommInit(comm);
+
+	// Mask that excludes the solid phase
+	Domain Mask(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
+
+	MPI_Barrier(comm);
+
+	Nx+=2; Ny+=2; Nz += 2;
+	//Nx = Ny = Nz;	// Cubic domain
+
+	int N = Nx*Ny*Nz;
+	int dist_mem_size = N*sizeof(double);
+
+	//.......................................................................
+	if (rank == 0)	printf("Read input media... \n");
+	//.......................................................................
+	
+	//.......................................................................
+	// Filenames used
+	char LocalRankString[8];
+	char LocalRankFilename[40];
+	char LocalRestartFile[40];
+	char tmpstr[10];
+	sprintf(LocalRankString,"%05d",rank);
+	sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
+	sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
+	
+//	printf("Local File Name =  %s \n",LocalRankFilename);
+	// .......... READ THE INPUT FILE .......................................
+//	char value;
+	char *id;
+	id = new char[N];
+	int sum = 0;
+	double sum_local;
+	double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
+	if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
+	double porosity, pore_vol;
+	//...........................................................................
+	if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
+
+	//.......................................................................
+	sprintf(LocalRankString,"%05d",rank);
+//	sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
+//	WriteLocalSolidID(LocalRankFilename, id, N);
+	sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
+	ReadBinaryFile(LocalRankFilename, Averages->SDs.data(), N);
+	MPI_Barrier(comm);
+	if (rank == 0) cout << "Domain set." << endl;
+	
+	//.......................................................................
+	// Assign the phase ID field based on the signed distance
+	//.......................................................................
+	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;
+				id[n] = 0;
+			}
+		}
+	}
+	sum=0;
+	pore_vol = 0.0;
+	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 (Averages->SDs(n) > 0.0){
+					id[n] = 2;	
+				}
+				// compute the porosity (actual interface location used)
+				if (Averages->SDs(n) > 0.0){
+					sum++;	
+				}
+			}
+		}
+	}
+
+	if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
+	sprintf(LocalRankFilename,"ID.%05i",rank);
+	size_t readID;
+	FILE *IDFILE = fopen(LocalRankFilename,"rb");
+	if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx");
+	readID=fread(id,1,N,IDFILE);
+	if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID (rank=%i) \n",rank);
+
+	fclose(IDFILE);
+
+	// Set up kstart, kfinish so that the reservoirs are excluded from averaging
+	int kstart,kfinish;
+	kstart = 1;
+	kfinish = Nz-1;
+	if (BoundaryCondition >  0 && Dm.kproc==0)		kstart = 4;
+	if (BoundaryCondition >  0 && Dm.kproc==nprocz-1)	kfinish = Nz-4;
+
+	// Compute the pore volume
+	sum_local = 0.0;
+	for ( k=kstart;k<kfinish;k++){
+		for ( j=1;j<Ny-1;j++){
+			for ( i=1;i<Nx-1;i++){
+				int n = k*Nx*Ny+j*Nx+i;
+				if (id[n] > 0){
+					sum_local += 1.0;
+				}
+			}
+		}
+	}
+	MPI_Allreduce(&sum_local,&pore_vol,1,MPI_DOUBLE,MPI_SUM,comm);
+	porosity = pore_vol*iVol_global;
+	if (rank==0) printf("Media porosity = %f \n",porosity);
+	//.........................................................
+	// If external boundary conditions are applied remove solid
+	if (BoundaryCondition >  0  && Dm.kproc == 0){
+		for (k=0; k<3; k++){
+			for (j=0;j<Ny;j++){
+				for (i=0;i<Nx;i++){
+					int n = k*Nx*Ny+j*Nx+i;
+					//id[n] = 1;
+					Averages->SDs(n) = max(Averages->SDs(n),1.0*(2.5-k));
+				}					
+			}
+		}
+	}
+	if (BoundaryCondition >  0  && Dm.kproc == nprocz-1){
+		for (k=Nz-3; k<Nz; k++){
+			for (j=0;j<Ny;j++){
+				for (i=0;i<Nx;i++){
+					int n = k*Nx*Ny+j*Nx+i;
+					//id[n] = 2;
+					Averages->SDs(n) = max(Averages->SDs(n),1.0*(k-Nz+2.5));
+				}					
+			}
+		}
+	}
+	//.........................................................
+	// don't perform computations at the eight corners
+	id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0;
+	id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
+	//.........................................................
+
+	// Initialize communication structures in averaging domain
+	for (i=0; i<Mask.Nx*Mask.Ny*Mask.Nz; i++) Mask.id[i] = id[i];
+	Mask.CommInit(comm);
+
+	//...........................................................................
+	if (rank==0)	printf ("Create ScaLBL_Communicator \n");
+	// Create a communicator for the device
+	ScaLBL_Communicator ScaLBL_Comm(Mask);
+
+	// set reservoirs
+	if (BoundaryCondition > 0){
+		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 (Dm.kproc==0 && k==0)			id[n]=1;
+					if (Dm.kproc==0 && k==1)			id[n]=1;
+					if (Dm.kproc==nprocz-1 && k==Nz-2)	id[n]=2;
+					if (Dm.kproc==nprocz-1 && k==Nz-1)	id[n]=2;
+					Mask.id[n] = id[n];
+				}
+			}
+		}
+	}
+
+	//...........device phase ID.................................................
+	if (rank==0)	printf ("Copy phase ID to device \n");
+	char *ID;
+	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
+	CopyToDevice(ID, id, N);
+	DeviceBarrier();
+	//...........................................................................
+
+	//...........................................................................
+	//				MAIN  VARIABLES ALLOCATED HERE
+	//...........................................................................
+	// LBM variables
+	if (rank==0)	printf ("Allocating distributions \n");
+	//......................device distributions.................................
+	double *f_even,*f_odd;
+	double *A_even,*A_odd,*B_even,*B_odd;
+	//...........................................................................
+	AllocateDeviceMemory((void **) &f_even, 10*dist_mem_size);	// Allocate device memory
+	AllocateDeviceMemory((void **) &f_odd, 9*dist_mem_size);	// Allocate device memory
+	AllocateDeviceMemory((void **) &A_even, 4*dist_mem_size);	// Allocate device memory
+	AllocateDeviceMemory((void **) &A_odd, 3*dist_mem_size);	// Allocate device memory
+	AllocateDeviceMemory((void **) &B_even, 4*dist_mem_size);	// Allocate device memory
+	AllocateDeviceMemory((void **) &B_odd, 3*dist_mem_size);	// Allocate device memory
+	//...........................................................................
+	double *Phi,*Den;
+	double *ColorGrad, *Velocity, *Pressure, *dvcSignDist;
+	//...........................................................................
+	AllocateDeviceMemory((void **) &Phi, dist_mem_size);
+	AllocateDeviceMemory((void **) &Pressure, dist_mem_size);
+	AllocateDeviceMemory((void **) &dvcSignDist, dist_mem_size);
+	AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
+	AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
+	AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size);
+	//...........................................................................
+
+	// Copy signed distance for device initialization
+	CopyToDevice(dvcSignDist, Averages->SDs.data(), dist_mem_size);
+	//...........................................................................
+
+	int logcount = 0; // number of surface write-outs
+	
+	//...........................................................................
+	//				MAIN  VARIABLES INITIALIZED HERE
+	//...........................................................................
+	//...........................................................................
+	//...........................................................................
+	if (rank==0)	printf("Setting the distributions, size = %i\n", N);
+	//...........................................................................
+	DeviceBarrier();
+	InitD3Q19(ID, f_even, f_odd, Nx, Ny, Nz);
+	InitDenColor(ID, Den, Phi, das, dbs, Nx, Ny, Nz);
+	DeviceBarrier();
+	//......................................................................
+
+	if (Restart == true){
+
+	  
+	  if (rank==0){
+	    printf("Reading restart file! \n");
+	    ifstream restart("Restart.txt");
+	    if (restart.is_open()){
+	      restart  >> timestep;
+	      printf("Restarting from timestep =%i \n",timestep);
+	    }
+	    else{
+	      printf("WARNING:No Restart.txt file, setting timestep=0 \n");
+	      timestep=0;
+	    }
+	  }
+	  MPI_Bcast(&timestep,1,MPI_INT,0,comm);
+
+		// Read in the restart file to CPU buffers
+	    double *cDen = new double[2*N];
+	    double *cDistEven = new double[10*N];
+	    double *cDistOdd = new double[9*N];
+		ReadCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
+		// Copy the restart data to the GPU
+		CopyToDevice(f_even,cDistEven,10*N*sizeof(double));
+		CopyToDevice(f_odd,cDistOdd,9*N*sizeof(double));
+		CopyToDevice(Den,cDen,2*N*sizeof(double));
+		DeviceBarrier();
+	    delete [] cDen;
+	    delete [] cDistEven;
+	    delete [] cDistOdd;
+		MPI_Barrier(comm);
+	}
+
+	//......................................................................
+	InitD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
+	InitD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
+	DeviceBarrier();
+	MPI_Barrier(comm);
+	//.......................................................................
+	// Once phase has been initialized, map solid to account for 'smeared' interface
+	for (i=0; i<N; i++)	Averages->SDs(i) -= (1.0);
+	// Make sure the id match for the two domains
+	for (i=0; i<N; i++)	Dm.id[i] = Mask.id[i];
+	//.......................................................................
+	// Finalize setup for averaging domain
+	//Averages->SetupCubes(Mask);
+	Averages->UpdateSolid();
+	//.......................................................................
+	
+	//*************************************************************************
+	// 		Compute the phase indicator field and reset Copy, Den
+	//*************************************************************************
+	ComputePhi(ID, Phi, Den, N);
+	//*************************************************************************
+	DeviceBarrier();
+	ScaLBL_Comm.SendHalo(Phi);
+	ScaLBL_Comm.RecvHalo(Phi);
+	DeviceBarrier();
+	MPI_Barrier(comm);
+	//*************************************************************************
+
+	if (rank==0 && BoundaryCondition==1){
+		printf("Setting inlet pressure = %f \n", din);
+		printf("Setting outlet pressure = %f \n", dout);
+	}
+	if (BoundaryCondition==1 && Mask.kproc == 0)	{
+		PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
+		ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+	}
+		
+	if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
+		PressureBC_outlet(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);
+	}
+
+	if (rank==0 && BoundaryCondition==2){
+		printf("Setting inlet velocity = %f \n", din);
+		printf("Setting outlet velocity = %f \n", dout);
+	}
+	if (BoundaryCondition==2 && Mask.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);
+		SetPhiSlice_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));
+		//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);
+	}
+
+	// Set dynamic pressure boundary conditions
+	double dp, slope;
+	dp = slope = 0.0;
+	if (BoundaryCondition==3){
+		slope = (dout-din)/timestepMax;
+		dp = din + timestep*slope;
+		if (rank==0) printf("Change in pressure / time =%f \n",slope);
+		// set the initial value
+		din  = 1.0+0.5*dp;
+		dout = 1.0-0.5*dp;
+		// set the initial boundary conditions
+		if (Mask.kproc == 0)	{
+			PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
+			ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+		}
+		if (Mask.kproc == nprocz-1){
+			PressureBC_outlet(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);
+		}
+	}
+
+	ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
+	ComputeVelocityD3Q19(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
+
+	//...........................................................................
+	// Copy the phase indicator field for the earlier timestep
+	DeviceBarrier();
+	CopyToHost(Averages->Phase_tplus.data(),Phi,N*sizeof(double));
+	//...........................................................................
+	//...........................................................................
+	// Copy the data for for the analysis timestep
+	//...........................................................................
+	// Copy the phase from the GPU -> CPU
+	//...........................................................................
+	DeviceBarrier();
+	ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
+	CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
+	CopyToHost(Averages->Press.data(),Pressure,N*sizeof(double));
+	CopyToHost(Averages->Vel_x.data(),&Velocity[0],N*sizeof(double));
+	CopyToHost(Averages->Vel_y.data(),&Velocity[N],N*sizeof(double));
+	CopyToHost(Averages->Vel_z.data(),&Velocity[2*N],N*sizeof(double));
+	//...........................................................................
+	
+	if (rank==0) printf("********************************************************\n");
+	if (rank==0)	printf("No. of timesteps: %i \n", timestepMax);
+
+	//.......create and start timer............
+	double starttime,stoptime,cputime;
+	DeviceBarrier();
+	MPI_Barrier(comm);
+	starttime = MPI_Wtime();
+	//.........................................
+	
+	err = 1.0; 	
+	double sat_w_previous = 1.01; // slightly impossible value!
+	if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol);
+
+    // Create the MeshDataStruct
+    fillHalo<double> fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
+    std::vector<IO::MeshDataStruct> meshData(1);
+    meshData[0].meshName = "domain";
+    meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );
+    std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
+    std::shared_ptr<IO::Variable> PressVar( new IO::Variable() );
+    std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
+    std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
+    PhaseVar->name = "phase";
+    PhaseVar->type = IO::VolumeVariable;
+    PhaseVar->dim = 1;
+    PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(PhaseVar);
+    PressVar->name = "Pressure";
+    PressVar->type = IO::VolumeVariable;
+    PressVar->dim = 1;
+    PressVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(PressVar);
+    SignDistVar->name = "SignDist";
+    SignDistVar->type = IO::VolumeVariable;
+    SignDistVar->dim = 1;
+    SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(SignDistVar);
+    BlobIDVar->name = "BlobID";
+    BlobIDVar->type = IO::VolumeVariable;
+    BlobIDVar->dim = 1;
+    BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(BlobIDVar);
+
+	//************ MAIN ITERATION LOOP ***************************************/
+    PROFILE_START("Loop");
+
+	while (timestep < timestepMax && err > tol ) {
+        //if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
+        PROFILE_START("Update");
+
+		//*************************************************************************
+		// Fused Color Gradient and Collision 
+		//*************************************************************************
+		ColorCollideOpt( ID,f_even,f_odd,Phi,ColorGrad,
+							 Velocity,Nx,Ny,Nz,rlxA,rlxB,alpha,beta,Fx,Fy,Fz);
+		//*************************************************************************
+
+		DeviceBarrier();
+		//*************************************************************************
+		// Pack and send the D3Q19 distributions
+		ScaLBL_Comm.SendD3Q19(f_even, f_odd);
+		//*************************************************************************
+
+		//*************************************************************************
+		// 		Carry out the density streaming step for mass transport
+		//*************************************************************************
+		MassColorCollideD3Q7(ID, A_even, A_odd, B_even, B_odd, Den, Phi,
+								ColorGrad, Velocity, beta, N, pBC);
+		//*************************************************************************
+
+		DeviceBarrier();
+		MPI_Barrier(comm);
+		//*************************************************************************
+		// 		Swap the distributions for momentum transport
+		//*************************************************************************
+		SwapD3Q19(ID, f_even, f_odd, Nx, Ny, Nz);
+		//*************************************************************************
+
+		DeviceBarrier();
+		MPI_Barrier(comm);
+		//*************************************************************************
+		// Wait for communications to complete and unpack the distributions
+		ScaLBL_Comm.RecvD3Q19(f_even, f_odd);
+		//*************************************************************************
+
+		DeviceBarrier();
+		//*************************************************************************
+		// Pack and send the D3Q7 distributions
+		ScaLBL_Comm.BiSendD3Q7(A_even, A_odd, B_even, B_odd);
+		//*************************************************************************
+
+		DeviceBarrier();
+		SwapD3Q7(ID, A_even, A_odd, Nx, Ny, Nz);
+		SwapD3Q7(ID, B_even, B_odd, Nx, Ny, Nz);
+
+		DeviceBarrier();
+		MPI_Barrier(comm);
+
+		//*************************************************************************
+		// Wait for communication and unpack the D3Q7 distributions
+		ScaLBL_Comm.BiRecvD3Q7(A_even, A_odd, B_even, B_odd);
+		//*************************************************************************
+
+		DeviceBarrier();
+		//..................................................................................
+		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 
+		//*************************************************************************
+		DeviceBarrier();
+		MPI_Barrier(comm);
+
+		ComputePhi(ID, Phi, Den, N);
+		//*************************************************************************
+		ScaLBL_Comm.SendHalo(Phi);
+		DeviceBarrier();
+		ScaLBL_Comm.RecvHalo(Phi);
+		//*************************************************************************
+
+		DeviceBarrier();
+		
+		// Pressure boundary conditions
+		if (BoundaryCondition==1 && Mask.kproc == 0)	{
+			PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
+			ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+		}
+		if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
+			PressureBC_outlet(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);
+		}
+
+		// Velocity boundary conditions
+		if (BoundaryCondition==2 && Mask.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);
+			SetPhiSlice_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));
+			//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);
+		}
+
+		if (BoundaryCondition==3){
+			// Increase the pressure difference
+			dp += slope;
+			din  = 1.0+0.5*dp;
+			dout = 1.0-0.5*dp;
+			// set the initial boundary conditions
+			if (Mask.kproc == 0)	{
+				PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
+				ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
+			}
+			if (Mask.kproc == nprocz-1){
+				PressureBC_outlet(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);
+			}
+		}
+
+		//...................................................................................
+
+		MPI_Barrier(comm);
+        PROFILE_STOP("Update");
+
+		// Timestep completed!
+		timestep++;
+  
+        // Run the analysis
+	//if (timestep > 5)
+		//...................................................................
+		if (timestep%1000 == 995){
+			//...........................................................................
+			// Copy the phase indicator field for the earlier timestep
+			DeviceBarrier();
+			CopyToHost(Averages.Phase_tplus,Phi,N*sizeof(double));
+	//		Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
+			//...........................................................................
+		}
+		if (timestep%1000 == 0){
+			//...........................................................................
+			// Copy the data for for the analysis timestep
+			//...........................................................................
+			// Copy the phase from the GPU -> CPU
+			//...........................................................................
+			DeviceBarrier();
+			ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
+			CopyToHost(Averages.Phase,Phi,N*sizeof(double));
+			CopyToHost(Averages.Press,Pressure,N*sizeof(double));
+			CopyToHost(Averages.Vel_x,&Velocity[0],N*sizeof(double));
+			CopyToHost(Averages.Vel_y,&Velocity[N],N*sizeof(double));
+			CopyToHost(Averages.Vel_z,&Velocity[2*N],N*sizeof(double));
+			MPI_Barrier(MPI_COMM_WORLD);
+		}
+		if (timestep%1000 == 5){
+			//...........................................................................
+			// Copy the phase indicator field for the later timestep
+			DeviceBarrier();
+			CopyToHost(Averages.Phase_tminus,Phi,N*sizeof(double));
+//			Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus);
+			//....................................................................
+            Averages.Initialize();
+            Averages.ComputeDelPhi();
+            Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn);
+            Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus);
+            Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus);
+            Averages.UpdateMeshValues();
+            Averages.ComputeLocal();
+            Averages.Reduce();
+            Averages.PrintAll(timestep);
+*/			//....................................................................
+		}
+
+		if (timestep%RESTART_INTERVAL == 0){
+			if (pBC){
+				//err = fabs(sat_w - sat_w_previous);
+				//sat_w_previous = sat_w;
+				if (rank==0) printf("Timestep %i: change in saturation since last checkpoint is %f \n", timestep, err);
+			}
+			else{
+				// Not clear yet
+			}
+			// Copy the data to the CPU
+			CopyToHost(cDistEven,f_even,10*N*sizeof(double));
+			CopyToHost(cDistOdd,f_odd,9*N*sizeof(double));
+			CopyToHost(cDen,Den,2*N*sizeof(double));
+			// Read in the restart file to CPU buffers
+			WriteCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
+		}
+
+
+        // Save the timers
+        if ( timestep%50==0 )
+            PROFILE_SAVE("lbpm_color_simulator",1);
+	}
+    PROFILE_STOP("Loop");
+	//************************************************************************
+	DeviceBarrier();
+	MPI_Barrier(comm);
+	stoptime = MPI_Wtime();
+	if (rank==0) printf("-------------------------------------------------------------------\n");
+	// Compute the walltime per timestep
+	cputime = (stoptime - starttime)/timestep;
+	// Performance obtained from each node
+	double MLUPS = double(Nx*Ny*Nz)/cputime/1000000;
+	
+	if (rank==0) printf("********************************************************\n");
+	if (rank==0) printf("CPU time = %f \n", cputime);
+	if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
+	MLUPS *= nprocs;
+	if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
+	if (rank==0) printf("********************************************************\n");
+	
+	// ************************************************************************
+/*	// Perform component averaging and write tcat averages
+	Averages->Initialize();
+	Averages->ComponentAverages();
+	Averages->SortBlobs();
+	Averages->PrintComponents(timestep);
+	// ************************************************************************
+
+	int NumberComponents_NWP = ComputeGlobalPhaseComponent(Mask.Nx-2,Mask.Ny-2,Mask.Nz-2,Mask.rank_info,Averages->PhaseID,1,Averages->Label_NWP);
+	printf("Number of non-wetting phase components: %i \n ",NumberComponents_NWP);
+	DeviceBarrier();
+	CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
+*/
+    
+/*	Averages->WriteSurfaces(0);
+
+	sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
+	FILE *PHASE;
+	PHASE = fopen(LocalRankFilename,"wb");
+	fwrite(Averages->SDn.data(),8,N,PHASE);
+	fclose(PHASE);
+	*/
+
+	/*	sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
+	FILE *PRESS;
+	PRESS = fopen(LocalRankFilename,"wb");
+	fwrite(Averages->Press.data(),8,N,PRESS);
+	fclose(PRESS);
+
+	CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
+	double * Grad;
+	Grad = new double [3*N];
+	CopyToHost(Grad,ColorGrad,3*N*sizeof(double));
+	sprintf(LocalRankFilename,"%s%s","ColorGrad.",LocalRankString);
+	FILE *GRAD;
+	GRAD = fopen(LocalRankFilename,"wb");
+	fwrite(Grad,8,3*N,GRAD);
+	fclose(GRAD);
+	*/
+    PROFILE_STOP("Main");
+    PROFILE_SAVE("lbpm_color_simulator",1);
+	// ****************************************************
+	MPI_Barrier(comm);
+  } // Limit scope so variables that contain communicators will free before MPI_Finialize
+  MPI_Comm_free(&comm);
+  MPI_Finalize();
+}
+
+
diff --git a/tests/lbpm_permeability_simulator.cpp b/tests/lbpm_permeability_simulator.cpp
index 5f693af0..b54eb4bd 100644
--- a/tests/lbpm_permeability_simulator.cpp
+++ b/tests/lbpm_permeability_simulator.cpp
@@ -229,7 +229,7 @@ int main(int argc, char **argv)
 		printf("Force(y) = %f \n", Fy);
 		printf("Force(z) = %f \n", Fz);
 		printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
-		printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
+		printf("Process grid = %i x %i x %i\n",nprocx,nprocy,nprocz);
 		printf("********************************************************\n");
 	}
 
@@ -247,7 +247,6 @@ int main(int argc, char **argv)
 	MPI_Barrier(comm);
 
 	Nx += 2;	Ny += 2;	Nz += 2;
-	//Nx = Ny = Nz;	// Cubic domain
 
 	int N = Nx*Ny*Nz;
 	int dist_mem_size = N*sizeof(double);
diff --git a/tests/lbpm_segmented_decomp.cpp b/tests/lbpm_segmented_decomp.cpp
index 2e30620d..63d4cb7e 100644
--- a/tests/lbpm_segmented_decomp.cpp
+++ b/tests/lbpm_segmented_decomp.cpp
@@ -33,6 +33,7 @@ int main(int argc, char **argv)
 	  printf("Solid Label: %i \n",SOLID);
 	  printf("NWP Label: %i \n",NWP);
 	}
+	{
     //.......................................................................
     // Reading the domain information file
     //.......................................................................
diff --git a/tests/lbpm_segmented_pp.cpp b/tests/lbpm_segmented_pp.cpp
index 0962874d..852e15f6 100644
--- a/tests/lbpm_segmented_pp.cpp
+++ b/tests/lbpm_segmented_pp.cpp
@@ -185,7 +185,7 @@ int main(int argc, char **argv)
     MPI_Comm comm = MPI_COMM_WORLD;
 	MPI_Comm_rank(comm,&rank);
 	MPI_Comm_size(comm,&nprocs);
-	
+	{	
     //.......................................................................
     // Reading the domain information file
     //.......................................................................
@@ -442,6 +442,7 @@ int main(int argc, char **argv)
 	Averages.SortBlobs();
 	Averages.PrintComponents(timestep);
 	*/
+	}
     MPI_Barrier(comm);
 	MPI_Finalize();
     return 0;
diff --git a/tests/lbpm_squaretube_pp.cpp b/tests/lbpm_squaretube_pp.cpp
index d4a9b794..d506f7c2 100644
--- a/tests/lbpm_squaretube_pp.cpp
+++ b/tests/lbpm_squaretube_pp.cpp
@@ -38,7 +38,9 @@ int main(int argc, char **argv)
 	MPI_Request req1[18],req2[18];
 	MPI_Status stat1[18],stat2[18];
 
-	int ORIENTATION;
+	int ORIENTATION=2; //default: the tube is aligned with Z axis
+	                   //ORIENTATION = 0: tube is aligned with X axis
+	                   //ORIENTATION = 1: tube is aligned with Y axis
 	double TubeWidth =15.0;
 	int BC;
 	int BubbleTop,BubbleBottom;
@@ -161,34 +163,66 @@ int main(int argc, char **argv)
 
 				if (ORIENTATION==0){
 					// square capillary tube aligned with the x direction
-					Averages.SDs(i,j,k) = TubeWidth/2 - fabs(i-0.5*Ny);
-					Averages.SDs(i,j,k) = min(Averages.SDs(i,j,k),TubeWidth/2-fabs(j-0.5*Nz));
+					Averages.SDs(i,j,k) = TubeWidth/2 - fabs(j-0.5*Ny);
+					Averages.SDs(i,j,k) = min(Averages.SDs(i,j,k),TubeWidth/2-fabs(k-0.5*Nz));
+					// Initialize phase positions
+					if (Averages.SDs(i,j,k) < 0.0){
+						id[n] = 0;
+					}
+					else if (Dm.iproc*Nx+k<BubbleBottom){
+						id[n] = 2;
+						sum++;
+					}
+					else if (Dm.iproc*Nx+k<BubbleTop){
+						id[n] = 1;
+						sum++;
+					}
+					else{
+						id[n] = 2;
+						sum++;
+					}
 				}
 				else if (ORIENTATION==1){
 					// square capillary tube aligned with the y direction
 					Averages.SDs(i,j,k) = TubeWidth/2 - fabs(i-0.5*Nx);
-					Averages.SDs(i,j,k) = min(Averages.SDs(i,j,k),TubeWidth/2-fabs(j-0.5*Nz));
+					Averages.SDs(i,j,k) = min(Averages.SDs(i,j,k),TubeWidth/2-fabs(k-0.5*Nz));
+					// Initialize phase positions
+					if (Averages.SDs(i,j,k) < 0.0){
+						id[n] = 0;
+					}
+					else if (Dm.jproc*Ny+k<BubbleBottom){
+						id[n] = 2;
+						sum++;
+					}
+					else if (Dm.jproc*Ny+k<BubbleTop){
+						id[n] = 1;
+						sum++;
+					}
+					else{
+						id[n] = 2;
+						sum++;
+					}
 				}
-				else {
+				else { //
 					// square capillary tube aligned with the z direction
 					Averages.SDs(i,j,k) = TubeWidth/2 - fabs(i-0.5*Nx);
 					Averages.SDs(i,j,k) = min(Averages.SDs(i,j,k),TubeWidth/2-fabs(j-0.5*Ny));
-				}
-				// Initialize phase positions
-				if (Averages.SDs(i,j,k) < 0.0){
-					id[n] = 0;
-				}
-				else if (Dm.kproc*Nz+k<BubbleBottom){
-					id[n] = 2;
-					sum++;
-				}
-				else if (Dm.kproc*Nz+k<BubbleTop){
-					id[n] = 1;
-					sum++;
-				}
-				else{
-					id[n] = 2;
-					sum++;
+					// Initialize phase positions
+					if (Averages.SDs(i,j,k) < 0.0){
+						id[n] = 0;
+					}
+					else if (Dm.kproc*Nz+k<BubbleBottom){
+						id[n] = 2;
+						sum++;
+					}
+					else if (Dm.kproc*Nz+k<BubbleTop){
+						id[n] = 1;
+						sum++;
+					}
+					else{
+						id[n] = 2;
+						sum++;
+					}
 				}
 			}
 		}
diff --git a/workflows/Preprocess/preprocess.py b/workflows/Preprocess/preprocess.py
new file mode 100644
index 00000000..088e6360
--- /dev/null
+++ b/workflows/Preprocess/preprocess.py
@@ -0,0 +1,123 @@
+#!/usr/bin/env python
+
+import numpy as np
+import sys
+import os
+import shutil
+from preprocess_utils import *
+
+# 1. this script take one input argument: experiment.csv
+# 2. Users need to put two things into this script
+#    2.1 surface tension in physical units
+#    2.2 physical depth
+# 3. Read csv + convert to LBM information
+# 4. write Color.in,-> this needs info from LBM pressure
+# 5. write Segmented.in -> this needs names of segmented image data
+# 6. write Domain.in -> this needs info on how users want to decompose the domain
+
+# Check if there is a proper command line argument
+if len(sys.argv) !=2:
+    sys.stderr.write('Usage: ' + sys.argv[0] + ' <Input experiment file>\n')
+    sys.exit()
+# end if
+
+#experiment_file = 'experiment.csv' # Give the name of the experiment file (e.g. *.csv)
+experiment_file = sys.argv[1] # Give the name of the experiment file (e.g. *.csv)
+seg_image_data_suffix = '_segmented.raw'# suffix of the segmented image data
+                                        # (should be consistent with what's in the segmenting script)
+                                        # TODO: It'd be good to find a better way to do this
+image_format = '.tiff'
+process_name = 'drainage'
+
+#### Users need to put information here ####
+ift = 24.0   # dyne/cm
+Depth = 8.8  # micron
+
+# A list of all default values in 'Color.in'
+#    Para['tau']         = 0.7
+#    Para['alpha']       = 0.005
+#    Para['beta']        = 0.95
+#    Para['phi_solid']   = -1.0
+#    Para['saturation']  = 0.0
+#    Para['Fx']          = 0.0
+#    Para['Fy']          = 0.0
+#    Para['Fz']          = 0.0
+#    Para['Restart']     = 0
+#    Para['pBC']         = 1
+#    Para['din']         = 1.001
+#    Para['dout']        = 0.999
+#    Para['maxtime']     = 100005
+#    Para['interval']    = 2000
+#    Para['tolerance']   = 1e-5
+
+# ***** Update any variables in 'Color.in', using names given in the key ***** # 
+alpha = 0.01
+
+# **************************************************************************** #
+
+# A list of all default values in 'Domain.in'
+#    Para['nprocx']       = 1
+#    Para['nprocy']       = 2
+#    Para['nprocz']       = 2
+#    Para['nx']           = 1
+#    Para['ny']           = 2
+#    Para['nz']           = 2
+#    Para['nspheres']     = 0 # deprecated
+#    Para['Lx']           = 10
+#    Para['Ly']           = 500
+#    Para['Lz']           = 500
+
+# ***** Update any variables in 'Domain.in', using names given in the key ***** # 
+
+
+# ***************************************************************************** #
+
+# A list of all default values in 'Segmented.in'
+#    Para['file_name']     = 'Micromodel_1_segmented.raw'
+#    Para['Nx']            = 10
+#    Para['Ny']            = 500
+#    Para['Nz']            = 500
+#    Para['xStart']        = 0
+#    Para['yStart']        = 0
+#    Para['zStart']        = 0
+
+# ***** Update any variables in 'Segmented.in', using names given in the key ***** # 
+
+
+# ******************************************************************************** #
+
+# Extract key parameters for LBM simulation from the experimental input *.csv file
+(Seg_data_name,din,dout)=get_LBM_parameters(experiment_file,seg_image_data_suffix,image_format,\
+                                            ift=ift,Depth=Depth)
+# Now 'name_for_Segmented_in' should match the name of segmented data files that are already generated
+
+# Write out 'Color.in', 'Domain.in' and 'Segmented.in' files
+cwd = os.getcwd()
+for k in range(Seg_data_name.size):
+    tag = k+1 # tag number for different folders to be created
+    dir_name = process_name+'_'+str(tag)
+    print "Creating folder : "+dir_name
+    if not os.path.exists(dir_name):
+        os.mkdir(dir_name)
+    #end if
+    
+    # Either move the corresponding '*.raw' data file to the folder 'dir_name'
+    #os.rename('./'+Seg_data_name[k],'./'+dir_name+'/'+Seg_data_name[k])
+    # Or copy the corresponding '*.raw' data file to the folder 'dir_name'
+    shutil.copy('./'+Seg_data_name[k],'./'+dir_name)
+    
+    # Change to the corresponding folder and write all input files 
+    os.chdir(dir_name)
+    write_Color_in_file(din=din[k],dout=dout[k],alpha=alpha)
+    write_Segment_in_file(Seg_data_name=Seg_data_name[k])
+    write_Domain_in_file()
+
+    os.chdir(cwd)
+#end for    
+
+
+
+
+
+
+
diff --git a/workflows/Preprocess/preprocess_utils.py b/workflows/Preprocess/preprocess_utils.py
new file mode 100644
index 00000000..da7d411a
--- /dev/null
+++ b/workflows/Preprocess/preprocess_utils.py
@@ -0,0 +1,272 @@
+#!/usr/bin/env python
+
+import os
+import numpy as np
+import csv
+
+def get_LBM_parameters(csv_file_name,base_name_suffix,image_format,ift=24,Depth=8.8,**kwargs):
+    # 'ift': dyne/cm
+    # 'Depth': micron
+    # Users need to provide the following information
+    # 1. surface tension in physical unit
+    # 2. physical depth  in physical unit
+    # 3. experimental file e.g. *.csv
+    # 4. Other optional information, which is summarised in a dictionary
+    #    'Para', including: 
+    #    Para['D']: characteristic length
+    #    Para['alpha']: LBM parameters, controlling the surface tension
+    #    Para['fitting']: LBM parameters, extracted from the bubble test
+
+
+    # Experimental definitions - units are converted in terms of N, m
+    micron=1e-6
+    cm=1e-2
+    dyne=1e-5
+    kPa=1e3
+    Para={'D':30.0} # characteristic length
+    #length=500*micron
+    
+    # LBM related parameters
+    # TODO: for parameters like 'alpha', you should find a way to 
+    #       communicate with other functions
+    #       It cannot be hard-coded here
+    Para['alpha'] = 0.01
+    Para['fitting'] = 5.796
+    # Check users' input arguments
+    for key in kwargs:
+        if key in Para.keys():
+            Para[key]=kwargs[key]
+        else:
+            print "Error: "+key+" is not a valid input !\n"
+            print "Error: LBM pressure boundaries are not set successfully !"
+            return
+        #end if
+    #end for
+    IFT = Para['alpha'] * Para['fitting']
+
+    # 'ift' is the surface tension from the experiment in physical unit
+    ift=ift*dyne/cm # This converts the unit of 'ift' to SI unit 
+    
+    # Process the experimental data
+    # 'EXP_data'  : original experimental data
+    #               It is a recorded numpy array, which means that its component
+    #               can be accessed by 'EXP_data.sw' or 'EXP_data['sw']'
+    # About 'EXP_data' see more information from the function 'read_csv'
+    EXP_data = read_csv(csv_file_name)
+    # A few notes for the experimental data
+    # 'Jwn': mean curvature in physical unit (e.g. 1/[micron])
+    # 'pwn': pressure difference in physical (e.g. kPa)
+
+    # Overall we need to map the measured physical pressures to get the principle radius of curvature R1 
+    # and associated curvature J1, and similarly R2 and J2 in the model's depth direction
+    # J1: principal curvature 1 
+    # J2: principal curvature 2 along the model's depth direction
+
+    # 'pc' is the dimensionless mean curvature scaled by the length scale 'D32'
+    # 'pc' is extracted from experimentally measured mean curvature 'Jwn'
+    pc = Para['D']*micron*EXP_data.Jwn*(1.0/micron)
+
+    # Alternatively, the dimensionless mean curvature can also be extracted from the 
+    # experimentally measured pressure difference (i.e. capillary pressure)
+    # 'pwn' is the dimensionless mean curvature scaled by the length scale 'D32' 
+    # pwn = Para['D']*micron*EXP_data.pwn*kPa/ift 
+
+    # Curvature is fixed in the micromodel "depth" direction
+    J2 = Para['D']*micron/(Depth*micron/2.0) 
+
+    # infer the curvature in the other direction
+    J1 = pc-J2
+
+    # determine the LBM pressure difference
+    dp=(J1+J2)*IFT/Para['D']
+    # determine the boundary pressure values for Color.in
+    din  = 1.0+0.5*dp 
+    dout = 1.0-0.5*dp
+    
+    # Generate names of segmented image data for Segmented.in file
+    # Need the input 'base_name_suffix' (e.g. '_segmented.raw')
+    data_name_for_Segmented_in = EXP_data.image_name.copy()
+    data_name_for_Segmented_in = data_name_for_Segmented_in.replace(image_format,base_name_suffix)
+
+    return (data_name_for_Segmented_in,din,dout)
+
+#end def
+
+def write_Domain_in_file(**kwargs):
+    # For most of the parameters in the Color.in file
+    # you wouldn't need to worry about them as by default 
+    # they are all hard-coded here 
+    Para={'nprocx':1}
+    Para['nprocy']=2
+    Para['nprocz']=2
+    Para['nx']=1
+    Para['ny']=2
+    Para['nz']=2
+    Para['nspheres']=0 # deprecated
+    Para['Lx']=10
+    Para['Ly']=500
+    Para['Lz']=500
+
+    for key in kwargs:
+        if key in Para.keys():
+            Para[key]=kwargs[key]
+        else:
+            print "Error: "+key+" is not a valid input !\n"
+            print "Error: Domain.in file is not writen."
+            return
+        #end if
+    #end for
+    
+    # Check if the dompositoin requirement is satisfied
+    if (Para['nx']*Para['nprocx']>Para['Lx']) or \
+       (Para['ny']*Para['nprocy']>Para['Ly']) or \
+       (Para['nz']*Para['nprocz']>Para['Lz']):
+            print "Error: The decomposed size does not match the total domain size!"
+            return
+    #end if
+
+    # write the Domain.in file
+    ParamFile = open("Domain.in","w")
+    ParamFile.write("%i "   % Para['nprocx'])
+    ParamFile.write("%i "   % Para['nprocy'])
+    ParamFile.write("%i\n" % Para['nprocz'])
+    ParamFile.write("%i "   % Para['nx'])
+    ParamFile.write("%i "   % Para['ny'])
+    ParamFile.write("%i\n" % Para['nz'])
+    ParamFile.write("%i\n" % Para['nspheres'])
+    ParamFile.write("%.1f "   % Para['Lx'])
+    ParamFile.write("%.1f "   % Para['Ly'])
+    ParamFile.write("%.1f\n" % Para['Lz'])
+    ParamFile.close()
+    print "A Domain.in file is created."
+#end def
+
+
+def write_Segment_in_file(**kwargs):
+    # For most of the parameters in the Color.in file
+    # you wouldn't need to worry about them as by default 
+    # they are all hard-coded here 
+    Para={'Seg_data_name':'Micromodel_1_segmented.raw'}
+    Para['Nx']=10
+    Para['Ny']=500
+    Para['Nz']=500
+    Para['xStart']=0
+    Para['yStart']=0
+    Para['zStart']=0
+
+    for key in kwargs:
+        if key in Para.keys():
+            Para[key]=kwargs[key]
+        else:
+            print "Error: "+key+" is not a valid input !\n"
+            print "Error: Segmented.in file is not writen."
+            return
+        #end if
+    #end for
+
+    ParamFile = open("Segmented.in","w")
+    ParamFile.write("%s\n" % Para['Seg_data_name'])
+    ParamFile.write("%i "  % Para['Nx'])
+    ParamFile.write("%i "  % Para['Ny'])
+    ParamFile.write("%i\n" % Para['Nz'])
+    ParamFile.write("%i "  % Para['xStart'])
+    ParamFile.write("%i "  % Para['yStart'])
+    ParamFile.write("%i\n" % Para['zStart'])
+    ParamFile.close()
+    print "A Segmented.in file is created."
+#end def
+
+def write_Color_in_file(**kwargs):
+    # For most of the parameters in the Color.in file
+    # you wouldn't need to worry about them as by default 
+    # they are all hard-coded here 
+    Para={'tau':0.7}
+    Para['alpha']=0.01
+    Para['beta']=0.95
+    Para['phi_solid']=-1.0
+    Para['saturation']=0.0
+    Para['Fx']=0.0
+    Para['Fy']=0.0
+    Para['Fz']=0.0
+    Para['Restart']=0
+    Para['pBC']=1
+    Para['din']=1.001
+    Para['dout']=0.999
+    Para['maxtime']=100005
+    Para['interval']=2000
+    Para['tolerance']=1e-5
+
+    for key in kwargs:
+        if key in Para.keys():
+            Para[key]=kwargs[key]
+        else:
+            print "Error: "+key+" is not a valid input !\n"
+            print "Error: Color.in file is not writen."
+            return
+        #end if
+    #end for
+    
+    # write the color.in file
+    ParamFile = open("Color.in","w")
+    ParamFile.write("%f\n" % Para['tau'])
+    ParamFile.write("%f " % Para['alpha'])
+    ParamFile.write("%f " % Para['beta'])
+    ParamFile.write("%f\n" % Para['phi_solid'])
+    ParamFile.write("%f\n" % Para['saturation'])
+    ParamFile.write("%f " % Para['Fx'])
+    ParamFile.write("%f " % Para['Fy'])
+    ParamFile.write("%f\n" % Para['Fz'])
+    ParamFile.write("%i " % Para['Restart'])
+    ParamFile.write("%i " % Para['pBC'])
+    ParamFile.write("%f " % Para['din'])
+    ParamFile.write("%f\n" % Para['dout'])
+    ParamFile.write("%i " % Para['maxtime'])
+    ParamFile.write("%i " % Para['interval'])
+    ParamFile.write("%e\n" % Para['tolerance'])
+    ParamFile.close()
+    print "A Color.in file is created."
+#end def
+
+def read_csv(csv_file_name):
+    #TODO Haven't thought about a better way of doing this
+    # Right now I just hard-code the possible variables from the experiment
+    # which means users are forced to prepare their *.csv file as listed here
+    image_name = []
+    sw = []
+    Jwn = []
+    pwn = []
+    with open(csv_file_name,"r") as f:
+        for line in f:
+            reader=csv.reader(f,delimiter=' ') #Note the header is skipped by default 
+            for row in reader:
+                image_name.append(row[0])
+                sw.append(float(row[1]))
+                Jwn.append(float(row[2]))
+                pwn.append(float(row[3]))
+            #end for
+        #end for
+    #end with
+    
+    # Convter the list to numpy array
+    image_name_array = np.asarray(image_name,dtype=str)
+    sw_array  = np.asarray(sw,dtype=np.float32)
+    Jwn_array = np.asarray(Jwn,dtype=np.float32)
+    pwn_array = np.asarray(pwn,dtype=np.float32)
+
+    # Restore the original shape of the experimental data
+    experiment_data = np.column_stack((image_name_array,sw_array,Jwn_array,pwn_array))
+    # Unfortunately the column stack will cast all float data to strings
+    experiment_data = np.core.records.fromrecords(experiment_data,names='image_name,sw,Jwn,pwn')
+    dt=experiment_data.dtype.descr
+    dt[1] = (dt[1][0],'float32')
+    dt[2] = (dt[2][0],'float32')
+    dt[3] = (dt[3][0],'float32')
+    experiment_data = experiment_data.astype(dt)
+    return experiment_data
+#end def
+
+
+
+
+
+
diff --git a/workflows/SegmentingImage/SegmentMicromodelTiff_Group.py b/workflows/SegmentingImage/SegmentMicromodelTiff_Group.py
new file mode 100644
index 00000000..b0d8b701
--- /dev/null
+++ b/workflows/SegmentingImage/SegmentMicromodelTiff_Group.py
@@ -0,0 +1,104 @@
+#!/usr/bin/env python
+
+import numpy as np
+import matplotlib.pylab as plt
+from glob import glob
+from PIL import Image # import Python Imaging Library (PIL)
+import sys
+from SegmentMicromodelTiff_utils import convert_image_to_array,read_input_parameters
+
+# Check if there is a proper command line argument
+if len(sys.argv) !=2:
+    sys.stderr.write('Usage: ' + sys.argv[0] + ' <Input parameter file>\n')
+    sys.exit()
+# end if
+
+input_parameter_file = sys.argv[1] # Give the name of the input parameter file 
+base_name_output = '_segmented.raw'# The string will be appended at the end of the 
+                                   # name of the input image as an output file name
+image_format = '.tiff' # format of experimental images
+
+# Read other input parameters
+Para = read_input_parameters(input_parameter_file)
+imin = Para['imin']
+imax = Para['imax']
+imin_b = Para['imin_b'] # greater than 'imin'
+imax_b = Para['imax_b'] # smaller than 'imax'
+top_layer = Para['top_layer']
+bot_layer = Para['bot_layer']
+DEPTH = Para['DEPTH']
+NY = imax - imin
+NZ = NY
+NX = DEPTH+top_layer+bot_layer
+
+# parameters for segmentation
+threshold_nw = Para['threshold_nw'] # NW phase: RGB values > threshold_nw
+threshold_s  = Para['threshold_s']  # solid:    RGB values < threshold_s
+# W phase: threshold_s <= RGB values <= threshold_nw
+
+# Load a group of image files with 'base_name' in the names of files 
+# e.g. 'Micromodel_1.tiff' etc.
+file_input_group = glob('*'+Para['base_name']+'*'+image_format) # need to match the image format
+
+# Process all imported experimental images
+if not file_input_group:
+    print 'Error: Input files cannot be found ! '
+else:
+    for ii in range(len(file_input_group)):
+        file_input_single = file_input_group[ii]
+        print "Processing image "+file_input_single+" now..."
+        print "------ Micromodel dimensions (NX, NY, NZ) are (%i, %i, %i) ------" % (NX,NY,NZ)
+        
+        # Get an array from the input .tiff image
+        im_array = convert_image_to_array(file_input_single,imin,imax)
+
+        # Initialize the segmented image 'ID'
+        # NOTE: 'ID' has the dimension (height, width, DEPTH)
+        ID = np.zeros((NZ,NY,NX),dtype=np.uint8)
+        # Map the picels to the 3D geometry
+        # 1. Identify the non-wetting phase
+        ID[im_array>threshold_nw,top_layer:top_layer+DEPTH] = 1
+        # 2. Identify the wetting phase
+        ID[np.logical_and(im_array>=threshold_s,im_array<=threshold_nw),\
+           top_layer:top_layer+DEPTH] = 2
+        # 3. Post boundary retreatment along y-axis. Note z-axis is always the flow direction
+        ID[:,:imin_b-imin,top_layer:top_layer+DEPTH]=0
+        ID[:,ID.shape[1]-(imax-imax_b):,top_layer:top_layer+DEPTH]=0
+
+        # Calculate the porosity
+        POROSITY = 1.0 - 1.0*ID[ID==0].size/ID.size
+        print "Porosity of the micromodel is: "+str(POROSITY) 
+        
+        # Write the segmeted data to the output
+        # set the output file names (need to chop off say the '.tiff' part)
+        output_file_name = file_input_single[:file_input_single.find('.')]+base_name_output
+        # Write out the segmented data (in binary format)
+        print "The segmented image is written to "+output_file_name+" now..."
+        ID.tofile(output_file_name)
+        # NOTE when you want to load the binary file
+        # Do as follows:
+        # ID_reload = np.fromfile(file_name,dtype=np.uint8)
+
+    #end for
+#end if
+
+    
+# In case you want to test the segmentation
+#plt.figure(1)
+#plt.subplot(1,2,1)
+#plt.title('original RGB figure')    
+#plt.pcolormesh(im_array);
+#plt.axis('equal')
+#plt.colorbar()
+#plt.subplot(1,2,2)
+#plt.title('Segmented image')
+# This will show the last-processed segmented image
+#cmap = plt.cm.get_cmap('hot',3) #Plot 3 discrete colors for NW, W and Solids
+#cax=plt.pcolormesh(ID[:,:,NX/2],cmap=cmap,vmin=-0.5,vmax=2.5);
+#cbar = plt.colorbar(cax,ticks=[0,1,2])
+#plt.axis('equal')
+#plt.show()
+
+
+
+
diff --git a/workflows/SegmentingImage/SegmentMicromodelTiff_utils.py b/workflows/SegmentingImage/SegmentMicromodelTiff_utils.py
new file mode 100644
index 00000000..92cf6662
--- /dev/null
+++ b/workflows/SegmentingImage/SegmentMicromodelTiff_utils.py
@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+
+import numpy as np
+from PIL import Image # import Python Imaging Library (PIL)
+
+def convert_image_to_array(file_name,imin,imax):
+    # file_name: the name of the input image (a string)
+    # imin: (left) boundary that defines the output array
+    # imax: (right) boundary that defines the output array
+    # Note: assume a grayscale image (whose R, G, and B values are the same)
+    # Return an numpy array with RGB values that has the same size as the wanted micromodel
+    im = Image.open(file_name).convert('RGB')
+    im_array = np.array(im) # 'im_array' has dimension (height, width, 3)
+    # Grayscale input image is assumed, thus its R, G, and B values are the same
+    im_array = im_array[:,:,0] # 2D array with dimension (height, width)
+    im_array = im_array[:,imin:imax] # Chop the original image to a desirable size
+    #im_array = np.flipud(im_array[:,imin:imax]) # You might want to flip the Z-axis
+    return im_array
+#end def
+
+def read_input_parameters(input_file_name):
+    # input_file_name is a string
+    # The *.in file has the following lines
+    # line 0: the base name of the experimental images
+    # line 1: imin imax
+    # line 2: imin_b imax_b
+    # line 3: top_layer bottom_layer DEPTH
+    # line 4: threshold_nw threshold_s
+    # Note:  'threshold_nw' means: NW phase: RGB values > threshold_nw
+    #        'threshold_s'  means: solid:    RGB values < threshold_s
+    f = open(input_file_name,'r') # read-only
+    lines = f.readlines()
+    output_file = {'base_name':lines[0].splitlines()[0]}
+    line1_array = np.fromstring(lines[1].splitlines()[0],dtype=np.int32,sep=' ')
+    line2_array = np.fromstring(lines[2].splitlines()[0],dtype=np.int32,sep=' ')
+    line3_array = np.fromstring(lines[3].splitlines()[0],dtype=np.int32,sep=' ')
+    line4_array = np.fromstring(lines[4].splitlines()[0],dtype=np.int32,sep=' ')
+    output_file['imin']=line1_array[0]
+    output_file['imax']=line1_array[1]
+    output_file['imin_b']=line2_array[0]
+    output_file['imax_b']=line2_array[1]
+    output_file['top_layer']=line3_array[0]
+    output_file['bot_layer']=line3_array[1]
+    output_file['DEPTH']=line3_array[2]
+    output_file['threshold_nw']=line4_array[0]
+    output_file['threshold_s'] =line4_array[1]
+    f.close()
+    return output_file
+#end def
+
+#def get_segmented_array(image_array,parameters,NX,NY,NZ):
+#    # 'image_array' is a numpy array of RGB values, with the same size as micromodel 
+#    # 'parameters' is a dictionary of input parameters
+#    
+#    # Initialize the segmented image 'ID'
+#    # NOTE: 'ID' has the dimension (DEPTH, height, width)
+#    ID = np.zeros((NX,NZ,NY),dtype=np.uint8)
+#    # Map the picels to the 3D geometry
+#    # 1. Identify the non-wetting phase
+#    ID[top_layer:top_layer+DEPTH,im_array>threshold_nw] = 1
+#    # 2. Identify the wetting phase
+#    ID[top_layer:top_layer+DEPTH,\
+#    np.logical_and(im_array>=threshold_s,im_array<=threshold_nw)] = 2
+#    # 3. Post boundary retreatment along y-axis
+#    ID[top_layer:top_layer+DEPTH,:,:imin_b-imin]=0
+#    ID[top_layer:top_layer+DEPTH,:,ID.shape[2]-(imax-imax_b):]=0
+#    return ID
+##end def