diff --git a/src_compressible_stag/main_driver.cpp b/src_compressible_stag/main_driver.cpp
index 10a92730..2789cbb3 100644
--- a/src_compressible_stag/main_driver.cpp
+++ b/src_compressible_stag/main_driver.cpp
@@ -1147,9 +1147,14 @@ void main_driver(const char* argv)
         }
         statsCount++;
         if (step%100 == 0) {
-            amrex::Print() << "Mean Density: " << ComputeSpatialMean(cu, 0) << " Mean Momentum (x):" << ComputeSpatialMean(cumom[0], 0) 
-                           <<  " Mean Momentum (y):" << ComputeSpatialMean(cumom[1], 0) <<  " Mean Momentum (z):" 
-                           << ComputeSpatialMean(cumom[2], 0) <<" Mean Energy:" << ComputeSpatialMean(cu, 4) << "\n";
+            amrex::Print() << "Mean Rho: "      << ComputeSpatialMean(cu, 0) 
+                           << " Mean Temp.:"    << ComputeSpatialMean(prim, 4)
+                           << " Mean Press.:"   << ComputeSpatialMean(prim, 5)
+                           << " Mean Mom. (x):" << ComputeSpatialMean(cumom[0], 0) 
+                           << " Mean Mom. (y):" << ComputeSpatialMean(cumom[1], 0) 
+                           << " Mean Mom. (z):" << ComputeSpatialMean(cumom[2], 0) 
+                           << " Mean En.:"      << ComputeSpatialMean(cu, 4) 
+                           << "\n";
         }
 
         // turbulence outputs
diff --git a/src_compressible_stag/timeStepStag.cpp b/src_compressible_stag/timeStepStag.cpp
index ac7ec926..8f0e93d1 100644
--- a/src_compressible_stag/timeStepStag.cpp
+++ b/src_compressible_stag/timeStepStag.cpp
@@ -761,10 +761,10 @@ void RK3stepStag(MultiFab& cu,
                                                     + grav[1]*(mompy(i,j+1,k)+mompy(i,j,k))
                                                     + grav[2]*(mompz(i,j,k+1)+mompz(i,j,k)) );
             if (turbForcing == 2) {
-                cup2_fab(i,j,k,4) += 0.5 * dt * (  (turbvf_x_o(i+1,j,k) + 0.25*(turbvf_x(i+1,j,k)-turbvf_x_o(i+1,j,k)))*mompx(i+1,j,k) + mompx(i,j,k)*(turbvf_x_o(i,j,k) + 0.25*(turbvf_x(i,j,k)-turbvf_x_o(i,j,k)))
+                cup2_fab(i,j,k,4) += 0.5 * 0.25 * dt * (  (turbvf_x_o(i+1,j,k) + 0.25*(turbvf_x(i+1,j,k)-turbvf_x_o(i+1,j,k)))*mompx(i+1,j,k) + mompx(i,j,k)*(turbvf_x_o(i,j,k) + 0.25*(turbvf_x(i,j,k)-turbvf_x_o(i,j,k)))
                                                  + (turbvf_y_o(i,j+1,k) + 0.25*(turbvf_y(i,j+1,k)-turbvf_y_o(i,j+1,k)))*mompy(i,j+1,k) + mompy(i,j,k)*(turbvf_y_o(i,j,k) + 0.25*(turbvf_y(i,j,k)-turbvf_y_o(i,j,k)))
                                                  + (turbvf_z_o(i,j,k+1) + 0.25*(turbvf_z(i,j,k+1)-turbvf_z_o(i,j,k+1)))*mompz(i,j,k+1) + mompz(i,j,k)*(turbvf_z_o(i,j,k) + 0.25*(turbvf_z(i,j,k)-turbvf_z_o(i,j,k))) );
-                cup2_fab(i,j,k,4) -= dt * energyp_in; // remove the average input energy from random turbulent forcing
+                cup2_fab(i,j,k,4) -= 0.25 * dt * energyp_in; // remove the average input energy from random turbulent forcing
             }
         });
     }
@@ -1035,10 +1035,10 @@ void RK3stepStag(MultiFab& cu,
                                                     + grav[1]*(momp2y(i,j+1,k)+momp2y(i,j,k))
                                                     + grav[2]*(momp2z(i,j,k+1)+momp2z(i,j,k)) );
             if (turbForcing == 2) {
-                cu_fab(i,j,k,4) += 0.5 * dt * (  (turbvf_x_o(i+1,j,k) + (2./3.)*(turbvf_x(i+1,j,k)-turbvf_x_o(i+1,j,k)))*momp2x(i+1,j,k) + momp2x(i,j,k)*(turbvf_x_o(i,j,k) + (2./3.)*(turbvf_x(i,j,k)-turbvf_x_o(i,j,k)))
+                cu_fab(i,j,k,4) += 0.5 * (2./3.) * dt * (  (turbvf_x_o(i+1,j,k) + (2./3.)*(turbvf_x(i+1,j,k)-turbvf_x_o(i+1,j,k)))*momp2x(i+1,j,k) + momp2x(i,j,k)*(turbvf_x_o(i,j,k) + (2./3.)*(turbvf_x(i,j,k)-turbvf_x_o(i,j,k)))
                                                + (turbvf_y_o(i,j+1,k) + (2./3.)*(turbvf_y(i,j+1,k)-turbvf_y_o(i,j+1,k)))*momp2y(i,j+1,k) + momp2y(i,j,k)*(turbvf_y_o(i,j,k) + (2./3.)*(turbvf_y(i,j,k)-turbvf_y_o(i,j,k)))
                                                + (turbvf_z_o(i,j,k+1) + (2./3.)*(turbvf_z(i,j,k+1)-turbvf_z_o(i,j,k+1)))*momp2z(i,j,k+1) + momp2z(i,j,k)*(turbvf_z_o(i,j,k) + (2./3.)*(turbvf_z(i,j,k)-turbvf_z_o(i,j,k))) );
-                cu_fab(i,j,k,4) -= dt * energyp2_in; // remove the average input energy from random turbulent forcing
+                cu_fab(i,j,k,4) -= (2./3.) * dt * energyp2_in; // remove the average input energy from random turbulent forcing
             }
         });
     }