diff --git a/src/hydro/hydro_cuda.cu b/src/hydro/hydro_cuda.cu
index 5651efbfb..1f1d07521 100644
--- a/src/hydro/hydro_cuda.cu
+++ b/src/hydro/hydro_cuda.cu
@@ -385,6 +385,7 @@ __global__ void Update_Conserved_Variables_3D(Real *dev_conserved, Real *Q_Lx, R
       printf("%3d %3d %3d Thread crashed in final update. %e %e %e %e %e\n", xid + x_off, yid + y_off, zid + z_off,
              dev_conserved[id], dtodx * (dev_F_x[imo] - dev_F_x[id]), dtody * (dev_F_y[jmo] - dev_F_y[id]),
              dtodz * (dev_F_z[kmo] - dev_F_z[id]), dev_conserved[4 * n_cells + id]);
+      Average_Cell_All_Fields(xid, yid, zid, nx, ny, nz, n_cells, n_fields, gamma, dev_conserved);
     }
   #endif  // DENSITY_FLOOR
     /*
@@ -653,7 +654,7 @@ __global__ void Average_Slow_Cells_3D(Real *dev_conserved, int nx, int ny, int n
           xid, yid, zid, 1. / max_dti, 1. / max_dti_slow, dev_conserved[id] * DENSITY_UNIT / 0.6 / MP, temp,
           speed * VELOCITY_UNIT * 1e-5, vx * VELOCITY_UNIT * 1e-5, vy * VELOCITY_UNIT * 1e-5, vz * VELOCITY_UNIT * 1e-5,
           cs);
-      Average_Cell_All_Fields(xid, yid, zid, nx, ny, nz, n_cells, n_fields, dev_conserved);
+      Average_Cell_All_Fields(xid, yid, zid, nx, ny, nz, n_cells, n_fields, gamma, dev_conserved);
     }
   }
 }
@@ -1153,23 +1154,94 @@ __device__ Real Average_Cell_Single_Field(int field_indx, int i, int j, int k, i
 }
 
 __device__ void Average_Cell_All_Fields(int i, int j, int k, int nx, int ny, int nz, int ncells, int n_fields,
-                                        Real *conserved)
+                                        Real gamma, Real *conserved)
 {
-  // Average Density
-  Average_Cell_Single_Field(0, i, j, k, nx, ny, nz, ncells, conserved);
-  // Average Momentum_x
-  Average_Cell_Single_Field(1, i, j, k, nx, ny, nz, ncells, conserved);
-  // Average Momentum_y
-  Average_Cell_Single_Field(2, i, j, k, nx, ny, nz, ncells, conserved);
-  // Average Momentum_z
-  Average_Cell_Single_Field(3, i, j, k, nx, ny, nz, ncells, conserved);
-  // Average Energy
-  Average_Cell_Single_Field(4, i, j, k, nx, ny, nz, ncells, conserved);
+  int id = i + (j)*nx + (k)*nx * ny;
+
+  Real d, mx, my, mz, E, P;
+  d  = conserved[grid_enum::density * ncells + id];
+  mx = conserved[grid_enum::momentum_x * ncells + id];
+  my = conserved[grid_enum::momentum_y * ncells + id];
+  mz = conserved[grid_enum::momentum_z * ncells + id];
+  E  = conserved[grid_enum::Energy * ncells + id];
+  P  = (E - (0.5 / d) * (mx * mx + my * my + mz * mz)) * (gamma - 1.0);
+
+  printf("%3d %3d %3d BC: d: %e  E:%e  P:%e  vx:%e  vy:%e  vz:%e\n", i, j, k, d, E, P, mx / d, my / d, mz / d);
+
+  int idn;
+  int N = 0;
+  Real d_av, vx_av, vy_av, vz_av, P_av;
+  d_av = vx_av = vy_av = vz_av = P_av = 0.0;
+  #ifdef SCALAR
+  Real scalar[NSCALARS], scalar_av[NSCALARS];
+  for (int n = 0; n < NSCALARS; n++) {  // NOLINT
+    scalar_av[n] = 0.0;
+  }
+  #endif
+
+  for (int kk = k - 1; kk <= k + 1; kk++) {
+    for (int jj = j - 1; jj <= j + 1; jj++) {
+      for (int ii = i - 1; ii <= i + 1; ii++) {
+        idn = ii + jj * nx + kk * nx * ny;
+        d   = conserved[grid_enum::density * ncells + idn];
+        mx  = conserved[grid_enum::momentum_x * ncells + idn];
+        my  = conserved[grid_enum::momentum_y * ncells + idn];
+        mz  = conserved[grid_enum::momentum_z * ncells + idn];
+        P   = (conserved[grid_enum::Energy * ncells + idn] - (0.5 / d) * (mx * mx + my * my + mz * mz)) * (gamma - 1.0);
+  #ifdef SCALAR
+        for (int n = 0; n < NSCALARS; n++) {  // NOLINT
+          scalar[n] = conserved[grid_enum::scalar * ncells + idn];
+        }
+  #endif
+        if (d > 0.0 && P > 0.0) {
+          d_av += d;
+          vx_av += mx;
+          vy_av += my;
+          vz_av += mz;
+          P_av += P / (gamma - 1.0);
+  #ifdef SCALAR
+          for (int n = 0; n < NSCALARS; n++) {  // NOLINT
+            scalar_av[n] += scalar[n];
+          }
+  #endif
+          N++;
+        }
+      }
+    }
+  }
 
+  P_av  = P_av / N;
+  vx_av = vx_av / d_av;
+  vy_av = vy_av / d_av;
+  vz_av = vz_av / d_av;
+  #ifdef SCALAR
+  for (int n = 0; n < NSCALARS; n++) {  // NOLINT
+    scalar_av[n] = scalar_av[n] / d_av;
+  }
+  #endif
+  d_av = d_av / N;
+
+  // replace cell values with new averaged values
+  conserved[id + ncells * grid_enum::density]    = d_av;
+  conserved[id + ncells * grid_enum::momentum_x] = d_av * vx_av;
+  conserved[id + ncells * grid_enum::momentum_y] = d_av * vy_av;
+  conserved[id + ncells * grid_enum::momentum_z] = d_av * vz_av;
+  conserved[id + ncells * grid_enum::Energy] =
+      P_av / (gamma - 1.0) + 0.5 * d_av * (vx_av * vx_av + vy_av * vy_av + vz_av * vz_av);
   #ifdef DE
-  // Average GasEnergy
-  Average_Cell_Single_Field(n_fields - 1, i, j, k, nx, ny, nz, ncells, conserved);
-  #endif  // DE
+  conserved[id + ncells * grid_enum::GasEnergy] = P_av / (gamma - 1.0);
+  #endif
+  #ifdef SCALAR
+  for (int n = 0; n < NSCALARS; n++) {  // NOLINT
+    conserved[id + ncells * grid_enum::scalar] = d_av * scalar_av[n];
+  }
+  #endif
+
+  d = d_av;
+  E = P_av / (gamma - 1.0) + 0.5 * d_av * (vx_av * vx_av + vy_av * vy_av + vz_av * vz_av);
+  P = P_av;
+
+  printf("%3d %3d %3d FC: d: %e  E:%e  P:%e  vx:%e  vy:%e  vz:%e\n", i, j, k, d, E, P, vx_av, vy_av, vz_av);
 }
 
 #endif  // CUDA
diff --git a/src/hydro/hydro_cuda.h b/src/hydro/hydro_cuda.h
index 016e3a84f..2de29eea3 100644
--- a/src/hydro/hydro_cuda.h
+++ b/src/hydro/hydro_cuda.h
@@ -109,7 +109,7 @@ __global__ void Select_Internal_Energy_2D(Real *dev_conserved, int nx, int ny, i
 __global__ void Select_Internal_Energy_3D(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields);
 
 __device__ void Average_Cell_All_Fields(int i, int j, int k, int nx, int ny, int nz, int ncells, int n_fields,
-                                        Real *conserved);
+                                        Real gamma, Real *conserved);
 
 __device__ Real Average_Cell_Single_Field(int field_indx, int i, int j, int k, int nx, int ny, int nz, int ncells,
                                           Real *conserved);