diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index 3ec620dc22..8079428225 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -73,3 +73,26 @@ Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
 #include "createMRF.H"
+
+
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        pimple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        pimple.dict(),
+        dimDensity,
+        0
+    )
+);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 83eeb8fc68..3a0ca9d042 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -84,9 +84,6 @@ else
     }
 }
 
-// Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
-
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -101,12 +98,20 @@ K = 0.5*magSqr(U);
 if (pressureControl.limit(p))
 {
     p.correctBoundaryConditions();
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
 else if (pimple.SIMPLErho())
 {
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
+else
+{
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
+}
+
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
 
 if (thermo.dpdt())
 {
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index e9b849bc05..546dfaf251 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -95,9 +95,6 @@ else
     }
 }
 
-// Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
-
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -112,12 +109,21 @@ K = 0.5*magSqr(U);
 if (pressureControl.limit(p))
 {
     p.correctBoundaryConditions();
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
+
 else if (pimple.SIMPLErho())
 {
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
+else
+{
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
+}
+
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
 
 if (thermo.dpdt())
 {
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index c097118753..0f53379eef 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -1,3 +1,7 @@
+
+dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+bool compressible = (compressibility.value() > SMALL);
+
 rho = thermo.rho();
 
 // Thermodynamic density needs to be updated by psi*d(p) after the
@@ -41,6 +45,12 @@ while (pimple.correctNonOrthogonal())
       - fvm::laplacian(rhorAUf, p_rgh)
     );
 
+    p_rghEqn.setReference
+    (
+        pRefCell,
+        compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+    );
+
     p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
@@ -63,12 +73,40 @@ while (pimple.correctNonOrthogonal())
 p = p_rgh + rho*gh;
 
 // Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
+//thermo.correctRho(psi*p - psip0);
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+if (p_rgh.needReference())
+{
+    if (!compressible)
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+    else
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
+            /compressibility;
+        thermo.correctRho(psi*p - psip0);
+        rho = thermo.rho();
+
+    }
+    p_rgh = p - rho*gh;
+}
+else
+{
+    thermo.correctRho(psi*p - psip0);
+}
+
+rho = thermo.rho();
 
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
 }
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index b9f275407f..64bc8558c8 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -88,3 +88,26 @@ dimensionedScalar totalVolume = sum(mesh.V());
 
 #include "createMRF.H"
 #include "createRadiationModel.H"
+
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        simple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        simple.dict(),
+        dimDensity,
+        0
+    )
+);
+
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index 647cf954bb..aa3ad0f163 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -1,7 +1,4 @@
 {
-    rho = thermo.rho();
-    rho.relax();
-
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
     volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
@@ -24,6 +21,8 @@
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);
 
     while (simple.correctNonOrthogonal())
     {
@@ -32,7 +31,12 @@
             fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );
+
         p_rghEqn.solve();
 
         if (simple.finalNonOrthogonalIter())
@@ -51,12 +55,9 @@
         }
     }
 
-    #include "continuityErrs.H"
-
     p = p_rgh + rho*gh;
 
-    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
-    bool compressible = (compressibility.value() > SMALL);
+    #include "continuityErrs.H"
 
     // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index 7fe544d6a5..5b8286ca4b 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -1,8 +1,10 @@
 {
+    /*
     rho = thermo.rho();
     rho = max(rho, rhoMin[i]);
     rho = min(rho, rhoMax[i]);
     rho.relax();
+    */
 
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -67,10 +69,22 @@
 
     // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity
-    if (closedVolume && compressible)
+    if (closedVolume)
     {
-        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
-            /compressibility;
+        if (!compressible)
+        {
+            p += dimensionedScalar
+            (
+                "p",
+                p.dimensions(),
+                pRefValue - getRefCellValue(p, pRefCell)
+            );
+        }
+        else
+        {
+            p += (initialMass - fvc::domainIntegrate(psi*p))
+                /compressibility;
+        }
         p_rgh = p - rho*gh;
     }
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index be3ec79a2c..87e372fcae 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -19,6 +19,9 @@ List<bool> frozenFlowFluid(fluidRegions.size(), false);
 PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
 PtrList<fv::options> fluidFvOptions(fluidRegions.size());
 
+List<label> refCellFluid(fluidRegions.size());
+List<scalar> refValueFluid(fluidRegions.size());
+
 // Populate fluid field pointer lists
 forAll(fluidRegions, i)
 {
@@ -248,4 +251,20 @@ forAll(fluidRegions, i)
     );
 
     turbulence[i].validate();
+
+    refCellFluid[i] = 0;
+    refValueFluid[i] = 0.0;
+
+    if (p_rghFluid[i].needReference())
+    {
+        setRefCell
+        (
+            thermoFluid[i].p(),
+            p_rghFluid[i],
+            pimpleDict,
+            refCellFluid[i],
+            refValueFluid[i]
+        );
+    }
+
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index d515c6a183..0c9195ee9d 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -4,6 +4,10 @@ bool compressible = (compressibility.value() > SMALL);
 
 rho = thermo.rho();
 
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
 volScalarField rAU("rAU", 1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
@@ -32,10 +36,6 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
       + fvc::div(phiHbyA)
     );
 
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution
-    const volScalarField psip0(psi*p);
-
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix p_rghEqn
@@ -44,6 +44,12 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );
+
         p_rghEqn.solve
         (
             mesh.solver
@@ -62,6 +68,9 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
         if (nonOrth == nNonOrthCorr)
         {
             phi = phiHbyA + p_rghEqn.flux();
+
+            p_rgh.relax();
+
             U = HbyA
               + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
@@ -73,15 +82,10 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
     p = p_rgh + rho*gh;
 
     // Thermodynamic density update
-    thermo.correctRho(psi*p - psip0);
+    //thermo.correctRho(psi*p - psip0);
 }
 
 
-// Update pressure time derivative if needed
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
-}
 
 // Solve continuity
 #include "rhoEqn.H"
