diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index 65cc137acb30ea5c325ae1df63e70971cd2b7c95..521628fa10bdbf1a7f11ead4c9a1fdc4181331f9 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -146,6 +146,10 @@ int main(int argc, char *argv[])
 
         rho = alpha1*rho1 + alpha2*rho2;
 
+        // Correct p_rgh for consistency with p and the updated densities
+        p_rgh = p - rho*gh;
+        p_rgh.correctBoundaryConditions();
+
         runTime.write();
 
         Info<< "ExecutionTime = "
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index 859fc9cc4e8aa240f7d280467dfa42017c67449a..fef73aa36593d93ce1270cf1a1ce5a9b9cb95297 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -116,6 +116,10 @@
 
     rho = alpha1*rho1 + alpha2*rho2;
 
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
+    p_rgh.correctBoundaryConditions();
+
     K = 0.5*magSqr(U);
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index a4a3fb2e2914150398b4ea3f6da31ee5e58e8b4e..78e46f466424fbc38f54ac607cb69e016e18e8a3 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -101,14 +101,16 @@
         }
     }
 
-    // p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-
     // Update densities from change in p_rgh
     rho1 += psi1*(p_rgh - p_rgh_0);
     rho2 += psi2*(p_rgh - p_rgh_0);
 
     rho = alpha1*rho1 + alpha2*rho2;
 
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
+    p_rgh.correctBoundaryConditions();
+
     K = 0.5*magSqr(U);
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
index 70ac37c1b80e29174dfa5c7447107873a05df60e..6526dacb589ce5f76bfaef87782a95ad4d4f45eb 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
@@ -95,9 +95,6 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            // p = max(p_rgh + mixture.rho()*gh, pMin);
-            // p_rgh = p - mixture.rho()*gh;
-
             phasei = 0;
             forAllIter
             (
@@ -125,6 +122,10 @@
     mixture.correctRho(p_rgh - p_rgh_0);
     rho = mixture.rho();
 
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
+    p_rgh.correctBoundaryConditions();
+
     K = 0.5*magSqr(U);
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index 3431a1db24e457c072ea63ed7adc874f37ef93ef..7e43e677492b10b5c07b9deb98fbe4732edef485 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -92,6 +92,11 @@ 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;
+
     PtrList<volVectorField> HbyAs(phases.size());
 
     forAll(phases, phasei)
@@ -121,9 +126,6 @@ while (pimple.correct())
             );
     }
 
-    // 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/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
index 6922869aa2d20667d2b2967869d129cdf44eda8f..9020cd6d5de4cc7c3a1f3bf14b4034bd9cf02fcc 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -85,6 +85,7 @@ 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
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index ddc262ba788014843b82547f64e416347565b6c4..b27c9f4646fdec6beb7d7c731def2287e866d43c 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -95,6 +95,7 @@ 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
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
index a64aeea388201177627be3b4924f0fafab355af5..2010e65ebc41843e14f9e08fdb103fc36debd3bc 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
@@ -83,6 +83,11 @@ while (pimple.correct())
     // Update continuity errors due to temperature changes
     #include "correctContErrs.H"
 
+    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);
     MRF.correctBoundaryFlux(U2, phi2);
@@ -113,9 +118,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/twoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
index bb78b2799d43a3e750514d137ee02fa05433d203..f8388b2d459ade8b784af305049d90876754cf4b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
@@ -95,6 +95,11 @@ while (pimple.correct())
     // Update continuity errors due to temperature changes
     #include "correctContErrs.H"
 
+    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));
 
@@ -114,7 +119,6 @@ while (pimple.correct())
         max(alphaf2, phase2.residualAlpha())*rAUf2
     );
 
-    volScalarField rho("rho", fluid.rho());
     surfaceScalarField ghSnGradRho
     (
         "ghSnGradRho",