diff --git a/common/Domain.cpp b/common/Domain.cpp
index fc797a8e..9cf9a1ca 100644
--- a/common/Domain.cpp
+++ b/common/Domain.cpp
@@ -1433,3 +1433,94 @@ void Domain::ReadFromFile(const std::string& Filename,const std::string& Datatyp
 	//Comm.barrier();
 	MPI_Barrier(Comm);
 }
+
+void Domain::AggregateLabels( const std::string& filename, DoubleArray &UserData ){
+	
+	int nx = Nx;
+	int ny = Ny;
+	int nz = Nz;
+	
+	int npx = nprocx();
+	int npy = nprocy();
+	int npz = nprocz();
+	
+	int ipx = iproc();
+	int ipy = jproc();
+	int ipz = kproc();		
+	
+	int nprocs = nprocx()*nprocy()*nprocz();
+		
+	int full_nx = npx*(nx-2);
+	int full_ny = npy*(ny-2);
+	int full_nz = npz*(nz-2);
+	int local_size = (nx-2)*(ny-2)*(nz-2);
+	unsigned long int full_size = long(full_nx)*long(full_ny)*long(full_nz);
+	
+	double *LocalID;
+	LocalID = new double [local_size];
+		
+	//printf("aggregate labels: local size=%i, global size = %i",local_size, full_size);
+	// assign the ID for the local sub-region
+	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 = k*nx*ny+j*nx+i;
+				double local_id_val = UserData(i,j,k); 
+				LocalID[(k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1] = local_id_val;
+			}
+		}
+	}
+	MPI_Barrier(Comm);
+
+	// populate the FullID 
+	if (rank() == 0){
+		double *FullID;
+		FullID = new double [full_size];
+		// first handle local ID for rank 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 x = i-1;
+					int y = j-1;
+					int z = k-1;
+					int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
+					unsigned long int n_full = z*long(full_nx)*long(full_ny) + y*long(full_nx) + x;
+					FullID[n_full] = LocalID[n_local];
+				}
+			}
+		}
+		// next get the local ID from the other ranks
+		for (int rnk = 1; rnk<nprocs; rnk++){
+			ipz = rnk / (npx*npy);
+			ipy = (rnk - ipz*npx*npy) / npx;
+			ipx = (rnk - ipz*npx*npy - ipy*npx); 
+			//printf("ipx=%i ipy=%i ipz=%i\n", ipx, ipy, ipz);
+			int tag = 15+rnk;
+			MPI_Recv(LocalID,local_size,MPI_DOUBLE,rnk,tag,Comm,MPI_STATUS_IGNORE);
+			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 x = i-1 + ipx*(nx-2);
+						int y = j-1 + ipy*(ny-2);
+						int z = k-1 + ipz*(nz-2);
+						int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
+						unsigned long int n_full = z*long(full_nx)*long(full_ny) + y*long(full_nx) + x;
+						FullID[n_full] = LocalID[n_local];
+					}
+				}
+			}
+		}
+		// write the output
+		FILE *OUTFILE = fopen(filename.c_str(),"wb");
+		fwrite(FullID,8,full_size,OUTFILE);
+		fclose(OUTFILE);
+	}
+	else{
+		// send LocalID to rank=0
+		int tag = 15+ rank();
+		int dstrank = 0;
+		MPI_Send(LocalID,local_size,MPI_DOUBLE,dstrank,tag,Comm);
+	}
+	MPI_Barrier(Comm);
+
+}
diff --git a/common/Domain.h b/common/Domain.h
index 9aa3e0ad..0b9a01ab 100755
--- a/common/Domain.h
+++ b/common/Domain.h
@@ -183,6 +183,7 @@ public: // Public variables (need to create accessors instead)
     void CommInit(); 
     int PoreCount();
     void AggregateLabels( const std::string& filename );
+    void AggregateLabels( const std::string& filename, DoubleArray &UserData );
 
 private:
 
diff --git a/models/IonModel.cpp b/models/IonModel.cpp
index bc58a649..0d977416 100644
--- a/models/IonModel.cpp
+++ b/models/IonModel.cpp
@@ -794,11 +794,14 @@ void ScaLBL_IonModel::getIonConcentration(int timestep){
         ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
         IonConcentration_LB_to_Phys(PhaseField);
 
-        FILE *OUTFILE;
-        sprintf(LocalRankFilename,"Ion%02i_Time_%i.%05i.raw",ic+1,timestep,rank);
-        OUTFILE = fopen(LocalRankFilename,"wb");
-        fwrite(PhaseField.data(),8,N,OUTFILE);
-        fclose(OUTFILE);
+        //FILE *OUTFILE;
+        //sprintf(LocalRankFilename,"Ion%02i_Time_%i.%05i.raw",ic+1,timestep,rank);
+        //OUTFILE = fopen(LocalRankFilename,"wb");
+        //fwrite(PhaseField.data(),8,N,OUTFILE);
+        //fclose(OUTFILE);
+
+        sprintf(OutputFilename,"Ion%02i_Time_%i.raw",ic+1,timestep);
+        Mask->AggregateLabels(OutputFilename,PhaseField);
     }
 
 }
diff --git a/models/IonModel.h b/models/IonModel.h
index 995898bd..0f527c3a 100644
--- a/models/IonModel.h
+++ b/models/IonModel.h
@@ -85,6 +85,7 @@ private:
     char LocalRankString[8];
     char LocalRankFilename[40];
     char LocalRestartFile[40];
+    char OutputFilename[200];
    
     //int rank,nprocs;
     void LoadParams(std::shared_ptr<Database> db0);