diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
index 3eda1b871e377fc01bd5b9245a82d0c9339b2bcf..6922869aa2d20667d2b2967869d129cdf44eda8f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -85,6 +85,10 @@ while (pimple.correct())
 {
     // Update continuity errors due to temperature changes
     fluid.correct();
+    volScalarField rho("rho", fluid.rho());
+
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
 
     // Correct fixed-flux BCs to be consistent with the velocity BCs
     MRF.correctBoundaryFlux(U1, phi1);
@@ -116,9 +120,6 @@ while (pimple.correct())
            *rho2*U2.oldTime()/runTime.deltaT()
         );
 
-    // Mean density for buoyancy force and p_rgh -> p
-    volScalarField rho("rho", fluid.rho());
-
     surfaceScalarField ghSnGradRho
     (
         "ghSnGradRho",
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index 01c120d481f513aed2fb58ac82cbd7179339ef84..ddc262ba788014843b82547f64e416347565b6c4 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -95,6 +95,10 @@ while (pimple.correct())
 {
     // Update continuity errors due to temperature changes
     fluid.correct();
+    volScalarField rho("rho", fluid.rho());
+
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
 
     surfaceScalarField rhof1(fvc::interpolate(rho1));
     surfaceScalarField rhof2(fvc::interpolate(rho2));
@@ -115,7 +119,6 @@ while (pimple.correct())
         max(alphaf2, phase2.residualAlpha())*rAUf2
     );
 
-    volScalarField rho("rho", fluid.rho());
     surfaceScalarField ghSnGradRho
     (
         "ghSnGradRho",
@@ -389,12 +392,6 @@ while (pimple.correct())
         }
     }
 
-    Info<< "min(p) = " << min(p_rgh + rho*gh).value() << endl;
-    if (min(p_rgh + rho*gh) < pMin)
-    {
-        Info<< "Clipping p" << endl;
-    }
-
     // Update and limit the static pressure
     p = max(p_rgh + rho*gh, pMin);