diff --git a/tests/lbpm_random_pp.cpp b/tests/lbpm_random_pp.cpp
index 52b4f8f1..53c18eb1 100644
--- a/tests/lbpm_random_pp.cpp
+++ b/tests/lbpm_random_pp.cpp
@@ -12,93 +12,6 @@
 #include <Array.h>
 #include <Domain.h>
 
-inline void GenerateResidual(char *ID, int Nx, int Ny, int Nz, double Saturation)
-{
-	//.......................................................................
-	int i,j,k,n,Number,N;
-	int x,y,z,ii,jj,kk;
-	int sizeX,sizeY,sizeZ;
-	int *SizeX, *SizeY, *SizeZ;
-
-#ifdef NORANDOM
-	srand(10009);
-#else
-	srand(time(NULL));
-#endif
-//	float bin;
-	//.......................................................................
-	N = Nx*Ny*Nz;
-
-	int bin, binCount;
-	ifstream Dist("BlobSize.in");
-	Dist >> binCount;
-//	printf("Number of blob sizes: %i \n",binCount);
-	SizeX = new int [binCount];
-	SizeY = new int [binCount];
-	SizeZ = new int [binCount];
-	for (bin=0; bin<binCount; bin++){
-		Dist >> SizeX[bin];
-		Dist >> SizeY[bin];
-		Dist >> SizeZ[bin];
-	//	printf("Blob %i dimension: %i x %i x %i \n",bin, SizeX[bin], SizeY[bin], SizeZ[bin]);
-	}
-	Dist.close();
-	//.......................................................................
-//	cout << "Generating blocks... " << endl;
-	// Count for the total number of oil nodes
-	int count = 0;
-	// Count the total number of non-solid nodes
-	int total = 0;
-	for (i=0;i<N;i++){
-		if (ID[i] != 0) total++;
-	}
-
-	float sat = 0.f;
-	Number = 0;		// number of features
-	while (sat < Saturation){
-		Number++;
-		// Randomly generate a point in the domain
-		x = Nx*float(rand())/float(RAND_MAX);
-		y = Ny*float(rand())/float(RAND_MAX);
-		z = Nz*float(rand())/float(RAND_MAX);
-
-		bin = int(floor(binCount*float(rand())/float(RAND_MAX)));
-		sizeX = SizeX[bin];
-		sizeY = SizeY[bin];
-		sizeZ = SizeZ[bin];
-
-//		cout << "Sampling from bin no. " << floor(bin) << endl;
-//		cout << "Feature size is: " << sizeX << "x" << sizeY << "x" << sizeZ << endl;
-
-		for (k=z;k<z+sizeZ;k++){
-			for (j=y;j<y+sizeY;j++){
-				for (i=x;i<x+sizeX;i++){
-					// Identify nodes in the domain (periodic BC)
-					ii = i;
-					jj = j;
-					kk = k;
-					if (ii < 1)			ii+=(Nx-2);
-					if (jj < 1)			jj+=(Ny-2);
-					if (kk < 1)			kk+=(Nz-2);
-					if (!(ii < Nx-1))		ii-=(Nx-2);
-					if (!(jj < Ny-1))		jj-=(Ny-2);
-					if (!(kk < Nz-1))		kk-=(Nz-2);
-
-					n = kk*Nx*Ny+jj*Nx+ii;
-
-					if (ID[n] == 2){
-						ID[n] = 1;
-						count++;
-					}
-				}
-			}
-		}
-		sat = float(count)/total;
-	}
-	//.......................................................................
-}
-
-
 inline void FlipID(char *ID, int N)
 {
 	for (int n=0; n<N; n++){
@@ -182,6 +95,9 @@ int main(int argc, char **argv)
 
 	char LocalRankFilename[40];
 
+	int BoundaryCondition=0;
+	Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
+
 	nx+=2; ny+=2; nz+=2;
 	int N = nx*ny*nz;
     char *id;
@@ -194,22 +110,147 @@ int main(int argc, char **argv)
     fread(SignDist.get(),8,N,DIST);
     fclose(DIST);
 
+	int count,countGlobal,totalGlobal;
+    count = 0;
     for (int k=0; k<nz; k++){
         for (int j=0; j<ny; j++){
             for (int i=0; i<nx; i++){
-	      n = k*nx*ny+j*nx+i;
+            	n = k*nx*ny+j*nx+i;
             	if (SignDist(i,j,k) < 0.0)  id[n] = 0;
-            	else 			    id[n] = 2;
+            	else{
+            		id[n] = 2;
+            		count++;
+            	}
             }
         }
     }
+	// total Global is the number of nodes in the pore-space
+	MPI_Allreduce(&count,&totalGlobal,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+
+    Dm.CommInit(MPI_COMM_WORLD);
+    int iproc = Dm.iproc;
+    int jproc = Dm.jproc;
+    int kproc = Dm.kproc;
+
+	int bin, binCount;
+	ifstream Dist("BlobSize.in");
+	Dist >> binCount;
+	int *SizeX, *SizeY, *SizeZ;
+//	printf("Number of blob sizes: %i \n",binCount);
+	SizeX = new int [binCount];
+	SizeY = new int [binCount];
+	SizeZ = new int [binCount];
+	for (bin=0; bin<binCount; bin++){
+		Dist >> SizeX[bin];
+		Dist >> SizeY[bin];
+		Dist >> SizeZ[bin];
+	//	printf("Blob %i dimension: %i x %i x %i \n",bin, SizeX[bin], SizeY[bin], SizeZ[bin]);
+	}
+	Dist.close();
 
 	// Generate the residual NWP
 	if (rank==0) printf("Initializing with NWP saturation = %f \n",Saturation);
-	GenerateResidual(id,nx,ny,nz,Saturation);
+	//	GenerateResidual(id,nx,ny,nz,Saturation);
+
+	int x,y,z;
+	int sizeX,sizeY,sizeZ;
+	int ii,jj,kk;
+	int Nx = nx;
+	int Ny = ny;
+	int Nz = nz;
+	float sat = 0.f;
+	int Number = 0;		// number of features
+	while (sat < Saturation){
+		if (rank==0){
+			Number++;
+			// Randomly generate a point in the domain
+			x = (Nx-2)*nprocx*float(rand())/float(RAND_MAX);
+			y = (Ny-2)*nprocy*float(rand())/float(RAND_MAX);
+			z = (Nz-2)*nprocz*float(rand())/float(RAND_MAX);
+
+			bin = int(floor(binCount*float(rand())/float(RAND_MAX)));
+			sizeX = SizeX[bin];
+			sizeY = SizeY[bin];
+			sizeZ = SizeZ[bin];
+		}
+		MPI_Bcast(&x,1,MPI_INT,0,MPI_COMM_WORLD);
+		MPI_Bcast(&y,1,MPI_INT,0,MPI_COMM_WORLD);
+		MPI_Bcast(&z,1,MPI_INT,0,MPI_COMM_WORLD);
+		MPI_Bcast(&sizeX,1,MPI_INT,0,MPI_COMM_WORLD);
+		MPI_Bcast(&sizeY,1,MPI_INT,0,MPI_COMM_WORLD);
+		MPI_Bcast(&sizeZ,1,MPI_INT,0,MPI_COMM_WORLD);
+
+		//if (rank==0) printf("Broadcast block at %i,%i,%i \n",x,y,z);
+
+		for (k=z;k<z+sizeZ;k++){
+			for (j=y;j<y+sizeY;j++){
+				for (i=x;i<x+sizeX;i++){
+
+					// Identify nodes in the domain (periodic BC)
+					ii = i;
+					jj = j;
+					kk = k;
+
+					if (ii>nprocx*(Nx-2)) ii-=nprocx*(Nx-2);
+					if (jj>nprocy*(Ny-2)) jj-=nprocy*(Ny-2);
+					if (kk>nprocz*(Nz-2)) kk-=nprocz*(Nz-2);
+
+					// Check if this is in the subdomain
+					if (ii < (iproc+1)*(Nx-2) && jj < (jproc+1)*(Ny-2) && kk < (kproc+1)*(Nz-2) &&
+							ii  > iproc*(Nx-2) && jj > jproc*(Ny-2) && kk > kproc*(Nz-2) ){
+
+						// Map from global to local coordinates
+						ii -= iproc*(Nx-2);
+						jj -= jproc*(Ny-2);
+						kk -= kproc*(Nz-2);
+
+						n = kk*Nx*Ny+jj*Nx+ii;
+
+						if (id[n] == 2){
+							id[n] = 1;
+							//count++;
+						}
+
+					}
+				}
+			}
+		}
+		count = 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++){
+				  n=k*Nx*Ny+j*Nx+i;
+					if (id[n] == 1){
+						count++;
+					}
+				}
+			}
+		}
+		MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+		sat = float(countGlobal)/totalGlobal;
+		//if (rank==0) printf("New count=%i\n",countGlobal);
+		//if (rank==0) printf("New saturation=%f\n",sat);
+	}
 
 	if (InitialWetting == 1)	FlipID(id,nx*ny*nz);
 
+
+		count = 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++){
+				  n=k*Nx*Ny+j*Nx+i;
+
+					if (id[n] == 1){
+						count++;
+					}
+				}
+			}
+		}
+		MPI_Allreduce(&count,&countGlobal,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+		sat = float(countGlobal)/totalGlobal;
+	if (rank==0) printf("Final saturation=%f\n",sat);
+
     sprintf(LocalRankFilename,"ID.%05i",rank);
     FILE *ID = fopen(LocalRankFilename,"wb");
     fwrite(id,1,N,ID);