diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
index d7a6b97ac1..360687a60d 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
@@ -45,14 +45,9 @@
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
 
-    dimensionedScalar rhoMax
+    dimensionedScalar pMin
     (
-        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMax")
-    );
-
-    dimensionedScalar rhoMin
-    (
-        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMin")
+        mesh.solutionDict().subDict("SIMPLE").lookup("pMin")
     );
 
     Info<< "Creating turbulence model\n" << endl;
@@ -70,16 +65,29 @@
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
     thermalPorousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);
 
     // nUCorrectors used for pressureImplicitPorosity
     int nUCorr = 0;
-    const bool pressureImplicitPorosity =
-    (
-        pZones.size()
-     && mesh.solutionDict().subDict("SIMPLE").readIfPresent
-        (
-            "nUCorrectors",
-            nUCorr
-        )
-     && (nUCorr > 0)
-    );
+    if (pZones.size())
+    {
+        // nUCorrectors for pressureImplicitPorosity
+        if (mesh.solutionDict().subDict("SIMPLE").found("nUCorrectors"))
+        {
+            nUCorr = readInt
+            (
+                mesh.solutionDict().subDict("SIMPLE").lookup("nUCorrectors")
+            );
+        }
+
+        if (nUCorr > 0)
+        {
+            pressureImplicitPorosity = true;
+            Info<< "Using pressure implicit porosity" << endl;
+        }
+        else
+        {
+            Info<< "Using pressure explicit porosity" << endl;
+        }
+    }
+
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
index ff0b91bcb7..9e5138de8a 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
@@ -5,8 +5,8 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
-      - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
+        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p, "div(U,p)"))
+      - p*fvc::div(phi/fvc::interpolate(rho))
     );
 
     pZones.addEnthalpySource(thermo, rho, hEqn);
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
index 0386899612..40b7cb3969 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 if (pressureImplicitPorosity)
 {
     U = trTU()&UEqn().H();
@@ -63,6 +58,8 @@ else
 
 U.correctBoundaryConditions();
 
+bound(p, pMin);
+
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
 if (closedVolume)
@@ -72,7 +69,5 @@ if (closedVolume)
 }
 
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
 rho.relax();
 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;