diff --git a/CMakeLists.txt b/CMakeLists.txt
index 31518d86..c9b7e66a 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -18,8 +18,7 @@ SET( TEST_MAX_PROCS 16 )
 
 
 # Initialize the project
-PROJECT( ${PROJ} LANGUAGES CXX )
-
+PROJECT( ${PROJ} LANGUAGES CXX)
 
 # Prevent users from building in place
 IF ("${CMAKE_CURRENT_SOURCE_DIR}" STREQUAL "${CMAKE_CURRENT_BINARY_DIR}" )
diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp
index 2a0b7169..dc932f91 100644
--- a/analysis/ElectroChemistry.cpp
+++ b/analysis/ElectroChemistry.cpp
@@ -53,19 +53,32 @@ void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poiss
 	Poisson.getElectricPotential(ElectricalPotential);
 	
 	/* local sub-domain averages */
-	double rho_avg_local[Ion.number_ion_species];
-	double rho_mu_avg_local[Ion.number_ion_species];
-	double rho_mu_fluctuation_local[Ion.number_ion_species];	
-	double rho_psi_avg_local[Ion.number_ion_species];
-	double rho_psi_fluctuation_local[Ion.number_ion_species];
+	double *rho_avg_local;
+	double *rho_mu_avg_local;
+	double *rho_mu_fluctuation_local;	
+	double *rho_psi_avg_local;
+	double *rho_psi_fluctuation_local;
 	/* global averages */
-	double rho_avg_global[Ion.number_ion_species];
-	double rho_mu_avg_global[Ion.number_ion_species];
-	double rho_mu_fluctuation_global[Ion.number_ion_species];
-	double rho_psi_avg_global[Ion.number_ion_species];
-	double rho_psi_fluctuation_global[Ion.number_ion_species];
+	double *rho_avg_global;
+	double *rho_mu_avg_global;
+	double *rho_mu_fluctuation_global;
+	double *rho_psi_avg_global;
+	double *rho_psi_fluctuation_global;
 	
-	for (int ion=0; ion<Ion.number_ion_species; ion++){
+	/* local sub-domain averages */
+	rho_avg_local = new double [Ion.number_ion_species];
+	rho_mu_avg_local = new double [Ion.number_ion_species];
+	rho_mu_fluctuation_local = new double [Ion.number_ion_species];	
+	rho_psi_avg_local = new double [Ion.number_ion_species];
+	rho_psi_fluctuation_local = new double [Ion.number_ion_species];
+	/* global averages */
+	rho_avg_global = new double [Ion.number_ion_species];
+	rho_mu_avg_global = new double [Ion.number_ion_species];
+	rho_mu_fluctuation_global = new double [Ion.number_ion_species];
+	rho_psi_avg_global = new double [Ion.number_ion_species];
+	rho_psi_fluctuation_global = new double [Ion.number_ion_species];
+	
+	for (size_t ion=0; ion<Ion.number_ion_species; ion++){
 		rho_avg_local[ion] = 0.0;
 		rho_mu_avg_local[ion] = 0.0;
 		rho_psi_avg_local[ion] = 0.0;
@@ -90,7 +103,7 @@ void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poiss
 		}
 	}
 	
