diff --git a/sd1d.cxx b/sd1d.cxx
index eb0fd70..2afbebf 100644
--- a/sd1d.cxx
+++ b/sd1d.cxx
@@ -382,6 +382,8 @@ class SD1D : public PhysicsModel {
 
     mesh->communicate(Ne, NVi, P);
 
+    Field3D P_solver = P; // Save the pressure in the solver
+
     // Floor small values
     P = floor(P, 1e-10);
     Ne = floor(Ne, 1e-10);
@@ -390,7 +392,7 @@ class SD1D : public PhysicsModel {
 
     Vi = NVi / Ne;
 
-    Field3D Te = 0.5 * P / Ne; // Assuming Te = Ti
+    Field3D Te = 0.5 * P / Nelim; // Assuming Te = Ti
 
     for (auto &i : Te.getRegion("RGN_NOBNDRY")) {
       if (Te[i] > 10.)
@@ -1157,7 +1159,7 @@ class SD1D : public PhysicsModel {
       if (rhs_explicit) {
         // Flux splitting with upwinding
         Field3D a = sqrt(gamma_sound * 2. * Te); // Local sound speed
-        ddt(NVi) = -FV::Div_par(NVi, Vi, a, bndry_flux_fix) // Momentum flow
+        ddt(NVi) = -FV::Div_par(NVi, Vi, a, false) // Momentum flow
                    - Grad_par(P);
 
         if (atomic) {
@@ -1271,6 +1273,10 @@ class SD1D : public PhysicsModel {
         if (ADpar > 0.0) {
           ddt(P) += ADpar * AddedDissipation(1.0, P, P, true);
         }
+
+        // Force P towards value in solver.
+        // Note: This only impacts solutions where Ne < floor
+        ddt(P) += 2. * Te * Ne - P_solver;
       }
     }