diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index f2f2c5c71c..8e2c3d836d 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -85,8 +85,6 @@ else
     }
 }
 
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
 
 // Explicitly relax pressure for momentum corrector
 p.relax();
@@ -102,6 +100,10 @@ if (pressureControl.limit(p))
 }
 
 thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
 rho = thermo.rho();
 
 // Correct rhoUf if the mesh is moving
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index 3e43dc3d9c..4b22369be9 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -44,6 +44,7 @@ Description
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index 82f63c93fb..22dc129a87 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@@ -35,6 +35,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
+pressureControl pressureControl(p, rho, pimple.dict(), false);
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
@@ -104,3 +105,7 @@ dimensionedScalar initialMass("initialMass", fvc::domainIntegrate(rho));
 #include "createMRF.H"
 #include "createRadiationModel.H"
 #include "createFvOptions.H"
+
+
+const dimensionedScalar rhoMax("rhoMax", dimDensity, GREAT, pimple.dict());
+const dimensionedScalar rhoMin("rhoMin", dimDensity, Zero, pimple.dict());
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index 9a2553a7e9..c36e096dbd 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -72,8 +72,7 @@ while (pimple.correctNonOrthogonal())
 
 p = p_rgh + rho*gh;
 
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
+pressureControl.limit(p);
 
 if (p_rgh.needReference())
 {
@@ -90,16 +89,20 @@ if (p_rgh.needReference())
     {
         p += (initialMass - fvc::domainIntegrate(psi*p))
             /compressibility;
-        thermo.correctRho(psi*p - psip0);
+        thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
         rho = thermo.rho();
         p_rgh = p - rho*gh;
+        p_rgh.correctBoundaryConditions();
     }
 }
 else
 {
-    thermo.correctRho(psi*p - psip0);
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
 }
 
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
 rho = thermo.rho();
 
 if (thermo.dpdt())
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
index df60bcf171..f1b67512be 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -51,6 +51,7 @@ Description
 #include "fvOptions.H"
 #include "coordinateSystem.H"
 #include "loopControl.H"
+#include "pressureControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index d0092f800d..1c9db4f887 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -25,6 +25,11 @@ PtrList<fv::options> fluidFvOptions(fluidRegions.size());
 List<label> pRefCellFluid(fluidRegions.size());
 List<scalar> pRefValueFluid(fluidRegions.size());
 
+PtrList<dimensionedScalar> rhoMinFluid(fluidRegions.size());
+PtrList<dimensionedScalar> rhoMaxFluid(fluidRegions.size());
+
+PtrList<pressureControl> pressureControls(fluidRegions.size());
+
 const uniformDimensionedVectorField& g = meshObjects::gravity::New(runTime);
 
 // Populate fluid field pointer lists
@@ -256,6 +261,24 @@ forAll(fluidRegions, i)
         fluidRegions[i].solutionDict().subDict("PIMPLE");
     pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
 
+    rhoMaxFluid.set
+    (
+        i,
+        new dimensionedScalar("rhoMax", dimDensity, GREAT, pimpleDict)
+    );
+
+    rhoMinFluid.set
+    (
+        i,
+        new dimensionedScalar("rhoMin", dimDensity, Zero, pimpleDict)
+    );
+
+    pressureControls.set
+    (
+        i,
+        new pressureControl(thermoFluid[i].p(), rhoFluid[i], pimpleDict, false)
+    );
+
     Info<< "    Adding MRF\n" << endl;
     MRFfluid.set
     (
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 9bc978b254..206c32eeb2 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -81,17 +81,9 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
     p = p_rgh + rho*gh;
 
-    // Thermodynamic density update
-    //thermo.correctRho(psi*p - psip0);
 }
 
-
-
-// Solve continuity
-#include "rhoEqn.H"
-
-// Update continuity errors
-#include "compressibleContinuityErrors.H"
+pressureControl.limit(p);
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
@@ -110,16 +102,20 @@ if (closedVolume)
     {
         p += (initialMass - fvc::domainIntegrate(psi*p))
             /compressibility;
-        thermo.correctRho(psi*p - psip0);
+        thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
         rho = thermo.rho();
         p_rgh = p - rho*gh;
+        p_rgh.correctBoundaryConditions();
     }
 }
 else
 {
-    thermo.correctRho(psi*p - psip0);
+    thermo.correctRho(psi*p - psip0,  rhoMin, rhoMax);
 }
 
+#include "rhoEqn.H"
+#include "compressibleContinuityErrors.H"
+
 rho = thermo.rho();
 
 // Update pressure time derivative if needed
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 42cfe36e8e..50935d3490 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -61,3 +61,8 @@
     const label pRefCell = pRefCellFluid[i];
     const scalar pRefValue = pRefValueFluid[i];
 
+    const dimensionedScalar rhoMax = rhoMaxFluid[i];
+    const dimensionedScalar rhoMin = rhoMinFluid[i];
+
+    const pressureControl& pressureControl = pressureControls[i];
+