diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 3a0ca9d042..ab89fcd811 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -111,8 +111,6 @@ else
     thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
 }
 
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
-
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index 546dfaf251..55be311749 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -112,7 +112,6 @@ if (pressureControl.limit(p))
     thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
-
 else if (pimple.SIMPLErho())
 {
     thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
@@ -123,8 +122,6 @@ else
     thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
 }
 
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
-
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index 0f53379eef..677fabffd5 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -72,9 +72,6 @@ while (pimple.correctNonOrthogonal())
 
 p = p_rgh + rho*gh;
 
-// Thermodynamic density update
-//thermo.correctRho(psi*p - psip0);
-
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -95,9 +92,8 @@ if (p_rgh.needReference())
             /compressibility;
         thermo.correctRho(psi*p - psip0);
         rho = thermo.rho();
-
+        p_rgh = p - rho*gh;
     }
-    p_rgh = p - rho*gh;
 }
 else
 {
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index aa3ad0f163..13ffff11b6 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -76,8 +76,8 @@
         {
             p += (initialMass - fvc::domainIntegrate(psi*p))
                 /fvc::domainIntegrate(psi);
+            p_rgh = p - rho*gh;
         }
-        p_rgh = p - rho*gh;
     }
 
     rho = thermo.rho();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index 5b8286ca4b..410fc7166d 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -84,8 +84,8 @@
         {
             p += (initialMass - fvc::domainIntegrate(psi*p))
                 /compressibility;
+            p_rgh = p - rho*gh;
         }
-        p_rgh = p - rho*gh;
     }
 
     rho = thermo.rho();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 0c9195ee9d..dc67ab09e2 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -112,8 +112,8 @@ if (closedVolume)
             /compressibility;
         thermo.correctRho(psi*p - psip0);
         rho = thermo.rho();
+        p_rgh = p - rho*gh;
     }
-    p_rgh = p - rho*gh;
 }
 else
 {
diff --git a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
index bb0364dc79..79a313b208 100644
--- a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
+++ b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
@@ -105,6 +105,7 @@ public:
 
             //- Add the given density correction to the density field.
             //  Used to update the density field following pressure solution
+            //  Limit thermo rho between rhoMin and rhoMax
             virtual void correctRho
             (
                 const volScalarField& deltaRho,
@@ -112,6 +113,14 @@ public:
                 const dimensionedScalar& rhoMax
             ) = 0;
 
+
+            //- 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;
+
             //- 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 b3e65677d9..3dafd068a1 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.C
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
@@ -110,6 +110,12 @@ void Foam::psiThermo::correctRho
 )
 {}
 
+void Foam::psiThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho
+)
+{}
+
 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 b2e05a9e9e..5b24276e0e 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.H
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
@@ -125,6 +125,13 @@ public:
                 const dimensionedScalar& rhoMax
             );
 
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho
+            );
+
             //- 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 bbc1a64ba5..e7b70a4fc9 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
@@ -185,6 +185,12 @@ void Foam::rhoThermo::correctRho
     rho_ = min(rho_, rhoMax);
 }
 
+void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)
+{
+    rho_ += deltaRho;
+}
+
+
 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 fd8bb6b81c..2728cf02e7 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
@@ -146,6 +146,10 @@ public:
                 const dimensionedScalar& rhoMax
             );
 
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            virtual void correctRho(const volScalarField& deltaRho);
+
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const;