diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index 690b53760d75b22618fc9cd0c1358c015d867854..122280cfac48ecd0264c6463fa711928868d7c52 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -44,9 +44,14 @@
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
 
-    dimensionedScalar pMin
+    dimensionedScalar rhoMax
     (
-        mesh.solutionDict().subDict("SIMPLE").lookup("pMin")
+        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMax")
+    );
+
+    dimensionedScalar rhoMin
+    (
+        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMin")
     );
 
     Info<< "Creating turbulence model\n" << endl;
diff --git a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
index e299d99f83c9b1e0c0ab95b35d015300f4f18d7b..57395e977ed3520730cb0efaea428e2bae2b5da9 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
@@ -5,8 +5,8 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho), p, "div(U,p)")
-      - p*fvc::div(phi/fvc::interpolate(rho))
+        fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
+      - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
     );
 
     hEqn.relax();
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 38922c99c9c5700518251a347c167b6848f3de62..329ebbd27fab99c901e1f727669ed2ce21ed5ebb 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -1,4 +1,7 @@
 rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
 
 volScalarField rUA = 1.0/UEqn().A();
 U = rUA*UEqn().H();
@@ -82,15 +85,9 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-rho = thermo.rho();
-rho.relax();
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
-
 U -= rUA*fvc::grad(p);
 U.correctBoundaryConditions();
 
-bound(p, pMin);
-
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
 if (closedVolume)
@@ -98,3 +95,9 @@ if (closedVolume)
     p += (initialMass - fvc::domainIntegrate(psi*p))
         /fvc::domainIntegrate(psi);
 }
+
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;