diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp
index 64d6ff62..414fa14f 100644
--- a/analysis/morphology.cpp
+++ b/analysis/morphology.cpp
@@ -344,8 +344,8 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, d
 	int nprocz = Dm->nprocz();
 	int rank = Dm->rank();
 	
-	DoubleArray phase(Nx,Ny,Nz);
-	IntArray phase_label(Nx,Ny,Nz);
+	DoubleArray phase(nx,ny,nz);
+	IntArray phase_label(nx,ny,nz);
 
 	int n;
 	double final_void_fraction;
@@ -569,11 +569,10 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, d
 		//......................................................................................
 		MPI_Allreduce(&LocalNumber,&GlobalNumber,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
 		
-		
-		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;
+		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){
 						phase(i,j,k) = 1.0;
 					}
@@ -586,13 +585,13 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, d
 		// Extract only the connected part
 		BlobIDstruct new_index;
 		double vF=0.0; double vS=0.0;
-		ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,Dm->rank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm);
+		ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm);
 		MPI_Barrier(comm);
 		
-		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;
+		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 && phase_label(i,j,k) > 1){
 						id[n] = 2;
 					}
@@ -601,10 +600,10 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, d
 		}
 
 		count = 0.f;
-		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;
+		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] == 2){
 						count+=1.0;
 					}