diff --git a/common/TwoPhase.h b/common/TwoPhase.h
index 6f8a3fb2..87ae7395 100644
--- a/common/TwoPhase.h
+++ b/common/TwoPhase.h
@@ -2,6 +2,55 @@
 #include "pmmc.h"
 #include "Domain.h"
 #include "Communication.h"
+#include <vector>
+
+#define BLOB_AVG_COUNT 26
+
+struct BlobContainer{
+
+    BlobContainer(){
+      NBLOBS=0;
+    }
+    ~BlobContainer(){
+    }
+    void Set(int size){
+      NBLOBS=size;
+      Data.resize(size*BLOB_AVG_COUNT);
+    }
+    int NBLOBS;
+    std::vector<int> Data;
+    
+    // if modified -- make sure to adjust COUNT so that
+    // there is enough memory to save all the averages
+    double Vn(int IDX){return Data[BLOB_AVG_COUNT*IDX];}
+    double pan(int IDX){return Data[BLOB_AVG_COUNT*IDX+1];}
+    double awn(int IDX){return Data[BLOB_AVG_COUNT*IDX+2];}
+    double ans(int IDX){return Data[BLOB_AVG_COUNT*IDX+3];}
+    double Jwn(int IDX){return Data[BLOB_AVG_COUNT*IDX+4];}
+    double Kwn(int IDX){return Data[BLOB_AVG_COUNT*IDX+5];}
+    double lwns(int IDX){return Data[BLOB_AVG_COUNT*IDX+6];}
+    double cwns(int IDX){return Data[BLOB_AVG_COUNT*IDX+7];}
+    double vanx(int IDX){return Data[BLOB_AVG_COUNT*IDX+8];}
+    double vany(int IDX){return Data[BLOB_AVG_COUNT*IDX+9];}
+    double vanz(int IDX){return Data[BLOB_AVG_COUNT*IDX+10];}
+    double vawnx(int IDX){return Data[BLOB_AVG_COUNT*IDX+11];}
+    double vawny(int IDX){return Data[BLOB_AVG_COUNT*IDX+12];}
+    double vawnz(int IDX){return Data[BLOB_AVG_COUNT*IDX+13];}
+    double Gwnxx(int IDX){return Data[BLOB_AVG_COUNT*IDX+14];}
+    double Gwnyy(int IDX){return Data[BLOB_AVG_COUNT*IDX+15];}
+    double Gwnzz(int IDX){return Data[BLOB_AVG_COUNT*IDX+16];}
+    double Gwnxy(int IDX){return Data[BLOB_AVG_COUNT*IDX+17];}
+    double Gwnxz(int IDX){return Data[BLOB_AVG_COUNT*IDX+18];}
+    double Gwnyz(int IDX){return Data[BLOB_AVG_COUNT*IDX+19];}
+    double Gnsxx(int IDX){return Data[BLOB_AVG_COUNT*IDX+20];}
+    double Gnsyy(int IDX){return Data[BLOB_AVG_COUNT*IDX+22];}
+    double Gnszz(int IDX){return Data[BLOB_AVG_COUNT*IDX+23];}
+    double Gnsxy(int IDX){return Data[BLOB_AVG_COUNT*IDX+23];}
+    double Gnsxz(int IDX){return Data[BLOB_AVG_COUNT*IDX+24];}
+    double Gnsyz(int IDX){return Data[BLOB_AVG_COUNT*IDX+25];}
+
+};
+
 
 class TwoPhase{
 
@@ -113,6 +162,7 @@ public:
 	DoubleArray Vel_x;		// Velocity
 	DoubleArray Vel_y;
 	DoubleArray Vel_z;
+	BlobContainer BlobAverages;
 	//...........................................................................
 	TwoPhase(Domain &dm) : Dm(dm){
 		Nx=dm.Nx; Ny=dm.Ny; Nz=dm.Nz;
@@ -443,8 +493,17 @@ void TwoPhase::ComputeLocal(){
 
 void TwoPhase::ComputeLocalBlob(){
         int i,j,k,n,label;
+	int nblobs_global;
 	double delphi;
 	int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};
+        // get the maximum label locally -- then compute number of global blobs
+	label=0; 
+	for (n=0; n<Nx*Ny*Nz; n++){
+	  if (label < BlobLabel.data[n]) label = BlobLabel.data[n];
+	}
+	MPI_Allreduce(&label,&nblobs_global,1,MPI_INT,MPI_MAX,Dm.Comm);
+	if (Dm.rank==0) printf("Number of blobs is %i \n",nblobs_global);
+	// Perform averaging
 
 	for (int c=0;c<ncubes;c++){
 		// Get cube from the list
@@ -653,15 +712,15 @@ void TwoPhase::PrintAll(int timestep){
 
 inline int TwoPhase::GetCubeLabel(int i, int j, int k){
 	int label;
+
 	label=BlobLabel(i,j,k);
 	label=max(label,BlobLabel(i+1,j,k));
-	label=max(label,BlobLabel(i,j+1,k));
+      	label=max(label,BlobLabel(i,j+1,k));
 	label=max(label,BlobLabel(i+1,j+1,k));
 	label=max(label,BlobLabel(i,j,k+1));
 	label=max(label,BlobLabel(i+1,j,k+1));
 	label=max(label,BlobLabel(i,j+1,k+1));
 	label=max(label,BlobLabel(i+1,j+1,k+1));
-
+	
 	return label;
 }
-
diff --git a/tests/lbpm_BlobAnalysis.cpp b/tests/lbpm_BlobAnalysis.cpp
index 64474f5c..2477d9f1 100644
--- a/tests/lbpm_BlobAnalysis.cpp
+++ b/tests/lbpm_BlobAnalysis.cpp
@@ -37,50 +37,6 @@ inline void ReadBlobFile(char *FILENAME, int *Data, int N)
     
 }
 
-struct BlobTwoPhase{
-    int COUNT; // number of averages to compute for each blob
-    BlobTwoPhase(int size){
-    	COUNT=26;
-        NBLOBS=size;
-        Data = new double [size*COUNT];
-    }
-    ~BlobTwoPhase(){
-        delete [] Data;
-    }
-    int NBLOBS;
-    double *Data;
-    
-    // if modified -- make sure to adjust COUNT so that
-    // there is enough memory to save all the averages
-    double Vn(int IDX){return Data[COUNT*IDX];}
-    double pan(int IDX){return Data[COUNT*IDX+1];}
-    double awn(int IDX){return Data[COUNT*IDX+2];}
-    double ans(int IDX){return Data[COUNT*IDX+3];}
-    double Jwn(int IDX){return Data[COUNT*IDX+4];}
-    double Kwn(int IDX){return Data[COUNT*IDX+5];}
-    double lwns(int IDX){return Data[COUNT*IDX+6];}
-    double cwns(int IDX){return Data[COUNT*IDX+7];}
-    double vanx(int IDX){return Data[COUNT*IDX+8];}
-    double vany(int IDX){return Data[COUNT*IDX+9];}
-    double vanz(int IDX){return Data[COUNT*IDX+10];}
-    double vawnx(int IDX){return Data[COUNT*IDX+11];}
-    double vawny(int IDX){return Data[COUNT*IDX+12];}
-    double vawnz(int IDX){return Data[COUNT*IDX+13];}
-    double Gwnxx(int IDX){return Data[COUNT*IDX+14];}
-    double Gwnyy(int IDX){return Data[COUNT*IDX+15];}
-    double Gwnzz(int IDX){return Data[COUNT*IDX+16];}
-    double Gwnxy(int IDX){return Data[COUNT*IDX+17];}
-    double Gwnxz(int IDX){return Data[COUNT*IDX+18];}
-    double Gwnyz(int IDX){return Data[COUNT*IDX+19];}
-    double Gnsxx(int IDX){return Data[COUNT*IDX+20];}
-    double Gnsyy(int IDX){return Data[COUNT*IDX+22];}
-    double Gnszz(int IDX){return Data[COUNT*IDX+23];}
-    double Gnsxy(int IDX){return Data[COUNT*IDX+23];}
-    double Gnsxz(int IDX){return Data[COUNT*IDX+24];}
-    double Gnsyz(int IDX){return Data[COUNT*IDX+25];}
-
-};
-
 int main(int argc, char **argv)
 {
 	//*****************************************
@@ -97,18 +53,19 @@ int main(int argc, char **argv)
 	int Nx,Ny,Nz,N,nspheres;
 	double Lx,Ly,Lz;
 
-	int BC;	// type of boundary condition applied: 0-periodic, 1-pressure/velocity
-	int nblobs_global; 	// number of blobs in the global system
+	int BC=0;	// type of boundary condition applied: 0-periodic, 1-pressure/velocity
+	int nblobs_global=0; 	// number of blobs in the global system
     
     // Get the global number of blobs from arguments
-    if (argc > 1){
+	/*    if (argc > 1){
         nblobs_global = atoi(argv[1]);
         if (rank==0)    printf("Number of global blobs is: %i \n",nblobs_global);
     }
     else{
         ERROR("Number of blobs was not specified");
     }
-
+	*/
+	
 	int *CubeList;
 
 	if (rank==0){
@@ -145,7 +102,9 @@ int main(int argc, char **argv)
 	//.................................................
 	Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
 	TwoPhase Averages(Dm);
-	BlobTwoPhase BlobAverages(nblobs_global);
+	//	BlobTwoPhase BlobAverages(nblobs_global);
+	Nx+=2;Ny+=2;Nz+=2;
+	N=Nx*Ny*Nz; // number of lattice points
 	//.......................................................................
 	// Filenames used
 	char LocalRankString[8];
@@ -167,6 +126,7 @@ int main(int argc, char **argv)
     ReadBlobFile(LocalRankFilename, Averages.BlobLabel.data, N);
     MPI_Barrier(MPI_COMM_WORLD);
     if (rank == 0) cout << "BlobLabel set." << endl;
+
     //.......................................................................
 	//copies of data needed to perform checkpointing from cpu
 	double *Den, *DistEven, *DistOdd;
@@ -180,6 +140,7 @@ int main(int argc, char **argv)
 	MPI_Barrier(MPI_COMM_WORLD);
 	//.........................................................................
 	// Populate the arrays needed to perform averaging
+	if (rank==0) printf("Populate arrays \n");
     for (int n=0; n<Nx*Ny*Nz; n++){
         double phi,da,db,press,vx,vy,vz;
         double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
@@ -219,16 +180,13 @@ int main(int argc, char **argv)
         Averages.Vel_y.data[n]=vy;
         Averages.Vel_z.data[n]=vz;
     }
+    delete [] DistEven;
+    delete [] DistOdd;
+    printf("Ready for averaging, rank=%i \n",rank);
     int label;
-	for (int k=1;k<Nz-1;k++){
-		for (int j=1;j<Ny-1;j++){
-			for (int i=1;i<Nx-1;i++){
-			  // Assign the label for the cube
-			  label = Averages.GetCubeLabel(i,j,k);
-			  
-			}
-		}
-	}
+    Averages.ComputeLocalBlob();
+    
+    printf("Exit, rank=%i \n",rank);
 
 	/*	Averages.Initialize();
 	Averages.ComputeDelPhi();