@@ -91,10 +95,35 @@ if (thermo.dpdt())
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
-if (closedVolume && compressible)
+if (closedVolume)
 {
-    p += (initialMass - fvc::domainIntegrate(thermo.rho()))
-        /compressibility;
-    rho = thermo.rho();
+    if (!compressible)
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+    else
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
+            /compressibility;
+        thermo.correctRho(psi*p - psip0);
+        rho = thermo.rho();
+    }
     p_rgh = p - rho*gh;
 }
+else
+{
+    thermo.correctRho(psi*p - psip0);
+}
+
+rho = thermo.rho();
+
+// Update pressure time derivative if needed
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 8d3a291bef..5edd27ac5e 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -33,3 +33,6 @@
     );
 
     const bool frozenFlow = frozenFlowFluid[i];
+
+    const label pRefCell = refCellFluid[i];
+    const scalar pRefValue = refValueFluid[i];
diff --git a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
index 995b1e9ad5..bb0364dc79 100644
--- a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
+++ b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
@@ -105,7 +105,12 @@ public:
 
             //- Add the given density correction to the density field.
             //  Used to update the density field following pressure solution
-            virtual void correctRho(const volScalarField& deltaRho) = 0;
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            ) = 0;
 
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const = 0;
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermo.C b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
index fff1f8dfd6..b3e65677d9 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.C
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
@@ -102,10 +102,14 @@ Foam::tmp<Foam::scalarField> Foam::psiThermo::rho(const label patchi) const
 }
 
 
-void Foam::psiThermo::correctRho(const Foam::volScalarField& deltaRho)
+void Foam::psiThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho,
+    const dimensionedScalar& rhoMin,
+    const dimensionedScalar& rhoMax
+)
 {}
 
-
 const Foam::volScalarField& Foam::psiThermo::psi() const
 {
     return psi_;
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermo.H b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
index 454eb5674f..b2e05a9e9e 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.H
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
@@ -118,7 +118,12 @@ public:
             //- Add the given density correction to the density field.
             //  Used to update the density field following pressure solution.
             //  For psiThermo does nothing.
-            virtual void correctRho(const volScalarField& deltaRho);
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            );
 
             //- Density [kg/m^3] - uses current value of pressure
             virtual tmp<volScalarField> rho() const;
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
index 2f2c5d893f..bbc1a64ba5 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
@@ -173,12 +173,18 @@ Foam::volScalarField& Foam::rhoThermo::rho()
 }
 
 
-void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)
+void Foam::rhoThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho,
+    const dimensionedScalar& rhoMin,
+    const dimensionedScalar& rhoMax
+)
 {
     rho_ += deltaRho;
+    rho_ = max(rho_, rhoMin);
+    rho_ = min(rho_, rhoMax);
 }
 
-
 const Foam::volScalarField& Foam::rhoThermo::psi() const
 {
     return psi_;
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
index a4abe62a53..fd8bb6b81c 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
@@ -138,7 +138,13 @@ public:
 
             //- Add the given density correction to the density field.
             //  Used to update the density field following pressure solution
-            virtual void correctRho(const volScalarField& deltaRho);
+            //  Limit thermo rho between rhoMin and rhoMax
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            );
 
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const;