-	for (int ion=0; ion<Ion.number_ion_species; ion++){
+	for (size_t ion=0; ion<Ion.number_ion_species; ion++){
 		rho_mu_fluctuation_local[ion] = 0.0;
 		rho_psi_fluctuation_local[ion] = 0.0;
 		/* Compute averages for each ion */
@@ -108,7 +121,7 @@ void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poiss
 	
 	if (Dm->rank()==0){	
 		fprintf(TIMELOG,"%i ",timestep); 
-		for (int ion=0; ion<Ion.number_ion_species; ion++){
+		for (size_t ion=0; ion<Ion.number_ion_species; ion++){
 			fprintf(TIMELOG,"%.8g ",rho_avg_global[ion]);
 			fprintf(TIMELOG,"%.8g ",rho_mu_avg_global[ion]);
 			fprintf(TIMELOG,"%.8g ",rho_psi_avg_global[ion]);
@@ -147,7 +160,7 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
     visData[0].mesh = std::make_shared<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
     auto ElectricPotential = std::make_shared<IO::Variable>();
     std::vector<shared_ptr<IO::Variable>> IonConcentration;
-    for (int ion=0; ion<Ion.number_ion_species; ion++){
+    for (size_t ion=0; ion<Ion.number_ion_species; ion++){
         IonConcentration.push_back(std::make_shared<IO::Variable>());
     }
     auto VxVar = std::make_shared<IO::Variable>();
@@ -163,8 +176,8 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
     }
 
     if (vis_db->getWithDefault<bool>( "save_concentration", true )){
-    	for (int ion=0; ion<Ion.number_ion_species; ion++){
-    		sprintf(VisName,"IonConcentration_%i",ion+1);
+    	for (size_t ion=0; ion<Ion.number_ion_species; ion++){
+    		sprintf(VisName,"IonConcentration_%zu",ion+1);
     		IonConcentration[ion]->name = VisName;
     		IonConcentration[ion]->type = IO::VariableType::VolumeVariable;
     		IonConcentration[ion]->dim = 1;
@@ -199,8 +212,8 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
     }
 
     if (vis_db->getWithDefault<bool>( "save_concentration", true )){
-    	for (int ion=0; ion<Ion.number_ion_species; ion++){
-    		sprintf(VisName,"IonConcentration_%i",ion+1);
+    	for (size_t ion=0; ion<Ion.number_ion_species; ion++){
+    		sprintf(VisName,"IonConcentration_%zu",ion+1);
     		IonConcentration[ion]->name = VisName;
     		ASSERT(visData[0].vars[1+ion]->name==VisName);
     		Array<double>& IonConcentrationData = visData[0].vars[1+ion]->data;
diff --git a/analysis/FreeEnergy.cpp b/analysis/FreeEnergy.cpp
index 6a641a95..b49f2ec5 100644
--- a/analysis/FreeEnergy.cpp
+++ b/analysis/FreeEnergy.cpp
@@ -47,8 +47,6 @@ void FreeEnergyAnalyzer::SetParams(){
 
 void FreeEnergyAnalyzer::Basic(ScaLBL_FreeLeeModel &LeeModel, int timestep){
 
-	int i,j,k;
-
 	if (Dm->rank()==0){	
 		fprintf(TIMELOG,"%i ",timestep); 
 		/*for (int ion=0; ion<Ion.number_ion_species; ion++){
@@ -78,7 +76,6 @@ void FreeEnergyAnalyzer::Basic(ScaLBL_FreeLeeModel &LeeModel, int timestep){
 void FreeEnergyAnalyzer::WriteVis( ScaLBL_FreeLeeModel &LeeModel, std::shared_ptr<Database> input_db, int timestep){
 	
 	auto vis_db =  input_db->getDatabase( "Visualization" );
-    char VisName[40];	
     
     std::vector<IO::MeshDataStruct> visData;
 	fillHalo<double> fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1);
diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp
index 324c45cb..4b023db6 100644
--- a/analysis/morphology.cpp
+++ b/analysis/morphology.cpp
@@ -120,6 +120,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
 	sendtag = recvtag = 7;
 
 	int ii,jj,kk;
+	int imin,jmin,kmin,imax,jmax,kmax;
 	int Nx = nx;
 	int Ny = ny;
 	int Nz = nz;
@@ -128,17 +129,13 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
 	double void_fraction_new=1.0; 
 	double void_fraction_diff_old = 1.0;
 	double void_fraction_diff_new = 1.0;
-
-	// Increase the critical radius until the target saturation is met
-	double deltaR=0.05; // amount to change the radius in voxel units
-	double Rcrit_old;
-
-	int imin,jmin,kmin,imax,jmax,kmax;
-
 	if (ErodeLabel == 1){
 		VoidFraction = 1.0 - VoidFraction;
 	}
 
+	// Increase the critical radius until the target saturation is met
+	double deltaR=0.05; // amount to change the radius in voxel units
+	double Rcrit_old = maxdistGlobal;
 	double Rcrit_new = maxdistGlobal;
 
 	while (void_fraction_new > VoidFraction)
@@ -406,6 +403,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
 	sendtag = recvtag = 7;
 */
 	int ii,jj,kk;
+	int imin,jmin,kmin,imax,jmax,kmax;
 	int Nx = nx;
 	int Ny = ny;
 	int Nz = nz;
@@ -417,10 +415,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
 
 	// Increase the critical radius until the target saturation is met
 	double deltaR=0.05; // amount to change the radius in voxel units
-	double Rcrit_old;
-
-	int imin,jmin,kmin,imax,jmax,kmax;
-
+	double Rcrit_old = maxdistGlobal;
 	double Rcrit_new = maxdistGlobal;
 	//if (argc>2){
 	//	Rcrit_new = strtod(argv[2],NULL);
@@ -457,7 +452,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
 						for (kk=kmin; kk<kmax; kk++){
 							for (jj=jmin; jj<jmax; jj++){
 								for (ii=imin; ii<imax; ii++){
-									int nn = kk*nx*ny+jj*nx+ii;
 									double dsq = double((ii-i)*(ii-i)+(jj-j)*(jj-j)+(kk-k)*(kk-k));
 									if (ID(ii,jj,kk) == 2 && dsq <= (Rcrit_new+1)*(Rcrit_new+1)){
 										LocalNumber+=1.0;
@@ -578,7 +572,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
 						// nwp
 						phase(i,j,k) = -1.0;
 					}
-					else{
+					else{i
 						// treat solid as WP since films can connect 
 						phase(i,j,k) = 1.0;
 					}
diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index a3fc5fa7..09b709f6 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -745,7 +745,6 @@ runAnalysis::runAnalysis(  ScaLBL_ColorModel &ColorModel)
 	
 	d_comm = ColorModel.Dm->Comm.dup();
 	d_Np = ColorModel.Np;
-	bool Regular = false;
 	
 	auto input_db = ColorModel.db;
     auto db     = input_db->getDatabase( "Analysis" );
diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp
index fa3ecc74..dcb99a47 100644
--- a/common/ScaLBL.cpp
+++ b/common/ScaLBL.cpp
@@ -1062,29 +1062,28 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
 			}
 		}
 	}
-	
+
 	int *bb_dist_tmp = new int [local_count];	
 	int *bb_interactions_tmp = new int [local_count];	
 	ScaLBL_AllocateDeviceMemory((void **) &bb_dist, sizeof(int)*local_count);
 	ScaLBL_AllocateDeviceMemory((void **) &bb_interactions, sizeof(int)*local_count);
-    int *fluid_boundary_tmp;
-    double *lattice_weight_tmp;
-    float *lattice_cx_tmp;
-    float *lattice_cy_tmp;
-    float *lattice_cz_tmp;
-    if(SlippingVelBC==true){
-        fluid_boundary_tmp = new int [local_count];
-        lattice_weight_tmp = new double [local_count];
-        lattice_cx_tmp = new float [local_count];
-        lattice_cy_tmp = new float [local_count];
-        lattice_cz_tmp = new float [local_count];
-	    ScaLBL_AllocateDeviceMemory((void **) &fluid_boundary, sizeof(int)*local_count);
-	    ScaLBL_AllocateDeviceMemory((void **) &lattice_weight, sizeof(double)*local_count);
-	    ScaLBL_AllocateDeviceMemory((void **) &lattice_cx, sizeof(float)*local_count);
-	    ScaLBL_AllocateDeviceMemory((void **) &lattice_cy, sizeof(float)*local_count);
-	    ScaLBL_AllocateDeviceMemory((void **) &lattice_cz, sizeof(float)*local_count);
-    }
-	
+	int *fluid_boundary_tmp;
+	double *lattice_weight_tmp;
+	float *lattice_cx_tmp;
+	float *lattice_cy_tmp;
+	float *lattice_cz_tmp;
+	/* allocate memory for bounce-back sites */
+	fluid_boundary_tmp = new int [local_count];
+	lattice_weight_tmp = new double [local_count];
+	lattice_cx_tmp = new float [local_count];
+	lattice_cy_tmp = new float [local_count];
+	lattice_cz_tmp = new float [local_count];
+	ScaLBL_AllocateDeviceMemory((void **) &fluid_boundary, sizeof(int)*local_count);
+	ScaLBL_AllocateDeviceMemory((void **) &lattice_weight, sizeof(double)*local_count);
+	ScaLBL_AllocateDeviceMemory((void **) &lattice_cx, sizeof(float)*local_count);
+	ScaLBL_AllocateDeviceMemory((void **) &lattice_cy, sizeof(float)*local_count);
+	ScaLBL_AllocateDeviceMemory((void **) &lattice_cz, sizeof(float)*local_count);
+
 	local_count=0;
 	for (k=1;k<Nz-1;k++){
 		for (j=1;j<Ny-1;j++){
@@ -1097,78 +1096,78 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
 					neighbor=Map(i-1,j,k);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i-1) + (j)*Nx + (k)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/18.0;
                             lattice_cx_tmp[local_count] = -1.0;
                             lattice_cy_tmp[local_count] =  0.0;
                             lattice_cz_tmp[local_count] =  0.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 2*Np;
 					}
 
 					neighbor=Map(i+1,j,k);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i+1) + (j)*Nx + (k)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/18.0;
                             lattice_cx_tmp[local_count] =  1.0;
                             lattice_cy_tmp[local_count] =  0.0;
                             lattice_cz_tmp[local_count] =  0.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++] = idx + 1*Np;
 					}
 
 					neighbor=Map(i,j-1,k);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i) + (j-1)*Nx + (k)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/18.0;
                             lattice_cx_tmp[local_count] =  0.0;
                             lattice_cy_tmp[local_count] = -1.0;
                             lattice_cz_tmp[local_count] =  0.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 4*Np;
 					}
 
 					neighbor=Map(i,j+1,k);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i) + (j+1)*Nx + (k)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/18.0;
                             lattice_cx_tmp[local_count] =  0.0;
                             lattice_cy_tmp[local_count] =  1.0;
                             lattice_cz_tmp[local_count] =  0.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 3*Np;
 					}
 
 					neighbor=Map(i,j,k-1);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i) + (j)*Nx + (k-1)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/18.0;
                             lattice_cx_tmp[local_count] =  0.0;
                             lattice_cy_tmp[local_count] =  0.0;
                             lattice_cz_tmp[local_count] = -1.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 6*Np;
 					}
 
 					neighbor=Map(i,j,k+1);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i) + (j)*Nx + (k+1)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/18.0;
                             lattice_cx_tmp[local_count] =  0.0;
                             lattice_cy_tmp[local_count] =  0.0;
                             lattice_cz_tmp[local_count] =  1.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 5*Np;
 					}
 				}
@@ -1186,156 +1185,156 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
 					neighbor=Map(i-1,j-1,k);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i-1) + (j-1)*Nx + (k)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] = -1.0;
                             lattice_cy_tmp[local_count] = -1.0;
                             lattice_cz_tmp[local_count] =  0.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 8*Np;
 					}
 
 					neighbor=Map(i+1,j+1,k);
 					if (neighbor==-1)	{
 						bb_interactions_tmp[local_count] = (i+1) + (j+1)*Nx + (k)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] =  1.0;
                             lattice_cy_tmp[local_count] =  1.0;
                             lattice_cz_tmp[local_count] =  0.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 7*Np;
 					}
 
 					neighbor=Map(i-1,j+1,k);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i-1) + (j+1)*Nx + (k)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] = -1.0;
                             lattice_cy_tmp[local_count] =  1.0;
                             lattice_cz_tmp[local_count] =  0.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 10*Np;
 					}
 
 					neighbor=Map(i+1,j-1,k);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i+1) + (j-1)*Nx + (k)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] =  1.0;
                             lattice_cy_tmp[local_count] = -1.0;
                             lattice_cz_tmp[local_count] =  0.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 9*Np;
 					}
 
 					neighbor=Map(i-1,j,k-1);
 					if (neighbor==-1) {
 						bb_interactions_tmp[local_count] = (i-1) + (j)*Nx + (k-1)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] = -1.0;
                             lattice_cy_tmp[local_count] =  0.0;
                             lattice_cz_tmp[local_count] = -1.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 12*Np;
 					}
 
 					neighbor=Map(i+1,j,k+1);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i+1) + (j)*Nx + (k+1)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] =  1.0;
                             lattice_cy_tmp[local_count] =  0.0;
                             lattice_cz_tmp[local_count] =  1.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 11*Np;
 					}
 
 					neighbor=Map(i-1,j,k+1);
 					if (neighbor==-1) {
 						bb_interactions_tmp[local_count] = (i-1) + (j)*Nx + (k+1)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] = -1.0;
                             lattice_cy_tmp[local_count] =  0.0;
                             lattice_cz_tmp[local_count] =  1.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 14*Np;
 					}
 
 					neighbor=Map(i+1,j,k-1);
 					if (neighbor==-1) {
 						bb_interactions_tmp[local_count] = (i+1) + (j)*Nx + (k-1)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] =  1.0;
                             lattice_cy_tmp[local_count] =  0.0;
                             lattice_cz_tmp[local_count] = -1.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 13*Np;
 					}
 
 					neighbor=Map(i,j-1,k-1);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i) + (j-1)*Nx + (k-1)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] =  0.0;
                             lattice_cy_tmp[local_count] = -1.0;
                             lattice_cz_tmp[local_count] = -1.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 16*Np;
 					}
 
 					neighbor=Map(i,j+1,k+1);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i) + (j+1)*Nx + (k+1)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] =  0.0;
                             lattice_cy_tmp[local_count] =  1.0;
                             lattice_cz_tmp[local_count] =  1.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 15*Np;
 					}
 
 					neighbor=Map(i,j-1,k+1);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i) + (j-1)*Nx + (k+1)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] =  0.0;
                             lattice_cy_tmp[local_count] = -1.0;
                             lattice_cz_tmp[local_count] =  1.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 18*Np;
 					}
 
 					neighbor=Map(i,j+1,k-1);
 					if (neighbor==-1){
 						bb_interactions_tmp[local_count] = (i) + (j+1)*Nx + (k-1)*Nx*Ny;
-                        if(SlippingVelBC==true){
+                        //if(SlippingVelBC==true){
                             fluid_boundary_tmp[local_count] = idx;
                             lattice_weight_tmp[local_count] = 1.0/36.0;
                             lattice_cx_tmp[local_count] =  0.0;
                             lattice_cy_tmp[local_count] =  1.0;
                             lattice_cz_tmp[local_count] = -1.0;
-                        } 
+                        //} 
 						bb_dist_tmp[local_count++]=idx + 17*Np;
 					}
 				}
@@ -1345,24 +1344,20 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
 	n_bb_d3q19 = local_count; // this gives the d3q19 distributions not part of d3q7 model
 	ScaLBL_CopyToDevice(bb_dist, bb_dist_tmp, local_count*sizeof(int));
 	ScaLBL_CopyToDevice(bb_interactions, bb_interactions_tmp, local_count*sizeof(int));
-    if(SlippingVelBC==true){
-	    ScaLBL_CopyToDevice(fluid_boundary, fluid_boundary_tmp, local_count*sizeof(int));
-	    ScaLBL_CopyToDevice(lattice_weight, lattice_weight_tmp, local_count*sizeof(double));
-	    ScaLBL_CopyToDevice(lattice_cx, lattice_cx_tmp, local_count*sizeof(float));
-	    ScaLBL_CopyToDevice(lattice_cy, lattice_cy_tmp, local_count*sizeof(float));
-	    ScaLBL_CopyToDevice(lattice_cz, lattice_cz_tmp, local_count*sizeof(float));
-    }
+	ScaLBL_CopyToDevice(fluid_boundary, fluid_boundary_tmp, local_count*sizeof(int));
+	ScaLBL_CopyToDevice(lattice_weight, lattice_weight_tmp, local_count*sizeof(double));
+	ScaLBL_CopyToDevice(lattice_cx, lattice_cx_tmp, local_count*sizeof(float));
+	ScaLBL_CopyToDevice(lattice_cy, lattice_cy_tmp, local_count*sizeof(float));
+	ScaLBL_CopyToDevice(lattice_cz, lattice_cz_tmp, local_count*sizeof(float));
 	ScaLBL_DeviceBarrier();
