diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
index d5c38e72ed64736462aa9c966f4967976ab8f4eb..d305d1a2614b9b8cb740c6ab146ddfd707410918 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
@@ -2,11 +2,6 @@ volScalarField rAUrel(1.0/UrelEqn().A());
 volVectorField HbyA("HbyA", Urel);
 HbyA = rAUrel*UrelEqn().H();
 
-if (pimple.nCorrPISO() <= 1)
-{
-    UrelEqn.clear();
-}
-
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
@@ -16,13 +11,28 @@ surfaceScalarField phiHbyA
 
 adjustPhi(phiHbyA, Urel, p);
 
+tmp<volScalarField> rAtUrel(rAUrel);
+
+if (pimple.consistent())
+{
+    rAtUrel = 1.0/max(1.0/rAUrel - UrelEqn().H1(), 0.1/rAUrel);
+    phiHbyA +=
+        fvc::interpolate(rAtUrel() - rAUrel)*fvc::snGrad(p)*mesh.magSf();
+    HbyA -= (rAUrel - rAtUrel())*fvc::grad(p);
+}
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UrelEqn.clear();
+}
+
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
 {
     // Pressure corrector
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(rAUrel, p) == fvc::div(phiHbyA)
+        fvm::laplacian(rAtUrel(), p) == fvc::div(phiHbyA)
     );
 
     pEqn.setReference(pRefCell, pRefValue);
@@ -40,6 +50,6 @@ while (pimple.correctNonOrthogonal())
 p.relax();
 
 // Momentum corrector
-Urel = HbyA - rAUrel*fvc::grad(p);
+Urel = HbyA - rAtUrel()*fvc::grad(p);
 Urel.correctBoundaryConditions();
 fvOptions.correct(Urel);
diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H
index 0abd31d143ecb99548981885a9134f95864d6a58..13118296c5920f4eac34ca0a4bbda42b8b4821bf 100644
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@@ -1,24 +1,34 @@
-surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
 
-if (pimple.nCorrPISO() <= 1)
-{
-    UEqn.clear();
-}
-
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (fvc::interpolate(HbyA) & mesh.Sf())
-  + rAUf*fvc::ddtCorr(U, phi)
+  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
 );
 
 MRF.makeRelative(phiHbyA);
 
 adjustPhi(phiHbyA, U, p);
 
+tmp<volScalarField> rAtU(rAU);
+
+if (pimple.consistent())
+{
+    rAtU = 1.0/max(1.0/rAU - UEqn().H1(), 0.1/rAU);
+    phiHbyA +=
+        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
+    HbyA -= (rAU - rAtU())*fvc::grad(p);
+}
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UEqn.clear();
+}
+
+surfaceScalarField rAUf("rAUf", fvc::interpolate(rAtU()));
+
 // Update the fixedFluxPressure BCs to ensure flux consistency
 setSnGrad<fixedFluxPressureFvPatchScalarField>
 (
@@ -53,6 +63,6 @@ while (pimple.correctNonOrthogonal())
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U = HbyA - rAU*fvc::grad(p);
+U = HbyA - rAtU()*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index c133a95727e69d6a733c03abe41e42cc434efced..5ad9b80020493835da52896afd949e5237b28cc3 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
@@ -1,18 +1,11 @@
-surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
 
-if (pimple.nCorrPISO() <= 1)
-{
-    UEqn.clear();
-}
-
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (fvc::interpolate(HbyA) & mesh.Sf())
-  + rAUf*fvc::ddtCorr(U, Uf)
+  + fvc::interpolate(rAU)*fvc::ddtCorr(U, Uf)
 );
 
 MRF.makeRelative(phiHbyA);
@@ -24,6 +17,23 @@ if (p.needReference())
     fvc::makeAbsolute(phiHbyA, U);
 }
 
+tmp<volScalarField> rAtU(rAU);
+
+if (pimple.consistent())
+{
+    rAtU = 1.0/max(1.0/rAU - UEqn().H1(), 0.1/rAU);
+    phiHbyA +=
+        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
+    HbyA -= (rAU - rAtU())*fvc::grad(p);
+}
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UEqn.clear();
+}
+
+surfaceScalarField rAUf("rAUf", fvc::interpolate(rAtU()));
+
 // Update the fixedFluxPressure BCs to ensure flux consistency
 setSnGrad<fixedFluxPressureFvPatchScalarField>
 (
@@ -57,7 +67,7 @@ while (pimple.correctNonOrthogonal())
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U = HbyA - rAU*fvc::grad(p);
+U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);