diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index a9d6f910..9b69441e 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -112,9 +112,9 @@ void SubPhase::Write(int timestep)
 		fprintf(SUBPHASE,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn); 
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.p, gwd.p, gnc.p, gnd.p);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
-		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx);
-		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny);
-		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Px, gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Py, gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny);
+		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Pz, gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.V, gwc.A, gwc.H, gwc.X);
 		fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
@@ -509,10 +509,6 @@ void SubPhase::Full(){
 	double vol_wc_bulk = 0.0;
 	double vol_nd_bulk = 0.0;
 	double vol_wd_bulk = 0.0;
-	double count_wc = 0.0;
-	double count_nc = 0.0;
-	double count_wd = 0.0;
-	double count_nd = 0.0;
 	for (k=kmin; k<kmax; k++){
 		for (j=jmin; j<Ny-1; j++){
 			for (i=imin; i<Nx-1; i++){
@@ -537,28 +533,27 @@ void SubPhase::Full(){
 					}
 					else if ( phi > 0.0){
 						if (morph_n->label(i,j,k) > 0 ){
-							count_nd += 1.0;
+							vol_nd_bulk += 1.0;
 							nd.p += Pressure(n);
 						}
 						else{
-							count_nc += 1.0;
+							vol_nc_bulk += 1.0;
 							nc.p += Pressure(n);
 						}
 					}
 					else{
 						// water region
 						if (morph_w->label(i,j,k) > 0 ){
-							count_wd += 1.0;
+							vol_wd_bulk += 1.0;
 							wd.p += Pressure(n);
 						}
 						else{
-							count_wc += 1.0;
+							vol_wc_bulk += 1.0;
 							wc.p += Pressure(n);
 						}
 					}
 					if ( phi > 0.0){
 						if (morph_n->label(i,j,k) > 0 ){
-							vol_nd_bulk += 1.0;
 							nd.M += nA*rho_n;						
 							nd.Px += nA*rho_n*ux;
 							nd.Py += nA*rho_n*uy;
@@ -566,7 +561,6 @@ void SubPhase::Full(){
 							nd.K += nA*rho_n*(ux*ux + uy*uy + uz*uz);
 						}
 						else{
-							vol_nc_bulk += 1.0;
 							nc.M += nA*rho_n;						
 							nc.Px += nA*rho_n*ux;
 							nc.Py += nA*rho_n*uy;
@@ -577,7 +571,6 @@ void SubPhase::Full(){
 					else{
 						// water region
 						if (morph_w->label(i,j,k) > 0 ){
-							vol_wd_bulk += 1.0;
 							wd.M += nB*rho_w;						
 							wd.Px += nB*rho_w*ux;
 							wd.Py += nB*rho_w*uy;
@@ -585,7 +578,6 @@ void SubPhase::Full(){
 							wd.K += nB*rho_w*(ux*ux + uy*uy + uz*uz);
 						}
 						else{
-							vol_wc_bulk += 1.0;
 							wc.M += nB*rho_w;						
 							wc.Px += nB*rho_w*ux;
 							wc.Py += nB*rho_w*uy;
@@ -597,14 +589,6 @@ void SubPhase::Full(){
 			}
 		}
 	}
-	count_wc=sumReduce( Dm->Comm, count_wc);
-	count_nc=sumReduce( Dm->Comm, count_nc);
-	count_wd=sumReduce( Dm->Comm, count_wd);
-	count_nd=sumReduce( Dm->Comm, count_nd);
-	gnd.p=sumReduce( Dm->Comm, nd.p) / count_nd;
-	gwd.p=sumReduce( Dm->Comm, wd.p) / count_wd;
-	gnc.p=sumReduce( Dm->Comm, nc.p) / count_nc;
-	gwc.p=sumReduce( Dm->Comm, wc.p) / count_wc;
 
 	gnd.M=sumReduce( Dm->Comm, nd.M);
 	gnd.Px=sumReduce( Dm->Comm, nd.Px);
diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index df7f478b..cd95a4b8 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -890,6 +890,7 @@ void runAnalysis::basic( std::shared_ptr<Database> db, SubPhase &Averages, const
     int timestep = db->getWithDefault<int>( "timestep", 0 );
 
     // Check which analysis steps we need to perform
+    int timestep = db->getWithDefault<int>( "timestep", 0 );
     auto type = computeAnalysisType( timestep );
     if ( type == AnalysisType::AnalyzeNone )
         return;