-    
+
 	delete [] bb_dist_tmp;
 	delete [] bb_interactions_tmp;
-    if(SlippingVelBC==true){
-	    delete [] fluid_boundary_tmp;
-	    delete [] lattice_weight_tmp;
-        delete [] lattice_cx_tmp;
-        delete [] lattice_cy_tmp;
-        delete [] lattice_cz_tmp;
-    }
+	delete [] fluid_boundary_tmp;
+	delete [] lattice_weight_tmp;
+	delete [] lattice_cx_tmp;
+	delete [] lattice_cy_tmp;
+	delete [] lattice_cz_tmp;
 }
 
 void ScaLBL_Communicator::SolidDirichletD3Q7(double *fq, double *BoundaryValue){
diff --git a/cpu/BGK.cpp b/cpu/BGK.cpp
index 0df5b1b0..530b71cd 100644
--- a/cpu/BGK.cpp
+++ b/cpu/BGK.cpp
@@ -15,7 +15,6 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz){
-	int n;
 	// conserved momemnts
 	double rho,ux,uy,uz,uu;
 	// non-conserved moments
@@ -127,14 +126,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int
 }
 
 extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz){
-	int n;
 	// conserved momemnts
 	double rho,ux,uy,uz,uu;
 	// non-conserved moments
 	double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
 	int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
 
-	int nread;
 	for (int n=start; n<finish; n++){
 		
 		// q=0
@@ -291,4 +288,4 @@ extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int star
 				rlx*0.02777777777777778*(rho - 3.0*(uy-uz) + 4.5*(uy-uz)*(uy-uz) - uu) - 0.08333333333*(Fy-Fz);
 
 	}
-}
\ No newline at end of file
+}
diff --git a/cpu/Color.cpp b/cpu/Color.cpp
index 7396d97c..52c8dc3c 100644
--- a/cpu/Color.cpp
+++ b/cpu/Color.cpp
@@ -935,17 +935,14 @@ extern "C" void ScaLBL_D3Q19_ColorCollide( char *ID, double *disteven, double *d
 extern "C" void ScaLBL_D3Q7_ColorCollideMass(char *ID, double *A_even, double *A_odd, double *B_even, double *B_odd, 
 		double *Den, double *Phi, double *ColorGrad, double *Velocity, double beta, int N, bool pBC)
 {
+	int n;
 	char id;
-
-	int idx,n,q,Cqx,Cqy,Cqz;
-	//	int sendLoc;
-
 	double f0,f1,f2,f3,f4,f5,f6;
 	double na,nb,nab;		// density values
 	double ux,uy,uz;	// flow velocity
 	double nx,ny,nz,C;	// color gradient components
 	double a1,a2,b1,b2;
-	double sp,delta;
+	double delta;
 	//double feq[6];		// equilibrium distributions
 	// Set of Discrete velocities for the D3Q19 Model
 	//int D3Q7[3][3]={{1,0,0},{0,1,0},{0,0,1}};
@@ -1271,7 +1268,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, do
 		double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 		double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){
 
-	int ijk,nn,n;
+	int ijk,nn;
 	double fq;
 	// conserved momemnts
 	double rho,jx,jy,jz;
@@ -1854,7 +1851,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double *di
 		double *Phi, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 		double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){
 	
-	int n,nn,ijk,nread;
+	int nn,ijk,nread;
 	int nr1,nr2,nr3,nr4,nr5,nr6;
 	int nr7,nr8,nr9,nr10;
 	int nr11,nr12,nr13,nr14;
@@ -2514,7 +2511,6 @@ extern "C" void ScaLBL_D3Q7_AAodd_Color(int *neighborList, int *Map, double *Aq,
 	double a1,b1,a2,b2,nAB,delta;
 	double C,nx,ny,nz; //color gradient magnitude and direction
 	double ux,uy,uz;
-	double phi;
 	// Instantiate mass transport distributions
 	// Stationary value - distribution 0
 	for (int n=start; n<finish; n++){
@@ -2547,7 +2543,6 @@ extern "C" void ScaLBL_D3Q7_AAodd_Color(int *neighborList, int *Map, double *Aq,
 		nB = Den[Np + n];
 
 		// compute phase indicator field
-		phi=(nA-nB)/(nA+nB);
 		nAB = 1.0/(nA+nB);
 		Aq[n] = 0.3333333333333333*nA;
 		Bq[n] = 0.3333333333333333*nB;
@@ -2698,7 +2693,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Color(int *Map, double *Aq, double *Bq, doubl
 extern "C" void ScaLBL_D3Q7_AAodd_PhaseField(int *neighborList, int *Map, double *Aq, double *Bq, 
 			double *Den, double *Phi, int start, int finish, int Np){
 
-	int idx,nread;
+        int idx,nread;
 	double fq,nA,nB;
 
 	for (int n=start; n<finish; n++){
@@ -2858,7 +2853,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_PhaseField(int *Map, double *Aq, double *Bq,
 }
 
 extern "C" void ScaLBL_D3Q19_Gradient(int *Map, double *phi, double *ColorGrad, int start, int finish, int Np, int Nx, int Ny, int Nz){
-	int idx,n,N,i,j,k,nn;
+	int idx,n,i,j,k,nn;
 	// distributions
 	double f1,f2,f3,f4,f5,f6,f7,f8,f9;
 	double f10,f11,f12,f13,f14,f15,f16,f17,f18;
diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index f492d4da..135104df 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -161,6 +161,20 @@ void ScaLBL_ColorModel::ReadParams(string filename){
 	else if (domain_db->keyExists( "BC" )){
 		BoundaryCondition = domain_db->getScalar<int>( "BC" );
 	}
+	if (domain_db->keyExists( "InletLayersPhase" )){
+		int inlet_layers_phase = domain_db->getScalar<int>( "InletLayersPhase" );
+		if (inlet_layers_phase == 2 ) {
+		  inletA = 0.0;
+		  inletB = 1.0;
+		}
+	}
+	if (domain_db->keyExists( "OutletLayersPhase" )){
+		int outlet_layers_phase = domain_db->getScalar<int>( "OutletLayersPhase" );
+		if (outlet_layers_phase == 1 ) {
+		  inletA = 1.0;
+		  inletB = 0.0;
+		}
+	}
 	
 	// Override user-specified boundary condition for specific protocols
 	auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
@@ -192,6 +206,21 @@ void ScaLBL_ColorModel::ReadParams(string filename){
 		}
 		domain_db->putScalar<int>( "BC", BoundaryCondition );
 	}  	
+	else if (protocol == "centrifuge"){
+		if (BoundaryCondition != 3 ){
+			BoundaryCondition = 3;
+			if (rank==0) printf("WARNING: protocol (centrifuge) supports only constant pressure boundary condition \n");
+		}
+		domain_db->putScalar<int>( "BC", BoundaryCondition );
+	} 
+	else if (protocol == "core flooding"){
+		if (BoundaryCondition != 4){
+			BoundaryCondition = 4;
+			if (rank==0) printf("WARNING: protocol (core flooding) supports only volumetric flux boundary condition \n");
+		}
+		domain_db->putScalar<int>( "BC", BoundaryCondition );
+	} 
+	
 }
 
 void ScaLBL_ColorModel::SetDomain(){
@@ -322,11 +351,12 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
 	if (WettingConvention == "SCAL"){
 	  for (size_t idx=0; idx<NLABELS; idx++) AffinityList[idx] *= -1.0;
 	}
-
-	double label_count[NLABELS];
-	double label_count_global[NLABELS];
+	
+	double * label_count;
+	double * label_count_global;
+	label_count = new double [NLABELS];
+	label_count_global = new double [NLABELS];
 	// Assign the labels
-
 	for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
 
 	for (int k=0;k<Nz;k++){
@@ -585,7 +615,7 @@ double ScaLBL_ColorModel::Run(int returntime){
 	if (analysis_db->keyExists( "tolerance" )){
 		tolerance = analysis_db->getScalar<double>( "tolerance" );
 	}
-
+	
 	runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map );
 	auto t1 = std::chrono::system_clock::now();
 	int CURRENT_TIMESTEP = 0;
@@ -894,32 +924,7 @@ void ScaLBL_ColorModel::Run(){
 	if (color_db->keyExists( "krA_morph_factor" )){
 		KRA_MORPH_FACTOR  = color_db->getScalar<double>( "krA_morph_factor" );
 	}
-	/* defaults for simulation protocols */
-	auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
-	if (protocol == "image sequence"){
-		// Get the list of images
-		USE_DIRECT = true;
-		ImageList = color_db->getVector<std::string>( "image_sequence");
-		IMAGE_INDEX = color_db->getWithDefault<int>( "image_index", 0 );
-		IMAGE_COUNT = ImageList.size();
-		morph_interval = 10000;
-		USE_MORPH = true;
-	}
-	else if (protocol == "seed water"){
-		morph_delta = -0.05;
-		seed_water = 0.01;
-		USE_SEED = true;
-		USE_MORPH = true;
-	}
-	else if (protocol == "open connected oil"){
-		morph_delta = -0.05;
-		USE_MORPH = true;
-		USE_MORPHOPEN_OIL = true;
-	}
-	else if (protocol == "shell aggregation"){
-		morph_delta = -0.05;
-		USE_MORPH = true;
-	}  
+
 	if (color_db->keyExists( "capillary_number" )){
 		capillary_number = color_db->getScalar<double>( "capillary_number" );
 		SET_CAPILLARY_NUMBER=true;
@@ -938,7 +943,6 @@ void ScaLBL_ColorModel::Run(){
 	if (analysis_db->keyExists( "seed_water" )){
 		seed_water = analysis_db->getScalar<double>( "seed_water" );
 		if (rank == 0) printf("Seed water in oil %f (seed_water) \n",seed_water);
-		ASSERT(protocol == "seed water");
 	}
 	if (analysis_db->keyExists( "morph_delta" )){
 		morph_delta = analysis_db->getScalar<double>( "morph_delta" );
@@ -968,7 +972,54 @@ void ScaLBL_ColorModel::Run(){
 	if (analysis_db->keyExists( "max_morph_timesteps" )){
 		MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" );
 	}
-
+	
+	/* defaults for simulation protocols */
+	auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
+	if (protocol == "image sequence"){
+		// Get the list of images
+		USE_DIRECT = true;
+		ImageList = color_db->getVector<std::string>( "image_sequence");
+		IMAGE_INDEX = color_db->getWithDefault<int>( "image_index", 0 );
+		IMAGE_COUNT = ImageList.size();
+		morph_interval = 10000;
+		USE_MORPH = true;
+		USE_SEED = false;
+	}
+	else if (protocol == "seed water"){
+		morph_delta = -0.05;
+		seed_water = 0.01;
+		USE_SEED = true;
+		USE_MORPH = true;
+	}
+	else if (protocol == "open connected oil"){
+		morph_delta = -0.05;
+		USE_SEED = false;
+		USE_MORPH = true;
+		USE_MORPHOPEN_OIL = true;
+	}
+	else if (protocol == "shell aggregation"){
+		morph_delta = -0.05;
+		USE_MORPH = true;
+		USE_SEED = false;
+	}  
+	else if (protocol == "fractional flow"){
+		USE_MORPH = false;
+		USE_SEED = false;
+	}
+	else if (protocol == "centrifuge"){
+		USE_MORPH = false;
+		USE_SEED = false;
+	} 
+	else if (protocol == "core flooding"){
+		USE_MORPH = false;
+		USE_SEED = false;
+		if (SET_CAPILLARY_NUMBER){
+			double MuB = rhoB*(tauB - 0.5)/3.0;
+			double IFT = 6.0*alpha;
+			double CrossSectionalArea = (double) (nprocx*(Nx-2)*nprocy*(Ny-2));
+			flux = Dm->Porosity()*CrossSectionalArea*IFT*capillary_number/MuB;
+		}
+	} 
 	if (rank==0){
 		printf("********************************************************\n");
 		if (protocol == "image sequence"){
@@ -1006,7 +1057,20 @@ void ScaLBL_ColorModel::Run(){
 			printf("     min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS);
 			printf("     max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS);
 			printf("     tolerance = %f \n",tolerance);
-			printf("     morph_delta = %f \n",morph_delta);
+		} 
+		else if (protocol == "centrifuge"){
+			printf("  using protocol =  centrifuge \n"); 
+			printf("     driving force = %f \n",Fz);
+			if (Fz < 0){
+				printf("      Component B displacing component A \n");				
+			}
+			else if (Fz > 0){
+				printf("      Component A displacing component B \n");				
+			}
+		} 
+		else if (protocol == "core flooding"){
+			printf("  using protocol = core flooding \n"); 
+			printf("     capillary number = %f \n", capillary_number);
 		} 
 		printf("No. of timesteps: %i \n", timestepMax);
 		fflush(stdout);
@@ -1120,7 +1184,7 @@ void ScaLBL_ColorModel::Run(){
 			double volB = Averages->gwb.V; 
 			double volA = Averages->gnb.V; 
 			volA /= Dm->Volume;
-			volB /= Dm->Volume;;
+			volB /= Dm->Volume;
 			//initial_volume = volA*Dm->Volume;
 			double vA_x = Averages->gnb.Px/Averages->gnb.M; 
 			double vA_y = Averages->gnb.Py/Averages->gnb.M; 
@@ -1944,12 +2008,11 @@ FlowAdaptor::~FlowAdaptor(){
 
 double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){	
 	
-  	  double MASS_FRACTION_CHANGE = 0.05;
-	  if (M.db->keyExists( "FlowAdaptor" )){
-		  auto flow_db = M.db->getDatabase( "FlowAdaptor" );
-		  MASS_FRACTION_CHANGE = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
-	  }
-
+  	  double MASS_FRACTION_CHANGE = 0.0005;
+  	  if (M.db->keyExists( "FlowAdaptor" )){
+  		  auto flow_db = M.db->getDatabase( "FlowAdaptor" );
+  		  MASS_FRACTION_CHANGE = flow_db->getWithDefault<double>( "mass_fraction_factor", 0.0005);
+  	  }
 	  int Np = M.Np;
 	  double dA, dB, phi;
 	  double vx,vy,vz;
@@ -1975,12 +2038,53 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 	  ScaLBL_CopyToHost(Vel_x, &M.Velocity[0], Np*sizeof(double));
 	  ScaLBL_CopyToHost(Vel_y, &M.Velocity[Np], Np*sizeof(double));
 	  ScaLBL_CopyToHost(Vel_z, &M.Velocity[2*Np], Np*sizeof(double));
-
+	  
+	  /* DEBUG STRUCTURES */
+	  int Nx = M.Nx;  int Ny = M.Ny;  int Nz = M.Nz;
+	  int N = Nx*Ny*Nz;
+	  double * DebugMassA, *DebugMassB;
+	  DebugMassA = new double[Np];
+	  DebugMassB = new double[Np];
+	  
 	  /* compute the total momentum */
 	  vax = vay = vaz = 0.0;
 	  vbx = vby = vbz = 0.0;
 	  mass_a = mass_b = 0.0;
-	  for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){
+	  
+	  double maxSpeed = 0.0;
+	  double localMaxSpeed = 0.0;
+	  for (int k=1; k<Nz-1; k++){
+		  for (int j=1; j<Ny-1; j++){
+			  for (int i=1; i<Nx-1; i++){
+				  int n=M.Map(i,j,k);
+				  double distance = M.Averages->SDs(i,j,k);
+				  if (!(n<0) ){
+					  dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
+					  dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
+					  phi = (dA - dB) / (dA + dB);
+					  Phase[n] = phi;
+					  mass_a += dA;
+					  mass_b += dB;
+					  if (phi > 0.0){
+						  vax += Vel_x[n];
+						  vay += Vel_y[n];
+						  vaz += Vel_z[n];
+					  }
+					  else {
+						  vbx += Vel_x[n];
+						  vby += Vel_y[n];
+						  vbz += Vel_z[n];
+					  }
+					  double speed =  sqrt(vax*vax + vay*vay + vaz*vaz + vbx*vbx + vby*vby + vbz*vbz);
+					  if (distance > 3.0 && speed > localMaxSpeed){
+						  localMaxSpeed = speed;
+					  }
+				  }
+			  }
+		  }
+	  }
+	  maxSpeed = M.Dm->Comm.maxReduce(localMaxSpeed);
+	  /*	  for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){
 		  dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
 		  dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
 		  phi = (dA - dB) / (dA + dB);
@@ -2016,27 +2120,65 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 				  vbz += Vel_z[n];
 			  }
 	   }
-	   mass_a_global = M.Dm->Comm.sumReduce(mass_a);
-	   mass_b_global = M.Dm->Comm.sumReduce(mass_b);
-	   vax_global = M.Dm->Comm.sumReduce(vax);
-	   vay_global = M.Dm->Comm.sumReduce(vay);
-	   vaz_global = M.Dm->Comm.sumReduce(vaz);
-	   vbx_global = M.Dm->Comm.sumReduce(vbx);
-	   vby_global = M.Dm->Comm.sumReduce(vby);
-	   vbz_global = M.Dm->Comm.sumReduce(vbz);	 
-	   
-	   double  total_momentum_A = sqrt(vax_global*vax_global+vay_global*vay_global+vaz_global*vaz_global);
-	   double  total_momentum_B = sqrt(vbx_global*vbx_global+vby_global*vby_global+vbz_global*vbz_global);
-	   
-	   /* compute the total mass change */
-	   double TOTAL_MASS_CHANGE = MASS_FRACTION_CHANGE*(mass_a_global + mass_b_global);
-	   if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_a_global )
-		   TOTAL_MASS_CHANGE = 0.1*mass_a_global;
-	   if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_b_global )
-		   TOTAL_MASS_CHANGE = 0.1*mass_b_global;
-	   
-	   double LOCAL_MASS_CHANGE = 0.0;
-	   for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){	    
+	   */
+	  mass_a_global = M.Dm->Comm.sumReduce(mass_a);
+	  mass_b_global = M.Dm->Comm.sumReduce(mass_b);
+	  vax_global = M.Dm->Comm.sumReduce(vax);
+	  vay_global = M.Dm->Comm.sumReduce(vay);
+	  vaz_global = M.Dm->Comm.sumReduce(vaz);
+	  vbx_global = M.Dm->Comm.sumReduce(vbx);
+	  vby_global = M.Dm->Comm.sumReduce(vby);
+	  vbz_global = M.Dm->Comm.sumReduce(vbz);	 
+
+	  double  total_momentum_A = sqrt(vax_global*vax_global+vay_global*vay_global+vaz_global*vaz_global);
+	  double  total_momentum_B = sqrt(vbx_global*vbx_global+vby_global*vby_global+vbz_global*vbz_global);
+	  /* compute the total mass change */
+	  double TOTAL_MASS_CHANGE = MASS_FRACTION_CHANGE*(mass_a_global + mass_b_global);
+	  if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_a_global )
+		  TOTAL_MASS_CHANGE = 0.1*mass_a_global;
+	  if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_b_global )
+		  TOTAL_MASS_CHANGE = 0.1*mass_b_global;
+
+	  double LOCAL_MASS_CHANGE = 0.0;
+	  for (int k=1; k<Nz-1; k++){
+		  for (int j=1; j<Ny-1; j++){
+			  for (int i=1; i<Nx-1; i++){
+				  int n=M.Map(i,j,k);
+				  if (!(n<0)){
+					  phi = Phase[n];
+					  vx = Vel_x[n];
+					  vy = Vel_y[n];
+					  vz = Vel_z[n];
+					  double local_momentum = sqrt(vx*vx+vy*vy+vz*vz);
+					  /* impose ceiling for spurious currents */
+					  if (local_momentum > maxSpeed) local_momentum =  maxSpeed;
+					  if (phi > 0.0){
+						  LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_A;
+						  Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE;
+						  Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  Aq_tmp[n+3*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE;
+					  }
+					  else{
+						  LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B;
+						  Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE;
+						  Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  Bq_tmp[n+3*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+						  DebugMassB[n] = LOCAL_MASS_CHANGE;
+					  }
+				  }
+			  }
+		  }
+	  }
+/*	   for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){	    
 		   phi = Phase[n];
 		   vx = Vel_x[n];
 		   vy = Vel_y[n];
@@ -2051,6 +2193,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 			   Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
 			   Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
 			   Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+			   DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE;
 		   }
 		   else{
 			   LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B;
@@ -2061,6 +2204,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 			   Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
 			   Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
 			   Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+			   DebugMassB[n] = LOCAL_MASS_CHANGE;
 		   }
 	   }
 
@@ -2079,6 +2223,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 			   Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
 			   Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
 			   Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
+			   DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE;
 		   }
 		   else{
 			   LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B;
@@ -2089,11 +2234,54 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 			   Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
 			   Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
 			   Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
+			   DebugMassB[n] = LOCAL_MASS_CHANGE;
 		   }
 	  }
-
+ 	 */
 	  if (M.rank == 0) printf("Update Fractional Flow: change mass of fluid B by %f \n",TOTAL_MASS_CHANGE/mass_b_global);
 
+	  /* Print out debugging info with mass update  */
+	  // initialize the array
+	  double value;
+	  char LocalRankFilename[40];
+	  DoubleArray regdata(Nx,Ny,Nz);
+	  regdata.fill(0.f);
+	  for (int k=0; k<Nz; k++){
+		  for (int j=0; j<Ny; j++){
+			  for (int i=0; i<Nx; i++){
+				  int idx=M.Map(i,j,k);
+				  if (!(idx<0)){
+					  value=DebugMassA[idx];
+					  regdata(i,j,k)=value;
+				  }
+			  }
+		  }
+	  }
+	  
+	  FILE *AFILE;
+	  sprintf(LocalRankFilename,"dA.%05i.raw",M.rank);
+	  AFILE = fopen(LocalRankFilename,"wb");
+	  fwrite(regdata.data(),8,N,AFILE);
+	  fclose(AFILE);
+
+	  regdata.fill(0.f);
+	  for (int k=0; k<Nz; k++){
+		  for (int j=0; j<Ny; j++){
+			  for (int i=0; i<Nx; i++){
+				  int idx=M.Map(i,j,k);
+				  if (!(idx<0)){
+					  value=DebugMassB[idx];
+					  regdata(i,j,k)=value;
+				  }
+			  }
+		  }
+	  }
+	  FILE *BFILE;
+	  sprintf(LocalRankFilename,"dB.%05i.raw",M.rank);
+	  BFILE = fopen(LocalRankFilename,"wb");
+	  fwrite(regdata.data(),8,N,BFILE);
+	  fclose(BFILE);
+
 	  // Need to initialize Aq, Bq, Den, Phi directly
 	  //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
 	  ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
@@ -2105,12 +2293,9 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
 void FlowAdaptor::Flatten(ScaLBL_ColorModel &M){	
 	
 	  int Np = M.Np;
-	  double dA, dB, phi;
-
-	  double mass_a, mass_b, mass_a_global, mass_b_global;
+	  double dA, dB;
 	  
 	  double *Aq_tmp, *Bq_tmp;
-	  double *Vel_x, *Vel_y, *Vel_z, *Phase;
 	  
 	  Aq_tmp = new double [7*Np];
 	  Bq_tmp = new double [7*Np];
@@ -2118,7 +2303,6 @@ void FlowAdaptor::Flatten(ScaLBL_ColorModel &M){
 	  ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double));
 	  ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
 
-
 	  for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){
 		  dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
 		  dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
@@ -2202,22 +2386,6 @@ double FlowAdaptor::MoveInterface(ScaLBL_ColorModel &M){
     	}
 	}
     ScaLBL_CopyToDevice( M.Phi, phi_t.data(), Nx*Ny*Nz* sizeof( double ) );	
-    
-    
-/*	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
-	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
-	if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
-		if (Dm->kproc()==0){
-			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
-			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
-			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
-		}
-		if (Dm->kproc() == nprocz-1){
-			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
-			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
-			ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
-		}
-	}
-	*/
+    return total_interface_sites;
 }
 
diff --git a/models/FreeLeeModel.cpp b/models/FreeLeeModel.cpp
index b1c0e99f..19a2afa8 100644
--- a/models/FreeLeeModel.cpp
+++ b/models/FreeLeeModel.cpp
@@ -63,7 +63,6 @@ void ScaLBL_FreeLeeModel::getData_RegularLayout(const double *data, DoubleArray
 	// Gets data (in optimized layout) from the HOST and stores in regular layout
     // Primarly for debugging
 	int i,j,k,idx;
-    int n;
 
 	// initialize the array
 	regdata.fill(0.f);
@@ -72,7 +71,6 @@ void ScaLBL_FreeLeeModel::getData_RegularLayout(const double *data, DoubleArray
 	for (k=0; k<Nz; k++){
 		for (j=0; j<Ny; j++){
 			for (i=0; i<Nx; i++){
-				n=k*Nx*Ny+j*Nx+i;
 				idx=Map(i,j,k);
 				if (!(idx<0)){
 					value=data[idx];
@@ -414,8 +412,10 @@ void ScaLBL_FreeLeeModel::AssignComponentLabels_ChemPotential_ColorGrad()
 		ERROR("Error: ComponentLabels and ComponentAffinity must be the same length! \n");
 	}
 
-	double label_count[NLABELS];
-	double label_count_global[NLABELS];
+	double *label_count;
+	double *label_count_global;
+	label_count = new double [NLABELS];
+	label_count_global = new double [NLABELS];
 
 	// Assign the labels
 	for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
@@ -738,75 +738,10 @@ void ScaLBL_FreeLeeModel::Initialize_SingleFluid(){
 	if (Restart == true){
         //TODO need to revise this function
         //remove the phase-related part
-
-
-
-//		if (rank==0){
-//			printf("Reading restart file! \n");
-//		}
-//
-//		// Read in the restart file to CPU buffers
-//		int *TmpMap;
-//		TmpMap = new int[Np];
-//		
-//		double *cPhi, *cDist, *cDen;
-//		cPhi = new double[N];
-//		cDen = new double[2*Np];
-//		cDist = new double[19*Np];
-//		ScaLBL_CopyToHost(TmpMap, dvcMap, Np*sizeof(int));
-//        //ScaLBL_CopyToHost(cPhi, Phi, N*sizeof(double));
-//    	
-//		ifstream File(LocalRestartFile,ios::binary);
-//		int idx;
-//		double value,va,vb;
-//		for (int n=0; n<Np; n++){
-//			File.read((char*) &va, sizeof(va));
-//			File.read((char*) &vb, sizeof(vb));
-//			cDen[n]    = va;
-//			cDen[Np+n] = vb;
-//		}
-//		for (int n=0; n<Np; n++){
-//			// Read the distributions
-//			for (int q=0; q<19; q++){
-//				File.read((char*) &value, sizeof(value));
-//				cDist[q*Np+n] = value;
-//			}
-//		}
-//		File.close();
-//		
-//		for (int n=0; n<ScaLBL_Comm->LastExterior(); n++){
-//			va = cDen[n];
-//			vb = cDen[Np + n];
-//			value = (va-vb)/(va+vb);
-//			idx = TmpMap[n];
-//			if (!(idx < 0) && idx<N)
-//				cPhi[idx] = value;
-//		}
-//		for (int n=ScaLBL_Comm->FirstInterior(); n<ScaLBL_Comm->LastInterior(); n++){
-//		  va = cDen[n];
-//		  vb = cDen[Np + n];
-//		  	value = (va-vb)/(va+vb);
-//		  	idx = TmpMap[n];
-//		  	if (!(idx < 0) && idx<N)
-//		  		cPhi[idx] = value;
-//		}
-//		
-//		// Copy the restart data to the GPU
-//		ScaLBL_CopyToDevice(Den,cDen,2*Np*sizeof(double));
-//		ScaLBL_CopyToDevice(gqbar,cDist,19*Np*sizeof(double));
-//		ScaLBL_CopyToDevice(Phi,cPhi,N*sizeof(double));
-//		ScaLBL_Comm->Barrier();
-//		comm.barrier();
-//
-//        if (rank==0)	printf ("Initializing phase and density fields on device from Restart\n");
-//        //TODO the following function is to be updated.
-//        //ScaLBL_FreeLeeModel_PhaseField_InitFromRestart(Den, hq, 0, ScaLBL_Comm->LastExterior(), Np);
-//        //ScaLBL_FreeLeeModel_PhaseField_InitFromRestart(Den, hq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 	}
 }
 
 double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
-	int nprocs=nprocx*nprocy*nprocz;
 	
 	int START_TIME = timestep;
 	int EXIT_TIME = min(returntime, timestepMax);
diff --git a/models/GreyscaleModel.cpp b/models/GreyscaleModel.cpp
index cb560ec2..d15cdf7f 100644
--- a/models/GreyscaleModel.cpp
+++ b/models/GreyscaleModel.cpp
@@ -165,7 +165,6 @@ void ScaLBL_GreyscaleModel::ReadInput(){
 	// Generate the signed distance map
 	// Initialize the domain and communication
 	Array<char> id_solid(Nx,Ny,Nz);
-	int count = 0;
 	// Solve for the position of the solid phase
 	for (int k=0;k<Nz;k++){
 		for (int j=0;j<Ny;j++){
@@ -182,7 +181,6 @@ void ScaLBL_GreyscaleModel::ReadInput(){
 	for (int k=0;k<Nz;k++){
 		for (int j=0;j<Ny;j++){
 			for (int i=0;i<Nx;i++){
-				int n=k*Nx*Ny+j*Nx+i;
 				// Initialize distance to +/- 1
 				SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
 			}
@@ -214,11 +212,13 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
 		ERROR("Error: ComponentLabels and PorosityList must be the same length! \n");
 	}
 
-	double label_count[NLABELS];
-	double label_count_global[NLABELS];
 	// Assign the labels
-
-	for (int idx=0; idx<NLABELS; idx++) label_count[idx]=0;
+	double *label_count;
+	double *label_count_global;
+	label_count = new double [NLABELS];
+	label_count_global = new double [NLABELS];
+	
+	for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
 
 	for (int k=0;k<Nz;k++){
 		for (int j=0;j<Ny;j++){
@@ -226,7 +226,7 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
 				int n = k*Nx*Ny+j*Nx+i;
 				VALUE=id[n];
 				// Assign the affinity from the paired list
-				for (unsigned int idx=0; idx < NLABELS; idx++){
+				for (size_t idx=0; idx < NLABELS; idx++){
 				      //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
 					if (VALUE == LabelList[idx]){
 						POROSITY=PorosityList[idx];
@@ -257,7 +257,7 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
 				int n = k*Nx*Ny+j*Nx+i;
 				VALUE=id[n];
 				// Assign the affinity from the paired list
-				for (unsigned int idx=0; idx < NLABELS; idx++){
+				for (size_t idx=0; idx < NLABELS; idx++){
 					//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
 					if (VALUE == LabelList[idx]){
 						PERMEABILITY=PermeabilityList[idx];
@@ -282,7 +282,7 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
 	// Set Dm to match Mask
 	for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i]; 
 	
-	for (int idx=0; idx<NLABELS; idx++)		label_count_global[idx]=Dm->Comm.sumReduce(  label_count[idx]);
+	for (size_t idx=0; idx<NLABELS; idx++)		label_count_global[idx]=Dm->Comm.sumReduce(  label_count[idx]);
     //Initialize a weighted porosity after considering grey voxels
     GreyPorosity=0.0;
 	for (unsigned int idx=0; idx<NLABELS; idx++){
@@ -613,7 +613,7 @@ void ScaLBL_GreyscaleModel::Run(){
             //double mass_loc,mass_glb;
             
             //parameters for domain average
-	        int64_t i,j,k,n,imin,jmin,kmin,kmax;
+	        int64_t imin,jmin,kmin,kmax;
             // If external boundary conditions are set, do not average over the inlet and outlet
             kmin=1; kmax=Nz-1;
             //In case user forgets to specify the inlet/outlet buffer layers for BC>0
@@ -706,23 +706,22 @@ void ScaLBL_GreyscaleModel::Run(){
 			//double absperm = h*h*mu*Mask->Porosity()*flow_rate / force_mag;
 			double absperm = h*h*mu*GreyPorosity*flow_rate / force_mag;
 
-            if (rank==0){
+			if (rank==0){
 				printf("     AbsPerm = %.5g [micron^2]\n",absperm);
-                bool WriteHeader=false;
-                FILE * log_file = fopen("Permeability.csv","r");
-                if (log_file != NULL)
-                    fclose(log_file);
-                else
-                    WriteHeader=true;
-                log_file = fopen("Permeability.csv","a");
-                if (WriteHeader)
-                    fprintf(log_file,"timestep Fx Fy Fz mu Vs As Hs Xs vax vay vaz AbsPerm \n",
-                            timestep,Fx,Fy,Fz,mu,h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz,absperm);
+				bool WriteHeader=false;
+				FILE * log_file = fopen("Permeability.csv","r");
+				if (log_file != NULL)
+					fclose(log_file);
+				else
+					WriteHeader=true;
+				log_file = fopen("Permeability.csv","a");
+				if (WriteHeader)
+					fprintf(log_file,"timestep Fx Fy Fz mu Vs As Hs Xs vax vay vaz AbsPerm \n");
 
-                fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, 
-                        h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm);
-                fclose(log_file);
-            }
+				fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, 
+						h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm);
+				fclose(log_file);
+			}
 		}
 
 		if (timestep%visualization_interval==0){
diff --git a/models/IonModel.cpp b/models/IonModel.cpp
index 136af02f..f86bc6d2 100644
--- a/models/IonModel.cpp
+++ b/models/IonModel.cpp
@@ -97,7 +97,7 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
         }
         else{
             time_conv.clear();
-            for (int i=0; i<tau.size();i++){
+            for (size_t i=0; i<tau.size();i++){
                 time_conv.push_back((tau[i]-0.5)/k2_inv*(h*h*1.0e-12)/Di[i]);
             }
         }
@@ -112,13 +112,13 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
 		    ERROR("Error: number_ion_species and IonDiffusivityList must be the same length! \n");
         }
         else{
-            for (int i=0; i<IonDiffusivity.size();i++){
+            for (size_t i=0; i<IonDiffusivity.size();i++){
                 IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
             }
         }
     }
     else {
-        for (int i=0; i<IonDiffusivity.size();i++){
+        for (size_t i=0; i<IonDiffusivity.size();i++){
             //convert ion diffusivity in physical unit to LB unit
             IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
         }
@@ -141,13 +141,13 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
 		    ERROR("Error: number_ion_species and IonConcentrationList must be the same length! \n");
         }
         else{
-            for (int i=0; i<IonConcentration.size();i++){
+            for (size_t i=0; i<IonConcentration.size();i++){
                 IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
             }
         }
     }
     else {
-        for (int i=0; i<IonConcentration.size();i++){
+        for (size_t i=0; i<IonConcentration.size();i++){
             IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
         }
     
@@ -186,7 +186,7 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
             else {
                 ERROR("Error: Non-periodic BCs are specified but InletValueList cannot be found! \n");
             }
-            for (unsigned int i=0;i<BoundaryConditionInlet.size();i++){
+            for (size_t i=0;i<BoundaryConditionInlet.size();i++){
                 switch (BoundaryConditionInlet[i]){
                     case 1://fixed boundary ion concentration [mol/m^3]
                        Cin[i] = Cin[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
@@ -220,7 +220,7 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
             else {
                 ERROR("Error: Non-periodic BCs are specified but OutletValueList cannot be found! \n");
             }
-            for (unsigned int i=0;i<BoundaryConditionOutlet.size();i++){
+            for (size_t i=0;i<BoundaryConditionOutlet.size();i++){
                 switch (BoundaryConditionOutlet[i]){
                     case 1://fixed boundary ion concentration [mol/m^3]
                        Cout[i] = Cout[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
@@ -307,7 +307,7 @@ void ScaLBL_IonModel::ReadParams(string filename){
         }
         else{
             time_conv.clear();
-            for (int i=0; i<tau.size();i++){
+            for (size_t i=0; i<tau.size();i++){
                 time_conv.push_back((tau[i]-0.5)/k2_inv*(h*h*1.0e-12)/Di[i]);
             }
         }
@@ -322,13 +322,13 @@ void ScaLBL_IonModel::ReadParams(string filename){
 		    ERROR("Error: number_ion_species and IonDiffusivityList must be the same length! \n");
         }
         else{
-            for (int i=0; i<IonDiffusivity.size();i++){
+            for (size_t i=0; i<IonDiffusivity.size();i++){
                 IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
             }
         }
     }
     else {
-        for (int i=0; i<IonDiffusivity.size();i++){
+        for (size_t i=0; i<IonDiffusivity.size();i++){
             //convert ion diffusivity in physical unit to LB unit
             IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
         }
@@ -351,13 +351,13 @@ void ScaLBL_IonModel::ReadParams(string filename){
 		    ERROR("Error: number_ion_species and IonConcentrationList must be the same length! \n");
         }
         else{
-            for (int i=0; i<IonConcentration.size();i++){
+            for (size_t i=0; i<IonConcentration.size();i++){
                 IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
             }
         }
     }
     else {
-        for (int i=0; i<IonConcentration.size();i++){
+        for (size_t i=0; i<IonConcentration.size();i++){
             IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
         }
     
@@ -396,7 +396,7 @@ void ScaLBL_IonModel::ReadParams(string filename){
             else {
                 ERROR("Error: Non-periodic BCs are specified but InletValueList cannot be found! \n");
             }
-            for (unsigned int i=0;i<BoundaryConditionInlet.size();i++){
+            for (size_t i=0;i<BoundaryConditionInlet.size();i++){
                 switch (BoundaryConditionInlet[i]){
                     case 1://fixed boundary ion concentration [mol/m^3]
                        Cin[i] = Cin[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
@@ -430,7 +430,7 @@ void ScaLBL_IonModel::ReadParams(string filename){
             else {
                 ERROR("Error: Non-periodic BCs are specified but OutletValueList cannot be found! \n");
             }
-            for (unsigned int i=0;i<BoundaryConditionOutlet.size();i++){
+            for (size_t i=0;i<BoundaryConditionOutlet.size();i++){
                 switch (BoundaryConditionOutlet[i]){
                     case 1://fixed boundary ion concentration [mol/m^3]
                        Cout[i] = Cout[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
@@ -558,10 +558,12 @@ void ScaLBL_IonModel::AssignSolidBoundary(double *ion_solid)
 		ERROR("Error: LB Ion Solver: SolidLabels and SolidValues must be the same length! \n");
 	}
 
-	double label_count[NLABELS];
-	double label_count_global[NLABELS];
-	// Assign the labels
 
+	// Assign the labels
+	double *label_count;
+	double *label_count_global;
+	label_count = new double [NLABELS];
+	label_count_global = new double [NLABELS];
 	for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
 
 	for (int k=0;k<Nz;k++){
@@ -571,7 +573,7 @@ void ScaLBL_IonModel::AssignSolidBoundary(double *ion_solid)
 				VALUE=Mask->id[n];
                 AFFINITY=0.f;
 				// Assign the affinity from the paired list
-				for (unsigned int idx=0; idx < NLABELS; idx++){
+				for (size_t idx=0; idx < NLABELS; idx++){
 				      //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
 					if (VALUE == LabelList[idx]){
 						AFFINITY=AffinityList[idx];
@@ -701,23 +703,23 @@ void ScaLBL_IonModel::Initialize(){
 		auto File_ion = ion_db->getVector<std::string>( "IonConcentrationFile" );
         double *Ci_host;
         Ci_host = new double[number_ion_species*Np];
-        for (int ic=0; ic<number_ion_species; ic++){
+        for (size_t ic=0; ic<number_ion_species; ic++){
             AssignIonConcentration_FromFile(&Ci_host[ic*Np],File_ion);
         }
 	    ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species*sizeof(double)*Np);
 	    comm.barrier();
-        for (int ic=0; ic<number_ion_species; ic++){
+        for (size_t ic=0; ic<number_ion_species; ic++){
             ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic*Np*7],&Ci[ic*Np],Np); 
         }
         delete [] Ci_host;
     }
     else{
-        for (int ic=0; ic<number_ion_species; ic++){
+        for (size_t ic=0; ic<number_ion_species; ic++){
             ScaLBL_D3Q7_Ion_Init(&fq[ic*Np*7],&Ci[ic*Np],IonConcentration[ic],Np); 
         }
     }
     if (rank==0)    printf ("LB Ion Solver: initializing charge density\n");
-	for (int ic=0; ic<number_ion_species; ic++){
+	for (size_t ic=0; ic<number_ion_species; ic++){
         ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
         ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np);
     }
@@ -734,35 +736,35 @@ void ScaLBL_IonModel::Initialize(){
           break;
     }
 
-    for (int i=0; i<number_ion_species;i++){
+    for (size_t i=0; i<number_ion_species;i++){
         switch (BoundaryConditionInlet[i]){
             case 0:
-                if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %i is periodic \n",i+1);
+                if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %zu is periodic \n",i+1);
                 break;
             case 1:
-                if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %i is concentration = %.5g [mol/m^3] \n",i+1,Cin[i]/(h*h*h*1.0e-18));
+                if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %zu is concentration = %.5g [mol/m^3] \n",i+1,Cin[i]/(h*h*h*1.0e-18));
                 break;
             case 2:
-                if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %i is (inward) flux = %.5g [mol/m^2/sec] \n",i+1,Cin[i]/(h*h*1.0e-12)/time_conv[i]);
+                if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %zu is (inward) flux = %.5g [mol/m^2/sec] \n",i+1,Cin[i]/(h*h*1.0e-12)/time_conv[i]);
                 break;
         }
         switch (BoundaryConditionOutlet[i]){
             case 0:
-                if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %i is periodic \n",i+1);
+                if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %zu is periodic \n",i+1);
                 break;
             case 1:
-                if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %i is concentration = %.5g [mol/m^3] \n",i+1,Cout[i]/(h*h*h*1.0e-18));
+                if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %zu is concentration = %.5g [mol/m^3] \n",i+1,Cout[i]/(h*h*h*1.0e-18));
                 break;
             case 2:
-                if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %i is (inward) flux = %.5g [mol/m^2/sec] \n",i+1,Cout[i]/(h*h*1.0e-12)/time_conv[i]);
+                if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %zu is (inward) flux = %.5g [mol/m^2/sec] \n",i+1,Cout[i]/(h*h*1.0e-12)/time_conv[i]);
                 break;
         }
     }
 
 	if (rank==0) printf("*****************************************************\n");
 	if (rank==0) printf("LB Ion Transport Solver: \n");
-    for (int i=0; i<number_ion_species;i++){
-	    if (rank==0) printf("      Ion %i: LB relaxation tau = %.5g\n", i+1,tau[i]);
+    for (size_t i=0; i<number_ion_species;i++){
+	    if (rank==0) printf("      Ion %zu: LB relaxation tau = %.5g\n", i+1,tau[i]);
 	    if (rank==0) printf("              Time conversion factor: %.5g [sec/lt]\n", time_conv[i]);
 	    if (rank==0) printf("              Internal iteration: %i [lt]\n", timestepMax[i]);
     }
@@ -777,7 +779,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
 
     //LB-related parameter
     vector<double> rlx;
-    for (unsigned int ic=0;ic<tau.size();ic++){
+    for (size_t ic=0;ic<tau.size();ic++){
         rlx.push_back(1.0/tau[ic]); 
     }
     
@@ -786,7 +788,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
 	//ScaLBL_Comm->Barrier(); comm.barrier();
     //auto t1 = std::chrono::system_clock::now();
 
-	for (int ic=0; ic<number_ion_species; ic++){
+	for (size_t ic=0; ic<number_ion_species; ic++){
         timestep=0;
         while (timestep < timestepMax[ic]) {
             //************************************************************************/
@@ -881,7 +883,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
     }
 
     //Compute charge density for Poisson equation
-	for (int ic=0; ic<number_ion_species; ic++){
+	for (size_t ic=0; ic<number_ion_species; ic++){
         ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
         ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np);
     }
@@ -901,7 +903,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
 	//if (rank==0) printf("********************************************************\n");
 }
 
-void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const int ic){
+void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const size_t ic){
 	//This function wirte out the data in a normal layout (by aggregating all decomposed domains)
 
 	ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],IonConcentration);
@@ -913,13 +915,13 @@ void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const i
 void ScaLBL_IonModel::getIonConcentration_debug(int timestep){
     //This function write out decomposed data
     DoubleArray PhaseField(Nx,Ny,Nz);
-	for (int ic=0; ic<number_ion_species; ic++){
+	for (size_t ic=0; ic<number_ion_species; ic++){
 	    ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],PhaseField);
         ScaLBL_Comm->Barrier(); comm.barrier();
         IonConcentration_LB_to_Phys(PhaseField);
 
         FILE *OUTFILE;
-        sprintf(LocalRankFilename,"Ion%02i_Time_%i.%05i.raw",ic+1,timestep,rank);
+        sprintf(LocalRankFilename,"Ion%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
         OUTFILE = fopen(LocalRankFilename,"wb");
         fwrite(PhaseField.data(),8,N,OUTFILE);
         fclose(OUTFILE);
@@ -986,7 +988,7 @@ double ScaLBL_IonModel::CalIonDenConvergence(vector<double> &ci_avg_previous){
     Ci_host = new double[Np];
     vector<double> error(number_ion_species,0.0);
 
-	for (int ic=0; ic<number_ion_species; ic++){
+	for (size_t ic=0; ic<number_ion_species; ic++){
 
 	    ScaLBL_CopyToHost(Ci_host,&Ci[ic*Np],Np*sizeof(double));
 		double count_loc=0;
diff --git a/models/IonModel.h b/models/IonModel.h
index 4b91ac8f..faf9c4ff 100644
--- a/models/IonModel.h
+++ b/models/IonModel.h
@@ -34,7 +34,7 @@ public:
 	void Create();
 	void Initialize();
 	void Run(double *Velocity, double *ElectricField);
-    void getIonConcentration(DoubleArray &IonConcentration, const int ic);
+    void getIonConcentration(DoubleArray &IonConcentration, const size_t ic);
     void getIonConcentration_debug(int timestep);
     void DummyFluidVelocity();
     void DummyElectricField();
@@ -51,7 +51,7 @@ public:
     double fluidVelx_dummy,fluidVely_dummy,fluidVelz_dummy;
     double Ex_dummy,Ey_dummy,Ez_dummy;
 	
-	int number_ion_species;
+	size_t number_ion_species;
 	vector<int> BoundaryConditionInlet;
 	vector<int> BoundaryConditionOutlet;
     vector<double> IonDiffusivity;//User input unit [m^2/sec]
diff --git a/tests/DataAggregator.cpp b/tests/DataAggregator.cpp
index 7fcbfb0f..137ebd61 100644
--- a/tests/DataAggregator.cpp
+++ b/tests/DataAggregator.cpp
@@ -9,8 +9,8 @@
 using namespace std;
 
 int main(int argc, char **argv){
-	
-  printf("Aggregating block data into single file \n");
+
+	printf("Aggregating block data into single file \n");
 	unsigned int Bx,By,Bz;
 	uint64_t Nx,Ny,Nz;
 	uint64_t N;
@@ -22,25 +22,25 @@ int main(int argc, char **argv){
 
 	//Read the block size
 	if (argc>9){
-	  printf("Input arguments accepted \n");
-	Nx = atoi(argv[1]);
-	Ny = atoi(argv[2]);
-	Nz = atoi(argv[3]);
-	x0 = atoi(argv[4]);
-	y0 = atoi(argv[5]);
-	z0 = atoi(argv[6]);
-	NX = atol(argv[7]);
-	NY = atol(argv[8]);
-	NZ = atol(argv[9]);
-	printf("Size %i X %i X %i \n",NX,NY,NZ);
-	fflush(stdout);
+		printf("Input arguments accepted \n");
+		Nx = atoi(argv[1]);
+		Ny = atoi(argv[2]);
+		Nz = atoi(argv[3]);
+		x0 = atoi(argv[4]);
+		y0 = atoi(argv[5]);
+		z0 = atoi(argv[6]);
+		NX = atol(argv[7]);
+		NY = atol(argv[8]);
+		NZ = atol(argv[9]);
+		printf("Size %llu X %llu X %llu \n",(unsigned long long) NX, (unsigned long long) NY, (unsigned long long) NZ);
+		fflush(stdout);
 	}
 	else{
-	  printf("setting defaults \n");
-	  Nx=Ny=Nz=1024;
-	  x0=y0=z0=0;
-	  NX=NY=8640;
-	  NZ=6480;
+		printf("setting defaults \n");
+		Nx=Ny=Nz=1024;
+		x0=y0=z0=0;
+		NX=NY=8640;
+		NZ=6480;
 	}
 	//Bx = By = Bz = 9;
 	//Nx = Ny = Nz = 1024;
@@ -53,29 +53,28 @@ int main(int argc, char **argv){
 	if (By>8) By=8;
 	if (Bz>8) Bz=8;
 
-	printf("System size (output) is: %i x %i x %i \n",NX,NY,NZ);
-	printf("Block size (read) is: %i x %i x %i \n",Nx,Ny,Nz);
-	printf("Starting location (read) is: %i, %i, %i \n", x0,y0,z0);
-	printf("Block number (read): %i x %i x %i \n",Bx,By,Bz);
+	printf("System size (output) is: %llu x %llu x %llu \n",(unsigned long long) NX,(unsigned long long) NY, (unsigned long long) NZ);
+	printf("Block size (read) is: %llu x %llu x %llu \n",(unsigned long long) Nx,(unsigned long long) Ny,(unsigned long long) Nz);
+	printf("Starting location (read) is: %llu, %llu, %llu \n", (unsigned long long) x0,(unsigned long long) y0,(unsigned long long) z0);
+	printf("Block number (read): %llu x %llu x %llu \n",(unsigned long long) Bx,(unsigned long long) By,(unsigned long long) Bz);
 	fflush(stdout);
-	
+
 	// Filenames used
 	//char LocalRankString[8];
 	char LocalRankFilename[40];
 	char sx[2];
 	char sy[2];
 	char sz[2];
-	char tmpstr[10];
 
 	//sprintf(LocalRankString,"%05d",rank);
 	N = Nx*Ny*Nz;
 	N_full=NX*NY*NZ;
-	
+
 	char *id;
 	id = new char [N];
 	char *ID;
 	ID = new char [N_full];
-	
+
 	for (unsigned int  bz=0; bz<Bz; bz++){
 		for (unsigned int by=0; by<By; by++){
 			for (unsigned int bx=0; bx<Bx; bx++){
@@ -90,7 +89,7 @@ int main(int argc, char **argv){
 				FILE *IDFILE = fopen(LocalRankFilename,"rb");
 				readID=fread(id,1,N,IDFILE);
 				fclose(IDFILE);
-				printf("Loading data ... \n");
+				printf("Loading data: %zu bytes ... \n", readID);
 				// Unpack the data into the main array
 				for ( k=0;k<Nz;k++){
 					for ( j=0;j<Ny;j++){
@@ -99,7 +98,7 @@ int main(int argc, char **argv){
 							y = by*Ny + j;
 							z = bz*Nz + k;
 							if ( x<NX && y<NY && z<NZ){
-							  ID[z*NX*NY+y*NX+x] = id[k*Nx*Ny+j*Nx+i];
+								ID[z*NX*NY+y*NX+x] = id[k*Nx*Ny+j*Nx+i];
 							}
 						}
 					}
@@ -111,11 +110,11 @@ int main(int argc, char **argv){
 	// Compute porosity 
 	uint64_t count=0;
 	for (k=0; k<NZ; k++){
-	  for (j=0; j<NY; j++){
-	    for (i=0; i<NX; i++){
-	      if (ID[k*NX*NY+j*NX+i] < 215) count++; 
-	    }
-	  }
+		for (j=0; j<NY; j++){
+			for (i=0; i<NX; i++){
+				if (ID[k*NX*NY+j*NX+i] < 1) count++; 
+			}
+		}
 	}
 	printf("Porosity is %f \n",double(count)/double(NX*NY*NZ));
 
diff --git a/tests/GenerateSphereTest.cpp b/tests/GenerateSphereTest.cpp
index e28042d7..fa2fb885 100644
--- a/tests/GenerateSphereTest.cpp
+++ b/tests/GenerateSphereTest.cpp
@@ -137,20 +137,18 @@ inline void MorphOpen(DoubleArray SignDist, char *id, Domain &Dm, int nx, int ny
 	int Nx = nx;
 	int Ny = ny;
 	int Nz = nz;
+	int imin,jmin,kmin,imax,jmax,kmax;
+    
 	
 	double sw_old=1.0;
 	double sw_new=1.0; 
-    double sw_diff_old = 1.0;
-    double sw_diff_new = 1.0;
+	double sw_diff_old = 1.0;
+	double sw_diff_new = 1.0;
 
 	// Increase the critical radius until the target saturation is met
 	double deltaR=0.05; // amount to change the radius in voxel units
-	double Rcrit_old;
-	double Rcrit_new;
-
-	int imin,jmin,kmin,imax,jmax,kmax;
-    
-	Rcrit_new = maxdistGlobal;
+	double Rcrit_new = maxdistGlobal;
+	double Rcrit_old = maxdistGlobal;
     while (sw_new > SW)
     {
         sw_diff_old = sw_diff_new;
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index 51ffe96e..56b6d3dd 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -9,8 +9,6 @@
 #include "common/Utilities.h"
 #include "models/ColorModel.h"
 
-//#define WRE_SURFACES
-
 /*
  * Simulator for two-phase flow in porous media
  * James E. McClure 2013-2014
@@ -69,30 +67,69 @@ int main( int argc, char **argv )
         ColorModel.Initialize(); // initializing the model will set initial conditions for variables
        
         if (SimulationMode == "development"){
-            double MLUPS=0.0;
-            int timestep = 0;
-            /* flow adaptor keys to control */
-            auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" );
-        	int MAX_STEADY_TIME = flow_db->getWithDefault<int>( "max_steady_timesteps", 1000000 );
-        	int SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 100000 );
-
-            int ANALYSIS_INTERVAL = ColorModel.timestepMax;
+        	double MLUPS=0.0;
+        	int timestep = 0;
+		bool ContinueSimulation = true;
+        	int ANALYSIS_INTERVAL = ColorModel.timestepMax;
+        	/* flow adaptor keys to control */
+        	int SKIP_TIMESTEPS = 0;
+        	int MAX_STEADY_TIME = 1000000;
+		double ENDPOINT_THRESHOLD = 0.1;
+        	double FRACTIONAL_FLOW_INCREMENT = 0.05;
+        	if (ColorModel.db->keyExists( "FlowAdaptor" )){
+        		auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" );
+        		MAX_STEADY_TIME = flow_db->getWithDefault<int>( "max_steady_timesteps", 1000000 );
+        		SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 100000 );
+        		FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
+        	        ENDPOINT_THRESHOLD = flow_db->getWithDefault<double>( "endpoint_threshold", 0.1);
+        	}
         	if (ColorModel.analysis_db->keyExists( "analysis_interval" )){
         		ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar<int>( "analysis_interval" );
         	}
         	/* Launch the simulation */
         	FlowAdaptor Adapt(ColorModel);
         	runAnalysis analysis(ColorModel);
-        	
-            while (ColorModel.timestep < ColorModel.timestepMax){
-            	/* this will run steady points */
-            	timestep += MAX_STEADY_TIME;
-            	MLUPS = ColorModel.Run(timestep);
-            	if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
-            	Adapt.UpdateFractionalFlow(ColorModel);
-            	
-            	/* apply timestep skipping algorithm to accelerate steady-state */
+		while (ContinueSimulation){
+            	    /* this will run steady points */
+                    timestep += MAX_STEADY_TIME;
+		    MLUPS = ColorModel.Run(timestep);
+            	    if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
+            	    if (ColorModel.timestep > ColorModel.timestepMax){
+		      ContinueSimulation = false;
+		    }
+            	/* update the fractional flow by adding mass */
             	int skip_time = 0;
+            	timestep = ColorModel.timestep;
+    			double SaturationChange = 0.0;
+    			double volB = ColorModel.Averages->gwb.V; 
+    			double volA = ColorModel.Averages->gnb.V; 
+    			double initialSaturation = volB/(volA + volB);
+
+			double vA_x = ColorModel.Averages->gnb.Px/ColorModel.Averages->gnb.M; 
+			double vA_y = ColorModel.Averages->gnb.Py/ColorModel.Averages->gnb.M; 
+			double vA_z = ColorModel.Averages->gnb.Pz/ColorModel.Averages->gnb.M; 
+			double vB_x = ColorModel.Averages->gwb.Px/ColorModel.Averages->gwb.M; 
+			double vB_y = ColorModel.Averages->gwb.Py/ColorModel.Averages->gwb.M; 
+			double vB_z = ColorModel.Averages->gwb.Pz/ColorModel.Averages->gwb.M; 			
+			double speedA = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z);
+			double speedB = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
+			if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false;
+		if (ContinueSimulation){
+            	while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){
+            		timestep += ANALYSIS_INTERVAL;
+                	Adapt.UpdateFractionalFlow(ColorModel);
+                	MLUPS = ColorModel.Run(timestep);
+        			double volB = ColorModel.Averages->gwb.V; 
+        			double volA = ColorModel.Averages->gnb.V;
+        			SaturationChange = volB/(volA + volB) - initialSaturation;
+            		skip_time += ANALYSIS_INTERVAL;
+		}
+            	if (rank==0) printf("  *********************************************************************  \n");
+            	if (rank==0) printf("   Updated fractional flow with saturation change = %f  \n", SaturationChange);
+            	if (rank==0) printf("  *********************************************************************  \n");
+		}
+            	/* apply timestep skipping algorithm to accelerate steady-state */
+            	/* skip_time = 0;
             	timestep = ColorModel.timestep;
             	while (skip_time < SKIP_TIMESTEPS){
             		timestep += ANALYSIS_INTERVAL;
@@ -100,10 +137,9 @@ int main( int argc, char **argv )
             		Adapt.MoveInterface(ColorModel);
             		skip_time += ANALYSIS_INTERVAL;
             	}
-            	//Adapt.Flatten(ColorModel);
-
+            	*/
             }
-            ColorModel.WriteDebug();
+		  //ColorModel.WriteDebug();
         }
 
         else
@@ -112,6 +148,7 @@ int main( int argc, char **argv )
         PROFILE_STOP( "Main" );
         auto file  = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );
         auto level = db->getWithDefault<int>( "TimerLevel", 1 );
+	NULL_USE(level);
         PROFILE_SAVE( file, level );
         // ****************************************************
 
diff --git a/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp b/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp
index 19d99b9c..30ec4fec 100644
--- a/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp
+++ b/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp
@@ -59,6 +59,7 @@ int main( int argc, char **argv )
         PROFILE_STOP("Main");
         auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_freelee_SingleFluidBGK_simulator" );
         auto level = db->getWithDefault<int>( "TimerLevel", 1 );
+	NULL_USE(level);
         PROFILE_SAVE( file,level );
         // ****************************************************
 
diff --git a/tests/lbpm_freelee_simulator.cpp b/tests/lbpm_freelee_simulator.cpp
index cdba5082..e2ca2d64 100644
--- a/tests/lbpm_freelee_simulator.cpp
+++ b/tests/lbpm_freelee_simulator.cpp
@@ -83,6 +83,7 @@ int main( int argc, char **argv )
         PROFILE_STOP("Main");
         auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_freelee_simulator" );
         auto level = db->getWithDefault<int>( "TimerLevel", 1 );
+	NULL_USE(level);
         PROFILE_SAVE( file,level );
         // ****************************************************