diff --git a/pmmc/BlobAnalysis.cpp b/pmmc/BlobAnalysis.cpp
new file mode 100644
index 00000000..65b17201
--- /dev/null
+++ b/pmmc/BlobAnalysis.cpp
@@ -0,0 +1,102 @@
+// Sequential blob analysis 
+// Reads parallel simulation data and performs connectivity analysis
+// and averaging on a blob-by-blob basis
+// James E. McClure 2014
+
+#include <iostream>
+#include <math.h>
+#include "pmmc.h"
+#include "Domain.h"
+
+using namespace std;
+
+inline void ReadFromAllRanks(){
+	
+}
+
+int main(int argc, char **argv)
+{
+	//.......................................................................
+	int nprocx,nprocy,nprocz,nprocs;
+	int Nx, Ny, Nz;
+	int nx,ny,nz;
+	int nspheres;
+	double Lx,Ly,Lz;
+	//.......................................................................
+	int i,j,k;
+	//.......................................................................
+	// Reading the domain information file
+	//.......................................................................
+	ifstream domain("Domain.in");
+	domain >> nprocx;
+	domain >> nprocy;
+	domain >> nprocz;
+	domain >> nx;
+	domain >> ny;
+	domain >> nz;
+	domain >> nspheres;
+	domain >> Lx;
+	domain >> Ly;
+	domain >> Lz;
+	//.......................................................................
+
+	nprocs = nprocx*nproy*nprocz;
+	printf("Number of MPI ranks: %i \n", nprocs);
+	Nx = (nx-2)*nprocx;
+	Ny = (ny-2)*nprocy;
+	Nz = (nz-2)*nprocz;
+	printf("Full domain size: %i x %i x %i  \n", Nx,Ny,Nz);
+	
+	// Filenames used
+	char LocalRankString[8];
+	char LocalRankFilename[40];
+	char LocalRestartFile[40];
+	char BaseFilename[20];
+	char tmpstr[10];
+	sprintf(LocalRankString,"%05d",rank);
+	sprintf(BaseFilename,"%s","Phase");
+	sprintf(LocalRankFilename,"%s%s",BaseFilename,LocalRankString);
+
+	                   
+	IntArray LocalBlobID(Nx,Ny,Nz);
+	DoubleArray Phase(Nx,Ny,Nz);
+	DoubleArray SignDist(Nx,Ny,Nz);
+	DoubleArray Press(Nx,Ny,Nz);
+	DoubleArray Vel_x(Nx,Ny,Nz);			// Velocity
+	DoubleArray Vel_y(Nx,Ny,Nz);
+	DoubleArray Vel_z(Nx,Ny,Nz);
+	DoubleArray MeanCurvature(Nx,Ny,Nz);
+	DoubleArray GaussCurvature(Nx,Ny,Nz);
+	DoubleArray SignDist_x(Nx,Ny,Nz);		// Gradient of the signed distance
+	DoubleArray SignDist_y(Nx,Ny,Nz);
+	DoubleArray SignDist_z(Nx,Ny,Nz);
+	DoubleArray Phase_x(Nx,Ny,Nz);			// Gradient of the phase indicator field
+	DoubleArray Phase_y(Nx,Ny,Nz);
+	DoubleArray Phase_z(Nx,Ny,Nz);
+	
+	double *Temp;
+	Temp = new double [3*nx*ny*nz];
+	
+	// read the files and populate main arrays
+	for (int proc=0; proc<nprocs; proc++){
+		
+		sprintf(LocalRankString,"%05d",proc);
+		sprintf(LocalRankFilename,"%s%s",BaseFilename,LocalRankString);
+		printf("Reading file %s \n",LocalRankFilename);
+		ReadBinaryFile(LocalRankFilename, Temp, nx*ny*nz);	
+		
+		for (k=1; k<nz-1; k++){
+			for (j=1; j<ny-1; j++){
+				for (i=1; i<nz-1; i++){
+					n = k*nx*ny+j*nx+i;
+					Phase(i-1,j-1,k-1) = Temp[n];
+				}
+			}
+		}
+		
+	}
+	
+	printf("Read %i ranks of %s \n",nprocs,BaseFilename);
+	
+}
+
diff --git a/tests/lb2_Color_wia_mpi.cpp b/tests/lb2_Color_wia_mpi.cpp
index 51fe6dcc..a29650ab 100644
--- a/tests/lb2_Color_wia_mpi.cpp
+++ b/tests/lb2_Color_wia_mpi.cpp
@@ -17,7 +17,7 @@
 
 //#define CBUB
 //#define WRITE_SURFACES
-#define USE_EXP_CONTACT_ANGLE
+//#define USE_EXP_CONTACT_ANGLE
 
 using namespace std;
 
@@ -1196,8 +1196,8 @@ int main(int argc, char **argv)
 	int kstart,kfinish;
 	kstart = 1;
 	kfinish = Nz-1;
-//	if (pBC && kproc==0)		kstart = 4;
-//	if (pBC && kproc==nprocz-1)	kfinish = Nz-4;
+	if (pBC && kproc==0)		kstart = 4;
+	if (pBC && kproc==nprocz-1)	kfinish = Nz-4;
 	for (k=kstart; k<kfinish; k++){
 		for (j=1; j<Ny-1; j++){
 			for (i=1; i<Nx-1; i++){
@@ -1208,6 +1208,7 @@ int main(int argc, char **argv)
 			}
 		}
 	}
+	ncubes = nc;
 
 
 	int logcount = 0; // number of surface write-outs