diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H
index 6fcf83624c621c96b621d78bb06f8597ff86d165..bc093581930e6d3a6ba3109a8e95d8709eb2cdef 100644
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/createFields.H
@@ -52,27 +52,7 @@ volScalarField& p = thermo.p();
 
 #include "compressibleCreatePhi.H"
 
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        pimple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        pimple.dict(),
-        dimDensity,
-        0
-    )
-);
+pressureControl pressureControl(p, rho, pimple.dict(), false);
 
 mesh.setFluxRequired(p.name());
 
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index ac7107acf066722b286275569b0a8b2b39af3e05..81ab6424b4c1a7bfac43df5538ff4d089da5dca0 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -1,7 +1,4 @@
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -87,19 +84,17 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
-    << " " << min(rho).value() << endl;
-
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
diff --git a/applications/solvers/combustion/reactingFoam/pcEqn.H b/applications/solvers/combustion/reactingFoam/pcEqn.H
index 713f443fc5dbed1e31266ae52d845518954f36e9..a1564b0f204ce1d212ae69aa039b29a38ef6676f 100644
--- a/applications/solvers/combustion/reactingFoam/pcEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pcEqn.H
@@ -1,7 +1,4 @@
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
 
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
@@ -109,19 +106,13 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-if (thermo.dpdt())
+if (pressureControl.limit(p))
 {
-    dpdt = fvc::ddt(p);
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
 }
 
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-
-if (!pimple.transonic())
+if (thermo.dpdt())
 {
-    rho.relax();
+    dpdt = fvc::ddt(p);
 }
-
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index a06eef279c3e9e0627a221eefaeb370502aa21fc..d3e1b19ad38b37b22bcef3b81fd336e0c9570bc7 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,6 +37,7 @@ Description
 #include "psiCombustionModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -114,6 +115,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
index 926e69f7c1f7933c402e0289b2e7b0f22cc25bb7..04296ddfc2c53d13dd5fa4f5433e91f6bcb7557f 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
@@ -1,76 +1,74 @@
-{
-    rho = thermo.rho();
+rho = thermo.rho();
 
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
-    surfaceScalarField phiHbyA
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
     (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-      + phig
-    );
-
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-    fvScalarMatrix p_rghDDtEqn
+        fvc::flux(rho*HbyA)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+    )
+  + phig
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+fvScalarMatrix p_rghDDtEqn
+(
+    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+  + fvc::div(phiHbyA)
+  ==
+    fvOptions(psi, p_rgh, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
+{
+    fvScalarMatrix p_rghEqn
     (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-      + fvc::div(phiHbyA)
-      ==
-        fvOptions(psi, p_rgh, rho.name())
+        p_rghDDtEqn
+      - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    while (pimple.correctNonOrthogonal())
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+    if (pimple.finalNonOrthogonalIter())
     {
-        fvScalarMatrix p_rghEqn
-        (
-            p_rghDDtEqn
-          - fvm::laplacian(rhorAUf, p_rgh)
-        );
-
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            // Calculate the conservative fluxes
-            phi = phiHbyA + p_rghEqn.flux();
-
-            // Explicitly relax pressure for momentum corrector
-            p_rgh.relax();
-
-            // Correct the momentum source with the pressure gradient flux
-            // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
-            U.correctBoundaryConditions();
-            fvOptions.correct(U);
-            K = 0.5*magSqr(U);
-        }
+        // Calculate the conservative fluxes
+        phi = phiHbyA + p_rghEqn.flux();
+
+        // Explicitly relax pressure for momentum corrector
+        p_rgh.relax();
+
+        // Correct the momentum source with the pressure gradient flux
+        // calculated from the relaxed pressure
+        U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
     }
+}
 
-    p = p_rgh + rho*gh;
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
+p = p_rgh + rho*gh;
 
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
 
-    #include "rhoEqn.H"
-    #include "compressibleContinuityErrs.H"
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
 }
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
index 687dffcc13db38d9b916c07df8cee4ee1b11c928..e56035f24a703420b00d2e539ee89f3a7af56090 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
@@ -52,6 +52,8 @@ volScalarField& p = thermo.p();
 
 #include "compressibleCreatePhi.H"
 
+pressureControl pressureControl(p, rho, pimple.dict(), false);
+
 mesh.setFluxRequired(p.name());
 
 
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
deleted file mode 100644
index eff014cf43d847917101add0d0be78cf34b64628..0000000000000000000000000000000000000000
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
+++ /dev/null
@@ -1,109 +0,0 @@
-{
-    rho = thermo.rho();
-
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
-
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
-    if (pimple.transonic())
-    {
-        surfaceScalarField phiHbyA
-        (
-            "phiHbyA",
-            (
-                fvc::flux(HbyA)
-              + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
-            )
-        );
-
-        MRF.makeRelative(phiHbyA);
-
-        surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
-
-        phiHbyA *= fvc::interpolate(rho);
-
-        fvScalarMatrix pDDtEqn
-        (
-            fvc::ddt(rho) + fvc::div(phiHbyA)
-          + correction(psi*fvm::ddt(p) + fvm::div(phid, p))
-        );
-
-        while (pimple.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                pDDtEqn
-              - fvm::laplacian(rhorAUf, p)
-             ==
-                fvOptions(psi, p, rho.name())
-            );
-
-            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-            if (pimple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
-        }
-    }
-    else
-    {
-        surfaceScalarField phiHbyA
-        (
-            "phiHbyA",
-            (
-                fvc::flux(rho*HbyA)
-              + rhorAUf*fvc::ddtCorr(rho, U, phi)
-            )
-        );
-
-        MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-        // Update the pressure BCs to ensure flux consistency
-        constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
-        fvScalarMatrix pDDtEqn
-        (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-          + fvc::div(phiHbyA)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
-
-        while (pimple.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                pDDtEqn
-              - fvm::laplacian(rhorAUf, p)
-            );
-
-            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-            if (pimple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
-        }
-    }
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
-
-    #include "rhoEqn.H"
-    #include "compressibleContinuityErrs.H"
-
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
-    K = 0.5*magSqr(U);
-
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-}
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C
index 55de274d839e85cbd63f37423b72a967b12bad17..5be897c6c7f4903ea4bf33d561c7b0a1ca48db46 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,6 +38,7 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -100,7 +101,14 @@ int main(int argc, char *argv[])
             // --- Pressure corrector loop
             while (pimple.correct())
             {
-                #include "pEqn.H"
+                if (pimple.consistent())
+                {
+                    #include "../../../compressible/rhoPimpleFoam/pcEqn.H"
+                }
+                else
+                {
+                    #include "../../../compressible/rhoPimpleFoam/pEqn.H"
+                }
             }
 
             if (pimple.turbCorr())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index 84f70f2cd02c8c7c0c8a406ea5324cf0dc480a0b..8079428225be4c497065ec40e6679df97d007252 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -42,6 +42,8 @@ volVectorField U
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
+mesh.setFluxRequired(p.name());
+
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -54,8 +56,6 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-mesh.setFluxRequired(p.name());
-
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
 (
@@ -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 5aab5362390ad6292ee1d9ff2caa4fabf5920ba2..80d5d9b7946945d031643b515ed54a930a953430 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -1,3 +1,12 @@
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -10,7 +19,7 @@ if (pimple.nCorrPISO() <= 1)
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    fvc::flux(rho*HbyA)
+    fvc::interpolate(rho)*fvc::flux(HbyA)
   + rhorAUf*fvc::ddtCorr(rho, U, phi)
 );
 
@@ -26,19 +35,20 @@ if (pimple.transonic())
         "phid",
         (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
-    phiHbyA -= fvc::interpolate(p)*phid;
+
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA) + fvm::div(phid, p)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
 
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
@@ -53,16 +63,17 @@ if (pimple.transonic())
 }
 else
 {
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -84,15 +95,14 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-pressureControl.limit(p);
-p.correctBoundaryConditions();
-rho = thermo.rho();
-
-if (!pimple.transonic())
+if (pressureControl.limit(p))
 {
-    rho.relax();
+    p.correctBoundaryConditions();
 }
 
+thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
+rho = thermo.rho();
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index afbc2851e4a3b5298ff153d25c4a3a0ecd580c2f..55be311749e7c83dbc5b22f48d530e31a147f3c9 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -1,3 +1,12 @@
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -11,7 +20,7 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (
-        fvc::flux(rho*HbyA)
+        fvc::interpolate(rho)*fvc::flux(HbyA)
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
     )
 );
@@ -33,21 +42,21 @@ if (pimple.transonic())
 
     phiHbyA +=
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-      - fvc::interpolate(p)*phid;
+      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA) + fvm::div(phid, p)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAtU, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
 
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
@@ -65,16 +74,17 @@ else
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAtU, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -96,13 +106,20 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-pressureControl.limit(p);
-p.correctBoundaryConditions();
-rho = thermo.rho();
-
-if (!pimple.transonic())
+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
 {
-    rho.relax();
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
 }
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index 0f855626189a25b93e5937bec958adaaeb9e608a..ed802ba6d2cd197ff21bd5248d4ef1b2ab9c5b37 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -1,3 +1,5 @@
+rho = thermo.rho();
+
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -27,13 +29,14 @@ if (pimple.transonic())
         "phid",
         (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
-    phiHbyA -= fvc::interpolate(p)*phid;
+
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, p)
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phiHbyA)
           + fvm::div(phid, p)
           - fvm::laplacian(rhorAUf, p)
@@ -59,7 +62,7 @@ else
         // Pressure corrector
         fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, p)
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phiHbyA)
           - fvm::laplacian(rhorAUf, p)
          ==
@@ -86,13 +89,10 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-pressureControl.limit(p);
-p.correctBoundaryConditions();
-rho = thermo.rho();
-
-if (!pimple.transonic())
+if (pressureControl.limit(p))
 {
-    rho.relax();
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
 }
 
 {
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index 0fbe7717e3569bec9026d98bd30d3ad51dc6a1ba..5d5ed3f5f4c3961b3e2e941e4a42278d73038bb5 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -122,6 +122,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index c7edbb9e877a550265257e4d0c2139e3a1b70e46..0f290ed273478ffc8756e6b01d395adfe3b5f1ec 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -1,109 +1,110 @@
-{
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-    tUEqn.clear();
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+tUEqn.clear();
+
+bool closedVolume = false;
 
-    bool closedVolume = false;
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+if (simple.transonic())
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+    );
 
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
-    if (simple.transonic())
+    while (simple.correctNonOrthogonal())
     {
-        surfaceScalarField phid
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rhorAUf, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
+        pEqn.setReference
         (
-            "phid",
-            (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+            pressureControl.refCell(),
+            pressureControl.refValue()
         );
-        phiHbyA -= fvc::interpolate(p)*phid;
 
-        while (simple.correctNonOrthogonal())
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
         {
-            fvScalarMatrix pEqn
-            (
-                fvc::div(phiHbyA)
-              + fvm::div(phid, p)
-              - fvm::laplacian(rhorAUf, p)
-              ==
-                fvOptions(psi, p, rho.name())
-            );
-
-            // Relax the pressure equation to ensure diagonal-dominance
-            pEqn.relax();
-
-            pEqn.setReference
-            (
-                pressureControl.refCell(),
-                pressureControl.refValue()
-            );
-
-            pEqn.solve();
-
-            if (simple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
+            phi = phiHbyA + pEqn.flux();
         }
     }
-    else
+}
+else
+{
+    closedVolume = adjustPhi(phiHbyA, U, p);
+
+    while (simple.correctNonOrthogonal())
     {
-        closedVolume = adjustPhi(phiHbyA, U, p);
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAUf, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
 
-        while (simple.correctNonOrthogonal())
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
         {
-            fvScalarMatrix pEqn
-            (
-                fvc::div(phiHbyA)
-              - fvm::laplacian(rhorAUf, p)
-              ==
-                fvOptions(psi, p, rho.name())
-            );
-
-            pEqn.setReference
-            (
-                pressureControl.refCell(),
-                pressureControl.refValue()
-            );
-
-            pEqn.solve();
-
-            if (simple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
+            phi = phiHbyA + pEqn.flux();
         }
     }
+}
 
+#include "incompressible/continuityErrs.H"
 
-    #include "incompressible/continuityErrs.H"
-
-    // Explicitly relax pressure for momentum corrector
-    p.relax();
+// Explicitly relax pressure for momentum corrector
+p.relax();
 
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
 
-    pressureControl.limit(p);
+bool pLimited = pressureControl.limit(p);
 
-    // For closed-volume cases adjust the pressure and density levels
-    // to obey overall mass continuity
-    if (closedVolume)
-    {
-        p += (initialMass - fvc::domainIntegrate(psi*p))
-            /fvc::domainIntegrate(psi);
-    }
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+if (closedVolume)
+{
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
+}
 
+if (pLimited || closedVolume)
+{
     p.correctBoundaryConditions();
+}
 
-    rho = thermo.rho();
+rho = thermo.rho();
 
-    if (!simple.transonic())
-    {
-        rho.relax();
-    }
+if (!simple.transonic())
+{
+    rho.relax();
 }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
index c7dc0f864d2997439cbba04ae2ef5a60053e0a98..51560aab5fd9c6162fcc4847e8134b16be27424a 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
@@ -1,3 +1,5 @@
+rho = thermo.rho();
+
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -5,7 +7,7 @@ tUEqn.clear();
 
 bool closedVolume = false;
 
-surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
 volScalarField rhorAtU("rhorAtU", rho*rAtU);
@@ -23,7 +25,7 @@ if (simple.transonic())
 
     phiHbyA +=
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-      - fvc::interpolate(p)*phid;
+      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
@@ -98,7 +100,7 @@ U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 
-pressureControl.limit(p);
+bool pLimited = pressureControl.limit(p);
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
@@ -108,9 +110,11 @@ if (closedVolume)
         /fvc::domainIntegrate(psi);
 }
 
-p.correctBoundaryConditions();
+if (pLimited || closedVolume)
+{
+    p.correctBoundaryConditions();
+}
 
-// Recalculate density from the relaxed pressure
 rho = thermo.rho();
 
 if (!simple.transonic())
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index 15c87f9c85eecf51d1a6d83e0633646ca1b14758..677fabffd5614bef92210a9a51b9585bbcbfacc1 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -1,115 +1,108 @@
-{
-    bool closedVolume = p_rgh.needReference();
 
-    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
-    bool compressible = (compressibility.value() > SMALL);
+dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+bool compressible = (compressibility.value() > SMALL);
 
-    rho = thermo.rho();
+rho = thermo.rho();
 
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p_rgh;
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
 
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
-    surfaceScalarField phiHbyA
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
     (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-      + phig
+        fvc::flux(rho*HbyA)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+    )
+  + phig
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+fvScalarMatrix p_rghDDtEqn
+(
+    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+  + fvc::div(phiHbyA)
+  ==
+    fvOptions(psi, p_rgh, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
+{
+    fvScalarMatrix p_rghEqn
+    (
+        p_rghDDtEqn
+      - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-    tmp<fvScalarMatrix> p_rghDDtEqn
+    p_rghEqn.setReference
     (
-        new fvScalarMatrix(p_rgh, dimMass/dimTime)
+        pRefCell,
+        compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
     );
 
-    if (compressible)
-    {
-        p_rghDDtEqn =
-        (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-        );
-    }
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-    while (pimple.correctNonOrthogonal())
+    if (pimple.finalNonOrthogonalIter())
     {
-        fvScalarMatrix p_rghEqn
-        (
-            p_rghDDtEqn()
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p_rgh)
-        ==
-            fvOptions(psi, p_rgh, rho.name())
-        );
-
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            // Calculate the conservative fluxes
-            phi = phiHbyA + p_rghEqn.flux();
-
-            // Explicitly relax pressure for momentum corrector
-            p_rgh.relax();
-
-            // Correct the momentum source with the pressure gradient flux
-            // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
-            U.correctBoundaryConditions();
-            fvOptions.correct(U);
-            K = 0.5*magSqr(U);
-        }
+        // Calculate the conservative fluxes
+        phi = phiHbyA + p_rghEqn.flux();
+
+        // Explicitly relax pressure for momentum corrector
+        p_rgh.relax();
+
+        // Correct the momentum source with the pressure gradient flux
+        // calculated from the relaxed pressure
+        U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
     }
+}
 
-    p = p_rgh + rho*gh;
+p = p_rgh + rho*gh;
 
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p_rgh;
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
 
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-
-    if (compressible)
+if (p_rgh.needReference())
+{
+    if (!compressible)
     {
-        #include "rhoEqn.H"
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
     }
-    #include "compressibleContinuityErrs.H"
-
-    if (closedVolume)
+    else
     {
-        if (!compressible)
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pRefCell)
-            );
-        }
-        else
-        {
-            p += (initialMass - fvc::domainIntegrate(thermo.rho()))
-                /compressibility;
-            rho = thermo.rho();
-        }
+        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();
 
-    Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
-        << endl;
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
 }
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index b9f275407f4ceb335c60de4158636dec4ded797c..64bc8558c8f2f1ddaee4d75bd6d1fed04273a22e 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 647cf954bb2d5168adc098f00764a368c5af1608..13ffff11b649b7b44889d2909c01ab9ff9a4e9d3 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
@@ -75,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 7fe544d6a57afe68850f81d39483f79c02f10660..410fc7166d14f31febb352aa18b0d841ea9e65ef 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,11 +69,23 @@
 
     // 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;
-        p_rgh = p - rho*gh;
+        if (!compressible)
+        {
+            p += dimensionedScalar
+            (
+                "p",
+                p.dimensions(),
+                pRefValue - getRefCellValue(p, pRefCell)
+            );
+        }
+        else
+        {
+            p += (initialMass - fvc::domainIntegrate(psi*p))
+                /compressibility;
+            p_rgh = p - rho*gh;
+        }
     }
 
     rho = thermo.rho();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index be3ec79a2cce235f969e379915a95a645f7bc9ce..87e372fcae3b5deef8b5df76ba8de942ade4ee93 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 09d90906c0237c79dc632492b6e3675db3e23f15..dc67ab09e20d88a3f05242417bed3905a4a4c439 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -1,114 +1,129 @@
-{
-    bool closedVolume = p_rgh.needReference();
-    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
-    bool compressible = (compressibility.value() > SMALL);
+bool closedVolume = p_rgh.needReference();
+dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+bool compressible = (compressibility.value() > SMALL);
+
+rho = thermo.rho();
 
-    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));
+volScalarField rAU("rAU", 1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
 
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
-    surfaceScalarField phiHbyA
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
     (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-      + phig
-    );
+        fvc::flux(rho*HbyA)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+    )
+  + phig
+);
 
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
-    tmp<fvScalarMatrix> p_rghDDtEqn
+{
+    fvScalarMatrix p_rghDDtEqn
     (
-        new fvScalarMatrix(p_rgh, dimMass/dimTime)
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+      + fvc::div(phiHbyA)
     );
 
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        if (compressible)
-        {
-            p_rghDDtEqn =
-            (
-                fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-             ==
-                fvOptions(psi, p_rgh, rho.name())
-            );
-        }
-
-        // Thermodynamic density needs to be updated by psi*d(p) after the
-        // pressure solution - done in 2 parts. Part 1:
-        thermo.rho() -= psi*p_rgh;
+        fvScalarMatrix p_rghEqn
+        (
+            p_rghDDtEqn
+          - fvm::laplacian(rhorAUf, p_rgh)
+        );
 
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-        {
-            fvScalarMatrix p_rghEqn
-            (
-                p_rghDDtEqn()
-              + fvc::div(phiHbyA)
-              - fvm::laplacian(rhorAUf, p_rgh)
-            );
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );
 
-            p_rghEqn.solve
+        p_rghEqn.solve
+        (
+            mesh.solver
             (
-                mesh.solver
+                p_rgh.select
                 (
-                    p_rgh.select
                     (
-                        (
-                           oCorr == nOuterCorr-1
-                        && corr == nCorr-1
-                        && nonOrth == nNonOrthCorr
-                        )
+                        oCorr == nOuterCorr-1
+                     && corr == nCorr-1
+                     && nonOrth == nNonOrthCorr
                     )
                 )
-            );
-
-            if (nonOrth == nNonOrthCorr)
-            {
-                phi = phiHbyA + p_rghEqn.flux();
-                U = HbyA
-                  + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
-                U.correctBoundaryConditions();
-                fvOptions.correct(U);
-                K = 0.5*magSqr(U);
-            }
-        }
+            )
+        );
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi = phiHbyA + p_rghEqn.flux();
+
+            p_rgh.relax();
 
-        // Second part of thermodynamic density update
-        thermo.rho() += psi*p_rgh;
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+            K = 0.5*magSqr(U);
+        }
     }
 
     p = p_rgh + rho*gh;
 
-    // Update pressure time derivative if needed
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
+    // Thermodynamic density update
+    //thermo.correctRho(psi*p - psip0);
+}
+
 
-    if (compressible)
-    {
-        // Solve continuity
-        #include "rhoEqn.H"
-    }
 
-    // Update continuity errors
-    #include "compressibleContinuityErrors.H"
+// Solve continuity
+#include "rhoEqn.H"
 
-    // For closed-volume cases adjust the pressure and density levels
-    // to obey overall mass continuity
-    if (closedVolume && compressible)
+// Update continuity errors
+#include "compressibleContinuityErrors.H"
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+if (closedVolume)
+{
+    if (!compressible)
     {
-        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+        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 8d3a291bef9cd98d110bd2ee8ffadb4b20854300..5edd27ac5e506079d003e29b1979d5606c81eb3d 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/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
index 544cccdae2ee420c8991c9894cc9ad6b5d1e492c..c7ea95929c950877d6c3f1ce512aa3659988672a 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,6 +42,7 @@ Description
 #include "radiationModel.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index f6c8ceec21d28dc6064fee4c129d1d875d724b02..f2cb2947186adfcee7a90c18cbdcc63c2e8ff7fb 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -1,73 +1,71 @@
-{
-    rho = thermo.rho();
+rho = thermo.rho();
 
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-    surfaceScalarField phiHbyA
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
     (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-    );
+        fvc::flux(rho*HbyA)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+    )
+);
 
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
 
-    fvScalarMatrix pDDtEqn
+fvScalarMatrix pDDtEqn
+(
+    fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+  + fvc::div(phiHbyA)
+ ==
+    parcels.Srho()
+  + fvOptions(psi, p, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
+{
+    fvScalarMatrix pEqn
     (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phiHbyA)
-     ==
-        parcels.Srho()
-      + fvOptions(psi, p, rho.name())
+        pDDtEqn
+      - fvm::laplacian(rhorAUf, p)
     );
 
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            pDDtEqn
-          - fvm::laplacian(rhorAUf, p)
-        );
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA + pEqn.flux();
-        }
+    if (pimple.finalNonOrthogonalIter())
+    {
+        phi = phiHbyA + pEqn.flux();
     }
+}
 
-    p.relax();
+p.relax();
 
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
 
-    #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
-    #include "compressibleContinuityErrs.H"
+#include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
+#include "compressibleContinuityErrs.H"
 
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
-    K = 0.5*magSqr(U);
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
 
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
 
-    rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
 
-    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
-}
+Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
index 055eff6f0b34ec7e30fb3da5d8ac9f6f96386d25..1ce5db0ec8f1fbd4758d126be362ebb0dd2d3354 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
@@ -1,57 +1,55 @@
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+tUEqn.clear();
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    fvc::interpolate(rho)*fvc::flux(HbyA)
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+while (simple.correctNonOrthogonal())
 {
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
-
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-    tUEqn.clear();
-    surfaceScalarField phiHbyA
+    fvScalarMatrix pEqn
     (
-        "phiHbyA",
-        fvc::interpolate(rho)*fvc::flux(HbyA)
+        fvc::div(phiHbyA)
+      - fvm::laplacian(rhorAUf, p)
+     ==
+        parcels.Srho()
+      + fvOptions(psi, p, rho.name())
     );
 
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+    pEqn.solve();
 
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
-    while (simple.correctNonOrthogonal())
+    if (simple.finalNonOrthogonalIter())
     {
-        fvScalarMatrix pEqn
-        (
-            fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            parcels.Srho()
-          + fvOptions(psi, p, rho.name())
-        );
-
-        pEqn.solve();
-
-        if (simple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA + pEqn.flux();
-        }
+        phi = phiHbyA + pEqn.flux();
     }
+}
 
-    p.relax();
+p.relax();
 
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
 
-    #include "compressibleContinuityErrs.H"
+#include "compressibleContinuityErrs.H"
 
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
 
-    rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
-    rho.relax();
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
 
-    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
-}
+Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index 7cd3875e85c6ed1c72b5e155ce687d3634cec71c..e84c19453ffcdfbb87c71099f8b536096a55fba3 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -104,8 +104,8 @@
     }
 
     // Update densities from change in p_rgh
-    rho1 += psi1*(p_rgh - p_rgh_0);
-    rho2 += psi2*(p_rgh - p_rgh_0);
+    mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
+    mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
 
     rho = alpha1*rho1 + alpha2*rho2;
 
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 46168e6b87f7bfb44bba6c1313f37bf2e72b1787..0616b80ea817e02d219348a6db3de57f73ac41d8 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -38,8 +38,8 @@ volScalarField& alpha2(mixture.alpha2());
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-volScalarField& rho1 = mixture.thermo1().rho();
-volScalarField& rho2 = mixture.thermo2().rho();
+const volScalarField& rho1 = mixture.thermo1().rho();
+const volScalarField& rho2 = mixture.thermo2().rho();
 
 volScalarField rho
 (
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 51c58ba53689d0b3582c9c0407c81d0136b90f01..93412c9dbaca55742f28d77bd15f74e354a15a93 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -104,8 +104,8 @@
     }
 
     // Update densities from change in p_rgh
-    rho1 += psi1*(p_rgh - p_rgh_0);
-    rho2 += psi2*(p_rgh - p_rgh_0);
+    mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
+    mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
 
     rho = alpha1*rho1 + alpha2*rho2;
 
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.C
index 5aee70ba33eb27e57bbce9b2625aa9bc88ca9e8c..ccd9cc3b5efa7b487f3d89fd0b7ff22099a8ff19 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.C
@@ -45,7 +45,7 @@ fixedIncidentRadiationFvPatchScalarField
 :
     fixedGradientFvPatchScalarField(p, iF),
     temperatureCoupledBase(patch(), "undefined", "undefined", "undefined-K"),
-    QrIncident_(p.size(), 0.0)
+    qrIncident_(p.size(), 0.0)
 {}
 
 
@@ -60,7 +60,7 @@ fixedIncidentRadiationFvPatchScalarField
 :
     fixedGradientFvPatchScalarField(psf, p, iF, mapper),
     temperatureCoupledBase(patch(), psf),
-    QrIncident_(psf.QrIncident_)
+    qrIncident_(psf.qrIncident_)
 {}
 
 
@@ -74,7 +74,7 @@ fixedIncidentRadiationFvPatchScalarField
 :
     fixedGradientFvPatchScalarField(p, iF),
     temperatureCoupledBase(patch(), dict),
-    QrIncident_("QrIncident", dict, p.size())
+    qrIncident_("qrIncident", dict, p.size())
 {
     if (dict.found("value") && dict.found("gradient"))
     {
@@ -99,7 +99,7 @@ fixedIncidentRadiationFvPatchScalarField
 :
     fixedGradientFvPatchScalarField(psf, iF),
     temperatureCoupledBase(patch(), psf),
-    QrIncident_(psf.QrIncident_)
+    qrIncident_(psf.qrIncident_)
 {}
 
 
@@ -111,7 +111,7 @@ fixedIncidentRadiationFvPatchScalarField
 :
     fixedGradientFvPatchScalarField(ptf),
     temperatureCoupledBase(patch(), ptf),
-    QrIncident_(ptf.QrIncident_)
+    qrIncident_(ptf.qrIncident_)
 {}
 
 
@@ -123,7 +123,7 @@ void Foam::radiation::fixedIncidentRadiationFvPatchScalarField::autoMap
 )
 {
     fixedGradientFvPatchScalarField::autoMap(m);
-    QrIncident_.autoMap(m);
+    qrIncident_.autoMap(m);
 }
 
 
@@ -141,7 +141,7 @@ void Foam::radiation::fixedIncidentRadiationFvPatchScalarField::rmap
             psf
         );
 
-    QrIncident_.rmap(thftpsf.QrIncident_, addr);
+    qrIncident_.rmap(thftpsf.qrIncident_, addr);
 }
 
 
@@ -165,7 +165,7 @@ void Foam::radiation::fixedIncidentRadiationFvPatchScalarField::updateCoeffs()
     gradient() =
         emissivity
        *(
-            QrIncident_
+            qrIncident_
           - physicoChemical::sigma.value()*pow4(*this)
         )/kappa(*this);
 
@@ -173,11 +173,11 @@ void Foam::radiation::fixedIncidentRadiationFvPatchScalarField::updateCoeffs()
 
     if (debug)
     {
-        scalar Qr = gSum(kappa(*this)*gradient()*patch().magSf());
+        scalar qr = gSum(kappa(*this)*gradient()*patch().magSf());
         Info<< patch().boundaryMesh().mesh().name() << ':'
             << patch().name() << ':'
             << this->internalField().name() << " -> "
-            << " radiativeFlux:" << Qr
+            << " radiativeFlux:" << qr
             << " walltemperature "
             << " min:" << gMin(*this)
             << " max:" << gMax(*this)
@@ -194,7 +194,7 @@ void Foam::radiation::fixedIncidentRadiationFvPatchScalarField::write
 {
     fixedGradientFvPatchScalarField::write(os);
     temperatureCoupledBase::write(os);
-    QrIncident_.writeEntry("QrIncident", os);
+    qrIncident_.writeEntry("qrIncident", os);
     writeEntry("value", os);
 }
 
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.H b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.H
index 79b48191fe3d97b582b2e303b4dd278884a65d4a..987adcb89caa359cc1b10afb038c1882e29b010a 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.H
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.H
@@ -34,19 +34,19 @@ Description
 
     the gradient heat flux is calculated as :
 
-        Qr = emissivity*(QrIncident - sigma_*T^4)
+        qr = emissivity*(qrIncident - sigma_*T^4)
 
     where:
 
     emissivity is the emissivity of the solid.
-    QrIncident is the specified fixed incident radiation.
+    qrIncident is the specified fixed incident radiation.
 
     Example usage:
 
     wall
     {
         type            fixedIncidentRadiation;
-        QrIncident      uniform 500;
+        qrIncident      uniform 500;
         kappa           solidThermo;
         KappaName       none;
     }
@@ -89,7 +89,7 @@ class fixedIncidentRadiationFvPatchScalarField
     // Private data
 
         //- Incident radiative heat flux
-        scalarField QrIncident_;
+        scalarField qrIncident_;
 
 
 public:
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C
index 6d9d8b171b1dbf6b04accf31b59abb0baca1d447..52de5b816bb2cb3258c9f67745f149714459033e 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C
@@ -150,18 +150,18 @@ void Foam::lumpedMassWallTemperatureFvPatchScalarField::updateCoeffs()
 
     if (debug)
     {
-        scalar qin(0);
-        scalar qout(0);
+        scalar Qin(0);
+        scalar Qout(0);
 
         forAll(q, facei)
         {
-            if (q[facei] > 0.0) //out the wall
+            if (q[facei] > 0.0) // out the wall
             {
-                qout += q[facei]*magSf[facei];
+                Qout += q[facei]*magSf[facei];
             }
-            else if (q[facei] < 0.0) //into the wall
+            else if (q[facei] < 0.0) // into the wall
             {
-                qin += q[facei]*magSf[facei];
+                Qin += q[facei]*magSf[facei];
             }
         }
 
@@ -173,8 +173,8 @@ void Foam::lumpedMassWallTemperatureFvPatchScalarField::updateCoeffs()
             << " min:" << gMin(*this)
             << " max:" << gMax(*this)
             << " avg:" << gAverage(*this)
-            << " Qin [W]:" << qin
-            << " Qout [W]:" << qout
+            << " Qin [W]:" << Qin
+            << " Qout [W]:" << Qout
             << endl;
     }
 
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C
index 6ed6e4f72afb060258182e8c36403d934863d02a..4221ad93de1d4b126d77672c84ba5fc0439b860f 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C
@@ -50,7 +50,7 @@ thermalBaffle1DFvPatchScalarField
     TName_("T"),
     baffleActivated_(true),
     thickness_(p.size()),
-    Qs_(p.size()),
+    qs_(p.size()),
     solidDict_(),
     solidPtr_(nullptr),
     qrPrevious_(p.size()),
@@ -74,7 +74,7 @@ thermalBaffle1DFvPatchScalarField
     TName_(ptf.TName_),
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_, mapper),
-    Qs_(ptf.Qs_, mapper),
+    qs_(ptf.qs_, mapper),
     solidDict_(ptf.solidDict_),
     solidPtr_(ptf.solidPtr_),
     qrPrevious_(ptf.qrPrevious_, mapper),
@@ -97,7 +97,7 @@ thermalBaffle1DFvPatchScalarField
     TName_("T"),
     baffleActivated_(dict.lookupOrDefault<bool>("baffleActivated", true)),
     thickness_(),
-    Qs_(p.size(), 0),
+    qs_(p.size(), 0),
     solidDict_(dict),
     solidPtr_(),
     qrPrevious_(p.size(), 0.0),
@@ -111,9 +111,9 @@ thermalBaffle1DFvPatchScalarField
         thickness_ = scalarField("thickness", dict, p.size());
     }
 
-    if (dict.found("Qs"))
+    if (dict.found("qs"))
     {
-        Qs_ = scalarField("Qs", dict, p.size());
+        qs_ = scalarField("qs", dict, p.size());
     }
 
     if (dict.found("qrPrevious"))
@@ -151,7 +151,7 @@ thermalBaffle1DFvPatchScalarField
     TName_(ptf.TName_),
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_),
-    Qs_(ptf.Qs_),
+    qs_(ptf.qs_),
     solidDict_(ptf.solidDict_),
     solidPtr_(ptf.solidPtr_),
     qrPrevious_(ptf.qrPrevious_),
@@ -173,7 +173,7 @@ thermalBaffle1DFvPatchScalarField
     TName_(ptf.TName_),
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_),
-    Qs_(ptf.Qs_),
+    qs_(ptf.qs_),
     solidDict_(ptf.solidDict_),
     solidPtr_(ptf.solidPtr_),
     qrPrevious_(ptf.qrPrevious_),
@@ -264,11 +264,11 @@ baffleThickness() const
 
 
 template<class solidType>
-tmp<scalarField> thermalBaffle1DFvPatchScalarField<solidType>::Qs() const
+tmp<scalarField> thermalBaffle1DFvPatchScalarField<solidType>::qs() const
 {
     if (this->owner())
     {
-         return Qs_;
+         return qs_;
     }
     else
     {
@@ -283,10 +283,10 @@ tmp<scalarField> thermalBaffle1DFvPatchScalarField<solidType>::Qs() const
             nbrPatch.template lookupPatchField<volScalarField, scalar>(TName_)
         );
 
-        tmp<scalarField> tQs(new scalarField(nbrField.Qs()));
-        scalarField& Qs = tQs.ref();
-        mapDist.distribute(Qs);
-        return tQs;
+        tmp<scalarField> tqs(new scalarField(nbrField.qs()));
+        scalarField& qs = tqs.ref();
+        mapDist.distribute(qs);
+        return tqs;
     }
 }
 
@@ -304,7 +304,7 @@ void thermalBaffle1DFvPatchScalarField<solidType>::autoMap
     if (this->owner())
     {
         thickness_.autoMap(m);
-        Qs_.autoMap(m);
+        qs_.autoMap(m);
     }
 }
 
@@ -324,7 +324,7 @@ void thermalBaffle1DFvPatchScalarField<solidType>::rmap
     if (this->owner())
     {
         thickness_.rmap(tiptf.thickness_, addr);
-        Qs_.rmap(tiptf.Qs_, addr);
+        qs_.rmap(tiptf.qs_, addr);
     }
 }
 
@@ -397,7 +397,7 @@ void thermalBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
 
         valueFraction() = alpha/(alpha + myKDelta);
 
-        refValue() = (KDeltaSolid*nbrTp + Qs()/2.0)/alpha;
+        refValue() = (KDeltaSolid*nbrTp + qs()/2.0)/alpha;
 
         if (debug)
         {
@@ -431,7 +431,7 @@ void thermalBaffle1DFvPatchScalarField<solidType>::write(Ostream& os) const
     if (this->owner())
     {
         baffleThickness()().writeEntry("thickness", os);
-        Qs()().writeEntry("Qs", os);
+        qs()().writeEntry("qs", os);
         solid().write(os);
     }
 
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.H b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.H
index f9491806fecfdde79494d021c10b96b524745e30..d93092d699fbbf3d8eed03d45c465e037b89887f 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.H
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.H
@@ -46,7 +46,7 @@ Usage
         samplePatch     <slavePatchName>;
 
         thickness       uniform 0.005;  // thickness [m]
-        Qs              uniform 100;    // heat flux [W/m2]
+        qs              uniform 100;    // heat flux [W/m2]
 
         qr              none;
         relaxation      1;
@@ -124,7 +124,7 @@ class thermalBaffle1DFvPatchScalarField
         mutable scalarField thickness_;
 
         //- Superficial heat source [W/m2]
-        mutable scalarField Qs_;
+        mutable scalarField qs_;
 
         //- Solid dictionary
         dictionary solidDict_;
@@ -147,8 +147,8 @@ class thermalBaffle1DFvPatchScalarField
         //- Return const solid thermo
         const solidType& solid() const;
 
-        //- Return Qs from master
-        tmp<scalarField> Qs() const;
+        //- Return qs from master
+        tmp<scalarField> qs() const;
 
         //- Return thickness from master
         tmp<scalarField> baffleThickness() const;
diff --git a/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C b/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C
index a672011e28aa1e884baf658a997d95608b2ae917..7a0cfe1a08814054a4c6ba1544e3d6220496ac8c 100644
--- a/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C
+++ b/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C
@@ -219,7 +219,7 @@ bool Foam::slidingInterface::projectPoints() const
             projectionAlgo_
         );
 
-//     Pout<< "USING N-SQAURED!!!" << endl;
+//     Pout<< "USING N-SQUARED!!!" << endl;
 //     List<objectHit> slavePointFaceHits =
 //         projectPointsNSquared<face, List, const pointField&>
 //         (
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
index aca22e769907b6bfaec86482e1cfd11f62f32335..bcd8ebd5f09ca99b719c2c1033715b937f51f9b3 100644
--- a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
@@ -38,24 +38,30 @@ Foam::pressureControl::pressureControl
 :
     refCell_(-1),
     refValue_(0),
-    pMax_("pMax", dimPressure, 0),
-    pMin_("pMin", dimPressure, GREAT)
+    pMax_("pMax", dimPressure, GREAT),
+    pMin_("pMin", dimPressure, 0),
+    limitMaxP_(false),
+    limitMinP_(false)
 {
     bool pLimits = false;
+    scalar pMax = -GREAT;
+    scalar pMin = GREAT;
 
     // Set the reference cell and value for closed domain simulations
     if (pRefRequired && setRefCell(p, dict, refCell_, refValue_))
     {
         pLimits = true;
 
-        pMax_.value() = refValue_;
-        pMin_.value() = refValue_;
+        pMax = refValue_;
+        pMin = refValue_;
     }
 
     if (dict.found("pMax") && dict.found("pMin"))
     {
         pMax_.value() = readScalar(dict.lookup("pMax"));
+        limitMaxP_ = true;
         pMin_.value() = readScalar(dict.lookup("pMin"));
+        limitMinP_ = true;
     }
     else
     {
@@ -73,8 +79,8 @@ Foam::pressureControl::pressureControl
                 pLimits = true;
                 rhoLimits = true;
 
-                pMax_.value() = max(pMax_.value(), max(pbf[patchi]));
-                pMin_.value() = min(pMin_.value(), min(pbf[patchi]));
+                pMax = max(pMax, max(pbf[patchi]));
+                pMin = min(pMin, min(pbf[patchi]));
 
                 rhoRefMax = max(rhoRefMax, max(rhobf[patchi]));
                 rhoRefMin = min(rhoRefMin, min(rhobf[patchi]));
@@ -84,8 +90,8 @@ Foam::pressureControl::pressureControl
         reduce(rhoLimits, andOp<bool>());
         if (rhoLimits)
         {
-            reduce(pMax_.value(), maxOp<scalar>());
-            reduce(pMin_.value(), minOp<scalar>());
+            reduce(pMax, maxOp<scalar>());
+            reduce(pMin, minOp<scalar>());
 
             reduce(rhoRefMax, maxOp<scalar>());
             reduce(rhoRefMin, minOp<scalar>());
@@ -94,6 +100,7 @@ Foam::pressureControl::pressureControl
         if (dict.found("pMax"))
         {
             pMax_.value() = readScalar(dict.lookup("pMax"));
+            limitMaxP_ = true;
         }
         else if (dict.found("pMaxFactor"))
         {
@@ -108,7 +115,8 @@ Foam::pressureControl::pressureControl
             }
 
             const scalar pMaxFactor(readScalar(dict.lookup("pMaxFactor")));
-            pMax_ *= pMaxFactor;
+            pMax_.value() = pMaxFactor*pMax;
+            limitMaxP_ = true;
         }
         else if (dict.found("rhoMax"))
         {
@@ -142,12 +150,14 @@ Foam::pressureControl::pressureControl
 
             dimensionedScalar rhoMax("rhoMax", dimDensity, dict);
 
-            pMax_ *= max(rhoMax.value()/rhoRefMax, 1);
+            pMax_.value() = max(rhoMax.value()/rhoRefMax, 1)*pMax;
+            limitMaxP_ = true;
         }
 
         if (dict.found("pMin"))
         {
             pMin_.value() = readScalar(dict.lookup("pMin"));
+            limitMinP_ = true;
         }
         else if (dict.found("pMinFactor"))
         {
@@ -162,7 +172,8 @@ Foam::pressureControl::pressureControl
             }
 
             const scalar pMinFactor(readScalar(dict.lookup("pMinFactor")));
-            pMin_ *= pMinFactor;
+            pMin_.value() = pMinFactor*pMin;
+            limitMinP_ = true;
         }
         else if (dict.found("rhoMin"))
         {
@@ -195,26 +206,64 @@ Foam::pressureControl::pressureControl
 
             dimensionedScalar rhoMin("rhoMin", dimDensity, dict);
 
-            pMin_ *= min(rhoMin.value()/rhoRefMin, 1);
+            pMin_.value() = min(rhoMin.value()/rhoRefMin, 1)*pMin;
+            limitMinP_ = true;
         }
     }
 
-    Info<< "pressureControl" << nl
-        << "    pMax/pMin " << pMax_.value() << " " << pMin_.value()
-        << nl << endl;
+    if (limitMaxP_ || limitMinP_)
+    {
+        Info<< "pressureControl" << nl;
+
+        if (limitMaxP_)
+        {
+            Info<< "    pMax " << pMax_.value() << nl;
+        }
+
+        if (limitMinP_)
+        {
+            Info<< "    pMin " << pMin_.value() << nl;
+        }
+
+        Info << endl;
+    }
 }
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::pressureControl::limit(volScalarField& p) const
+bool Foam::pressureControl::limit(volScalarField& p) const
 {
-    Info<< "pressureControl: p max/min "
-        << max(p).value() << " "
-        << min(p).value() << endl;
+    if (limitMaxP_ || limitMinP_)
+    {
+        if (limitMaxP_)
+        {
+            const scalar pMax = max(p).value();
 
-    p = max(p, pMin_);
-    p = min(p, pMax_);
+            if (pMax > pMax_.value())
+            {
+                Info<< "pressureControl: p max " << pMax << endl;
+                p = min(p, pMax_);
+            }
+        }
+
+        if (limitMinP_)
+        {
+            const scalar pMin = min(p).value();
+
+            if (pMin < pMin_.value())
+            {
+                Info<< "pressureControl: p min " << pMin << endl;
+                p = max(p, pMin_);
+            }
+        }
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
index 61f510ac60301d489b3ec26522d04e005b6545aa..fb6caa9ff32f972371c1ad6edcf4dd25ee24dc19 100644
--- a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
@@ -66,6 +66,12 @@ class pressureControl
         //- Pressure upper-limit
         dimensionedScalar pMin_;
 
+        //- Pressure lower-limit
+        bool limitMaxP_;
+
+        //- Pressure upper-limit
+        bool limitMinP_;
+
 
 public:
 
@@ -89,8 +95,8 @@ public:
         //- Return the pressure reference level
         inline scalar refValue() const;
 
-        //- Limit the pressure
-        void limit(volScalarField& p) const;
+        //- Limit the pressure if necessary and return true if so
+        bool limit(volScalarField& p) const;
 };
 
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
index e500ab75872fa4676acf2ca686e929150e1b644f..6d031f05b6dacf5e9a87b44158710b6a79c00e90 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,6 +45,7 @@ void Foam::pimpleControl::read()
     solveFlow_ = pimpleDict.lookupOrDefault<Switch>("solveFlow", true);
     nCorrPIMPLE_ = pimpleDict.lookupOrDefault<label>("nOuterCorrectors", 1);
     nCorrPISO_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
+    SIMPLErho_ = pimpleDict.lookupOrDefault<Switch>("SIMPLErho", false);
     turbOnFinalIterOnly_ =
         pimpleDict.lookupOrDefault<Switch>("turbOnFinalIterOnly", true);
 }
@@ -128,6 +129,7 @@ Foam::pimpleControl::pimpleControl(fvMesh& mesh, const word& dictName)
     nCorrPIMPLE_(0),
     nCorrPISO_(0),
     corrPISO_(0),
+    SIMPLErho_(false),
     turbOnFinalIterOnly_(true),
     converged_(false)
 {
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
index 033ce84d2f54a6029712f74f7c49c584fc96402d..5e87b20ea33832a19914a830f794341584ae1a9a 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,6 +81,10 @@ protected:
             //- Current PISO corrector
             label corrPISO_;
 
+            //- Flag to indicate whether to update density in SIMPLE
+            //  rather than PISO mode
+            bool SIMPLErho_;
+
             //- Flag to indicate whether to only solve turbulence on final iter
             bool turbOnFinalIterOnly_;
 
@@ -128,6 +132,10 @@ public:
             //- Current PISO corrector index
             inline label corrPISO() const;
 
+            //- Flag to indicate whether to update density in SIMPLE
+            //  rather than PISO mode
+            inline bool SIMPLErho() const;
+
 
         // Solution control
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
index f6c6e0ce35fbd7bd100455ed0035f46ba28cb34b..deed2958c0e4ba023fd4f626e1dd9717ee198a83 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,12 @@ inline Foam::label Foam::pimpleControl::corrPISO() const
 }
 
 
+inline bool Foam::pimpleControl::SIMPLErho() const
+{
+    return SIMPLErho_;
+}
+
+
 inline bool Foam::pimpleControl::correct()
 {
     corrPISO_++;
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C
index 305743076c27f4c6c8873b2ba4e1bb2018bf98c2..0933f4bda8ea4e27998b90b1fe72d04f4cc749f6 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C
@@ -92,6 +92,7 @@ Foam::fv::correctedSnGrad<Type>::correction
         )
     );
     GeometricField<Type, fvsPatchField, surfaceMesh>& ssf = tssf.ref();
+    ssf.setOriented();
 
     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
diff --git a/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H b/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H
index 8893abb3e57c15324de1d6abae9f8bfb72ac9396..9e276fa01e4c70b04873a7b222307a93f6b7d02c 100644
--- a/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H
+++ b/src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H
@@ -48,7 +48,7 @@ Description
     \endvartable
 
 
-    The distribution inside the hear exchanger is given by:
+    The distribution inside the heat exchanger is given by:
     \f[
         Q_c = \frac{V_c |U_c| (T_c - T_{ref})}{\sum(V_c |U_c| (T_c - T_{ref}))}
     \f]
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C
index ec9236d802e5487659a49a08ca0a2691954cf500..d8d91f045a9954db6a10861fe8d47ad8938d73c3 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C
@@ -56,11 +56,11 @@ primaryRadiation::primaryRadiation
 )
 :
     filmRadiationModel(typeName, film, dict),
-    QinPrimary_
+    qinPrimary_
     (
         IOobject
         (
-            "Qin", // same name as Qin on primary region to enable mapping
+            "qin", // same name as qin on primary region to enable mapping
             film.time().timeName(),
             film.regionMesh(),
             IOobject::NO_READ,
@@ -83,8 +83,8 @@ primaryRadiation::~primaryRadiation()
 
 void primaryRadiation::correct()
 {
-    // Transfer Qin from primary region
-    QinPrimary_.correctBoundaryConditions();
+    // Transfer qin from primary region
+    qinPrimary_.correctBoundaryConditions();
 }
 
 
@@ -108,10 +108,10 @@ tmp<volScalarField> primaryRadiation::Shs()
     );
 
     scalarField& Shs = tShs.ref();
-    const scalarField& QinP = QinPrimary_;
+    const scalarField& qinP = qinPrimary_;
     const scalarField& alpha = filmModel_.alpha();
 
-    Shs = QinP*alpha;
+    Shs = qinP*alpha;
 
     return tShs;
 }
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.H b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.H
index b93561d49e83bc8249c7d847a3701ceae9005727..0334c8de1a6ed0f052e75daa08395b3c7ecfcc2f 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.H
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.H
@@ -61,7 +61,7 @@ private:
     // Private data
 
         //- Incident radiative flux mapped from  the primary region / [kg/s3]
-        volScalarField QinPrimary_;
+        volScalarField qinPrimary_;
 
 
     // Private member functions
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
index f333752847f435abf238e5af01eea0e2c0c97156..d210d5e7b0d5f5f2a839b94a9d2a0302c99e704a 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
@@ -57,11 +57,11 @@ standardRadiation::standardRadiation
 )
 :
     filmRadiationModel(typeName, film, dict),
-    QinPrimary_
+    qinPrimary_
     (
         IOobject
         (
-            "Qin", // same name as Qin on primary region to enable mapping
+            "qin", // same name as qin on primary region to enable mapping
             film.time().timeName(),
             film.regionMesh(),
             IOobject::NO_READ,
@@ -101,7 +101,7 @@ standardRadiation::~standardRadiation()
 void standardRadiation::correct()
 {
     // Transfer qr from primary region
-    QinPrimary_.correctBoundaryConditions();
+    qinPrimary_.correctBoundaryConditions();
 }
 
 
@@ -125,14 +125,14 @@ tmp<volScalarField> standardRadiation::Shs()
     );
 
     scalarField& Shs = tShs.ref();
-    const scalarField& QinP = QinPrimary_;
+    const scalarField& qinP = qinPrimary_;
     const scalarField& delta = filmModel_.delta();
     const scalarField& alpha = filmModel_.alpha();
 
-    Shs = beta_*QinP*alpha*(1.0 - exp(-kappaBar_*delta));
+    Shs = beta_*qinP*alpha*(1.0 - exp(-kappaBar_*delta));
 
     // Update net qr on local region
-    qrNet_.primitiveFieldRef() = QinP - Shs;
+    qrNet_.primitiveFieldRef() = qinP - Shs;
     qrNet_.correctBoundaryConditions();
 
     return tShs;
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.H b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.H
index 7375ef79ec9e3e75bfd8b01ae807f27565c3d97f..3562166202eecb5769889b20c586139d087802ee 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.H
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.H
@@ -60,7 +60,7 @@ private:
     // Private data
 
         //- Radiative incident flux mapped from  the primary region / [kg/s3]
-        volScalarField QinPrimary_;
+        volScalarField qinPrimary_;
 
         //- Remaining radiative flux after removing local contribution
         volScalarField qrNet_;
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H
index e770d8b56cb93d8c76dacc73cec9f9a0ab65ece7..3bb38b436779066bd7f26d6524b3fac550124b67 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H
@@ -367,10 +367,10 @@ public:
         // Derived fields (calculated on-the-fly)
 
             //- Return the convective heat energy from film to wall
-            inline tmp<scalarField> Qconvw(const label patchi) const;
+            inline tmp<scalarField> qconvw(const label patchi) const;
 
             //- Return the convective heat energy from primary region to film
-            inline tmp<scalarField> Qconvp(const label patchi) const;
+            inline tmp<scalarField> qconvp(const label patchi) const;
 
 
         // Evolution
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
index bf9f8e60c76354990b6254da648fefda51be5a1d..38fdd485400b2194299f691f17ac67e07ea09ca5 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
@@ -155,7 +155,7 @@ inline const filmRadiationModel& thermoSingleLayer::radiation() const
 }
 
 
-inline tmp<scalarField> thermoSingleLayer::Qconvw(const label patchi) const
+inline tmp<scalarField> thermoSingleLayer::qconvw(const label patchi) const
 {
     const scalarField htc(htcw_->h()().boundaryField()[patchi]);
     const scalarField& Tp = T_.boundaryField()[patchi];
@@ -165,7 +165,7 @@ inline tmp<scalarField> thermoSingleLayer::Qconvw(const label patchi) const
 }
 
 
-inline tmp<scalarField> thermoSingleLayer::Qconvp(const label patchi) const
+inline tmp<scalarField> thermoSingleLayer::qconvp(const label patchi) const
 {
     const scalarField htc(htcs_->h()().boundaryField()[patchi]);
     const scalarField& Tp = T_.boundaryField()[patchi];
diff --git a/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.C b/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.C
index 5ebe085f2d15247122b043b780267235c619a842..61475876004f4e38a642830b52d16a1cf10b335f 100644
--- a/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.C
+++ b/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.C
@@ -97,7 +97,7 @@ void thermalBaffle::solveEnergy()
     volScalarField alpha("alpha", thermo_->alpha());
 
 
-    //If region is one-dimension variable thickness
+    // If region is one-dimension variable thickness
     if (oneD_ && !constantThickness_)
     {
         // Scale K and rhoCp and fill Q in the internal baffle region.
@@ -112,7 +112,7 @@ void thermalBaffle::solveEnergy()
                 const label cellId = cells[i];
 
                 Q[cellId] =
-                    Qs_.boundaryField()[patchi][localFacei]
+                    qs_.boundaryField()[patchi][localFacei]
                    /thickness_[localFacei];
 
                 rho[cellId] *= delta_.value()/thickness_[localFacei];
@@ -167,11 +167,11 @@ thermalBaffle::thermalBaffle
     nNonOrthCorr_(readLabel(solution().lookup("nNonOrthCorr"))),
     thermo_(solidThermo::New(regionMesh(), dict)),
     h_(thermo_->he()),
-    Qs_
+    qs_
     (
         IOobject
         (
-            "Qs",
+            "qs",
             regionMesh().time().timeName(),
             regionMesh(),
             IOobject::READ_IF_PRESENT,
@@ -227,11 +227,11 @@ thermalBaffle::thermalBaffle
     nNonOrthCorr_(readLabel(solution().lookup("nNonOrthCorr"))),
     thermo_(solidThermo::New(regionMesh())),
     h_(thermo_->he()),
-    Qs_
+    qs_
     (
         IOobject
         (
-            "Qs",
+            "qs",
             regionMesh().time().timeName(),
             regionMesh(),
             IOobject::READ_IF_PRESENT,
@@ -289,13 +289,13 @@ void thermalBaffle::init()
     if (oneD_ && !constantThickness_)
     {
         label patchi = intCoupledPatchIDs_[0];
-        const label Qsb = Qs_.boundaryField()[patchi].size();
+        const label qsb = qs_.boundaryField()[patchi].size();
 
-        if (Qsb!= thickness_.size())
+        if (qsb!= thickness_.size())
         {
             FatalErrorInFunction
-                << "the boundary field of Qs is "
-                << Qsb << " and " << nl
+                << "the boundary field of qs is "
+                << qsb << " and " << nl
                 << "the field 'thickness' is " << thickness_.size() << nl
                 << exit(FatalError);
         }
diff --git a/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.H b/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.H
index 48768dca9313a2398e233f98091df79f04823eaa..9016284b6052a49e8b25bf5ab7d5ddddc504e3de 100644
--- a/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.H
+++ b/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.H
@@ -94,7 +94,7 @@ protected:
         // Source term fields
 
             //- Surface energy source  / [J/m2/s]
-            volScalarField Qs_;
+            volScalarField qs_;
 
             //- Volumetric energy source  / [J/m3/s]
             volScalarField Q_;
diff --git a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
index e18cda73431f6485264520edd12ffb4e872f8acc..79a313b20828ef0cd87662e3979f8b28a50ed919 100644
--- a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
+++ b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,24 @@ public:
 
         // Access to thermodynamic state variables
 
+            //- 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,
+                const dimensionedScalar& rhoMin,
+                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 630402280a448259cd8ca3d233d95d2cdd28222d..3dafd068a13dc564c294d92d155602bff355c2ec 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.C
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,6 +102,20 @@ Foam::tmp<Foam::scalarField> Foam::psiThermo::rho(const label patchi) const
 }
 
 
+void Foam::psiThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho,
+    const dimensionedScalar& rhoMin,
+    const dimensionedScalar& rhoMax
+)
+{}
+
+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 76404a026b36f3668d94da12296b0d0ab93b8c7e..5b24276e0ea735281da5ad4d75a343bf67e05a5c 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.H
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,6 +115,23 @@ public:
 
         // Fields derived from thermodynamic state variables
 
+            //- 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,
+                const dimensionedScalar& rhoMin,
+                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 ada7716c4a8ee0b69c8bd722b2c36e3eec1b5e7d..e7b70a4fc920f0f00b56992dd74618bfeb8d1e3e 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -173,6 +173,24 @@ Foam::volScalarField& Foam::rhoThermo::rho()
 }
 
 
+void Foam::rhoThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho,
+    const dimensionedScalar& rhoMin,
+    const dimensionedScalar& rhoMax
+)
+{
+    rho_ += deltaRho;
+    rho_ = max(rho_, rhoMin);
+    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 680094b90f0a6e28b6903961e8507c20489272a0..2728cf02e74e1ffba593f643128daf6fcfab764a 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -136,6 +136,20 @@ public:
             //- Return non-const access to the local density field [kg/m^3]
             virtual volScalarField& rho();
 
+            //- 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,
+                const dimensionedScalar& rhoMin,
+                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;
 
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
index 72b0803f1c281b6d37ce50e869daf5647fde1ccd..70423f1a722fc6adcb5d02356a1aa19f45c26ba8 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
@@ -184,18 +184,18 @@ updateCoeffs()
 
     const scalarField& emissivity = temissivity();
 
-    scalarField& Qem = ray.Qem().boundaryFieldRef()[patchi];
-    scalarField& Qin = ray.Qin().boundaryFieldRef()[patchi];
+    scalarField& qem = ray.qem().boundaryFieldRef()[patchi];
+    scalarField& qin = ray.qin().boundaryFieldRef()[patchi];
 
     const vector& myRayId = dom.IRay(rayId).d();
 
     // Use updated Ir while iterating over rays
-    // avoids to used lagged Qin
-    scalarField Ir = dom.IRay(0).Qin().boundaryField()[patchi];
+    // avoids to used lagged qin
+    scalarField Ir = dom.IRay(0).qin().boundaryField()[patchi];
 
     for (label rayI=1; rayI < dom.nRay(); rayI++)
     {
-        Ir += dom.IRay(rayI).Qin().boundaryField()[patchi];
+        Ir += dom.IRay(rayI).qin().boundaryField()[patchi];
     }
 
     if (solarLoad_)
@@ -221,7 +221,7 @@ updateCoeffs()
                 )/pi;
 
             // Emmited heat flux from this ray direction
-            Qem[faceI] = refValue()[faceI]*nAve[faceI];
+            qem[faceI] = refValue()[faceI]*nAve[faceI];
         }
         else
         {
@@ -231,7 +231,7 @@ updateCoeffs()
             refValue()[faceI] = 0.0; //not used
 
             // Incident heat flux on this ray direction
-            Qin[faceI] = Iw[faceI]*nAve[faceI];
+            qin[faceI] = Iw[faceI]*nAve[faceI];
         }
     }
 
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
index 102a8f87c5a9218fc625ea4a9805c5fdb4167166..8bf1e455f669ee8ea2b869102cdcfad8286acbb6 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
@@ -177,17 +177,17 @@ updateCoeffs()
 
     const scalarField& emissivity = temissivity();
 
-    scalarField& Qem = ray.Qem().boundaryFieldRef()[patchi];
-    scalarField& Qin = ray.Qin().boundaryFieldRef()[patchi];
+    scalarField& qem = ray.qem().boundaryFieldRef()[patchi];
+    scalarField& qin = ray.qin().boundaryFieldRef()[patchi];
 
     // Use updated Ir while iterating over rays
-    // avoids to used lagged Qin
+    // avoids to used lagged qin
     /*
-    scalarField Ir = dom.IRay(0).Qin().boundaryField()[patchi];
+    scalarField Ir = dom.IRay(0).qin().boundaryField()[patchi];
 
     for (label rayI=1; rayI < dom.nRay(); rayI++)
     {
-        Ir += dom.IRay(rayI).Qin().boundaryField()[patchi];
+        Ir += dom.IRay(rayI).qin().boundaryField()[patchi];
     }
     */
 
@@ -227,7 +227,7 @@ updateCoeffs()
                 )/pi;
 
             // Emmited heat flux from this ray direction (sum over lambdaId)
-            Qem[facei] += refValue()[facei]*nAve[facei];
+            qem[facei] += refValue()[facei]*nAve[facei];
         }
         else
         {
@@ -237,7 +237,7 @@ updateCoeffs()
             refValue()[facei] = 0.0; //not used
 
             // Incident heat flux on this ray direction (sum over lambdaId)
-            Qin[facei] += Iw[facei]*nAve[facei];
+            qin[facei] += Iw[facei]*nAve[facei];
         }
     }
 
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
index 71cca246dc21c930b51910215ef50ecb413bce79..23fd128691b6e4ea775302e4619f269dfe2adaf1 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
@@ -243,31 +243,31 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         mesh_,
         dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    Qem_
+    qem_
     (
         IOobject
         (
-            "Qem",
+            "qem",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qem", dimMass/pow3(dimTime), 0.0)
     ),
-    Qin_
+    qin_
     (
         IOobject
         (
-            "Qin",
+            "qin",
             mesh_.time().timeName(),
             mesh_,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qin", dimMass/pow3(dimTime), 0.0)
     ),
     a_
     (
@@ -336,31 +336,31 @@ Foam::radiation::fvDOM::fvDOM
         mesh_,
         dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    Qem_
+    qem_
     (
         IOobject
         (
-            "Qem",
+            "qem",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qem", dimMass/pow3(dimTime), 0.0)
     ),
-    Qin_
+    qin_
     (
         IOobject
         (
-            "Qin",
+            "qin",
             mesh_.time().timeName(),
             mesh_,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qin", dimMass/pow3(dimTime), 0.0)
     ),
     a_
     (
@@ -509,16 +509,16 @@ void Foam::radiation::fvDOM::updateG()
 {
     G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
     qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
-    Qem_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
-    Qin_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
+    qem_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
+    qin_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
 
     forAll(IRay_, rayI)
     {
         IRay_[rayI].addIntensity();
         G_ += IRay_[rayI].I()*IRay_[rayI].omega();
         qr_.boundaryFieldRef() += IRay_[rayI].qr().boundaryField();
-        Qem_.boundaryFieldRef() += IRay_[rayI].Qem().boundaryField();
-        Qin_.boundaryFieldRef() += IRay_[rayI].Qin().boundaryField();
+        qem_.boundaryFieldRef() += IRay_[rayI].qem().boundaryField();
+        qin_.boundaryFieldRef() += IRay_[rayI].qin().boundaryField();
     }
 }
 
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H
index 84d12a057e8a587cfc441be7f4eb738109f71fbf..3c00d3b2872808eb14e16c46ebc35d2ae2b26863 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H
@@ -107,10 +107,10 @@ class fvDOM
         volScalarField qr_;
 
          //- Emmited radiative heat flux [W/m2]
-        volScalarField Qem_;
+        volScalarField qem_;
 
         //- Incidet radiative heat flux [W/m2]
-        volScalarField Qin_;
+        volScalarField qin_;
 
         //- Total absorption coefficient [1/m]
         volScalarField a_;
@@ -255,10 +255,10 @@ public:
             inline const volScalarField& qr() const;
 
             //- Const access to incident radiative heat flux field
-            inline const volScalarField& Qin() const;
+            inline const volScalarField& qin() const;
 
             //- Const access to emitted radiative heat flux field
-            inline const volScalarField& Qem() const;
+            inline const volScalarField& qem() const;
 
             //- Const access to black body
             inline const blackBodyEmission& blackBody() const;
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOMI.H b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOMI.H
index dc64c7b76c218cf53883dff3d9910f3fc7eeafb0..64e3aa7027b3241fbb468438bcf235f388543341 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOMI.H
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOMI.H
@@ -91,15 +91,15 @@ inline const Foam::volScalarField& Foam::radiation::fvDOM::qr() const
     return qr_;
 }
 
-inline const Foam::volScalarField& Foam::radiation::fvDOM::Qin() const
+inline const Foam::volScalarField& Foam::radiation::fvDOM::qin() const
 {
-    return Qin_;
+    return qin_;
 }
 
 
-inline const Foam::volScalarField& Foam::radiation::fvDOM::Qem() const
+inline const Foam::volScalarField& Foam::radiation::fvDOM::qem() const
 {
-    return Qem_;
+    return qem_;
 }
 
 
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
index adf5f54912b675c5c089e340c038562a3fd118bb..84a70d45e897da379063d732efd99ac2e868d488 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
@@ -80,31 +80,31 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
         mesh_,
         dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    Qin_
+    qin_
     (
         IOobject
         (
-            "Qin" + name(rayId),
+            "qin" + name(rayId),
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qin", dimMass/pow3(dimTime), 0.0)
     ),
-    Qem_
+    qem_
     (
         IOobject
         (
-            "Qem" + name(rayId),
+            "qem" + name(rayId),
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qem", dimMass/pow3(dimTime), 0.0)
     ),
     d_(Zero),
     dAve_(Zero),
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
index 2dc4eb411346fb521432c4941994f196f2cffa9f..7cd2713b662f165db53979b6ac9b72c0042a5c04 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
@@ -83,10 +83,10 @@ private:
         volScalarField qr_;
 
         //- Incident radiative heat flux on boundary
-        volScalarField Qin_;
+        volScalarField qin_;
 
         //- Emitted radiative heat flux on boundary
-        volScalarField Qem_;
+        volScalarField qem_;
 
         //- Direction
         vector d_;
@@ -182,16 +182,16 @@ public:
             inline volScalarField& qr();
 
             //- Return non-const access to the boundary incident heat flux
-            inline volScalarField& Qin();
+            inline volScalarField& qin();
 
             //- Return non-const access to the boundary emmited heat flux
-            inline volScalarField& Qem();
+            inline volScalarField& qem();
 
             //- Return const access to the boundary incident heat flux
-            inline const volScalarField& Qin() const;
+            inline const volScalarField& qin() const;
 
             //- Return const access to the boundary emmited heat flux
-            inline const volScalarField& Qem() const;
+            inline const volScalarField& qem() const;
 
             //- Return direction
             inline const vector& d() const;
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRayI.H b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRayI.H
index d307d890bc0233154148a7c670fc60a4921dffc9..0356aa782ec331e2b9773590235d4b03118aa4c4 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRayI.H
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRayI.H
@@ -43,28 +43,28 @@ inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::qr()
 }
 
 inline const Foam::volScalarField& Foam::radiation::
-radiativeIntensityRay::Qin() const
+radiativeIntensityRay::qin() const
 {
-    return Qin_;
+    return qin_;
 }
 
 
-inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::Qin()
+inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::qin()
 {
-    return Qin_;
+    return qin_;
 }
 
 
 inline const Foam::volScalarField& Foam::radiation::
-radiativeIntensityRay::Qem() const
+radiativeIntensityRay::qem() const
 {
-    return Qem_;
+    return qem_;
 }
 
 
-inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::Qem()
+inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::qem()
 {
-    return Qem_;
+    return qem_;
 }
 
 
diff --git a/src/thermophysicalModels/radiation/radiationModels/radiationModel/radiationModel.C b/src/thermophysicalModels/radiation/radiationModels/radiationModel/radiationModel.C
index 1f097e7855781935ba2c9a305b5fd9f64ec5304b..b8de0f96bcd4219702abd3fcf0043ccbaf7d881b 100644
--- a/src/thermophysicalModels/radiation/radiationModels/radiationModel/radiationModel.C
+++ b/src/thermophysicalModels/radiation/radiationModels/radiationModel/radiationModel.C
@@ -43,7 +43,7 @@ namespace Foam
 }
 
 const Foam::word Foam::radiation::radiationModel::externalRadHeatFieldName_ =
-    "QrExt";
+    "qrExt";
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
index 866865fbcb2ffe263beb1806f05c6388c815141c..97069ab988c1a9fc628cff168783a4559205e129 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
@@ -119,7 +119,7 @@ void Foam::radiation::solarLoad::updateDirectHitRadiation
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
     const scalarField& V = mesh_.V();
-    volScalarField::Boundary& QrBf = Qr_.boundaryFieldRef();
+    volScalarField::Boundary& qrBf = qr_.boundaryFieldRef();
 
     forAll(hitFacesId, i)
     {
@@ -135,7 +135,7 @@ void Foam::radiation::solarLoad::updateDirectHitRadiation
 
             for (label bandI = 0; bandI < nBands_; bandI++)
             {
-                QrBf[patchID][localFaceI] +=
+                qrBf[patchID][localFaceI] +=
                     (qPrim & n[localFaceI])
                   * spectralDistribution_[bandI]
                   * absorptivity_[patchID][bandI]()[localFaceI];
@@ -167,7 +167,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
     const scalarField& V = mesh_.V();
-    volScalarField::Boundary& QrBf = Qr_.boundaryFieldRef();
+    volScalarField::Boundary& qrBf = qr_.boundaryFieldRef();
 
     switch(solarCalc_.sunLoadModel())
     {
@@ -231,7 +231,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
                     {
                         for (label bandI = 0; bandI < nBands_; bandI++)
                         {
-                            QrBf[patchID][faceI] +=
+                            qrBf[patchID][faceI] +=
                                 (Ed + Er)
                               * spectralDistribution_[bandI]
                               * absorptivity_[patchID][bandI]()[faceI];
@@ -269,7 +269,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
                     {
                         for (label bandI = 0; bandI < nBands_; bandI++)
                         {
-                            QrBf[patchID][faceI] +=
+                            qrBf[patchID][faceI] +=
                                 solarCalc_.diffuseSolarRad()
                               * spectralDistribution_[bandI]
                               * absorptivity_[patchID][bandI]()[faceI];
@@ -555,9 +555,9 @@ void Foam::radiation::solarLoad::calculateQdiff
         }
     }
 
-    volScalarField::Boundary& QsBf = QsecondRad_.boundaryFieldRef();
+    volScalarField::Boundary& qsBf = qsecondRad_.boundaryFieldRef();
 
-    // Fill QsecondRad_
+    // Fill qsecondRad_
     label compactId = 0;
     forAll(includePatches_, i)
     {
@@ -566,7 +566,7 @@ void Foam::radiation::solarLoad::calculateQdiff
 
         if (pp.size() > 0)
         {
-            scalarField& Qrp = QsBf[patchID];
+            scalarField& qrp = qsBf[patchID];
 
             const labelList& coarsePatchFace =
                 coarseMesh_->patchFaceMap()[patchID];
@@ -578,7 +578,7 @@ void Foam::radiation::solarLoad::calculateQdiff
                 forAll(fineFaces, k)
                 {
                     label faceI = fineFaces[k];
-                    Qrp[faceI] = localqDiffusive[compactId];
+                    qrp[faceI] = localqDiffusive[compactId];
                 }
                 compactId ++;
             }
@@ -587,15 +587,15 @@ void Foam::radiation::solarLoad::calculateQdiff
 
     const scalarField& V = mesh_.V();
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    volScalarField::Boundary& QrBf = Qr_.boundaryFieldRef();
+    volScalarField::Boundary& qrBf = qr_.boundaryFieldRef();
 
     forAllConstIter(labelHashSet, includePatches, iter)
     {
         const label patchID = iter.key();
-        const scalarField& qSecond = QsecondRad_.boundaryField()[patchID];
+        const scalarField& qSecond = qsecondRad_.boundaryField()[patchID];
         if (includeMappedPatchBasePatches[patchID])
         {
-            QrBf[patchID] += qSecond;
+            qrBf[patchID] += qSecond;
         }
         else
         {
@@ -630,31 +630,31 @@ Foam::radiation::solarLoad::solarLoad(const volScalarField& T)
         )
     ),
     coarseMesh_(),
-    Qr_
+    qr_
     (
         IOobject
         (
-            "Qr",
+            "qr",
             mesh_.time().timeName(),
             mesh_,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    QsecondRad_
+    qsecondRad_
     (
         IOobject
         (
-            "QsecondRad",
+            "qsecondRad",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("QsecondRad", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qsecondRad", dimMass/pow3(dimTime), 0.0)
     ),
     hitFaces_(),
     Ru_
@@ -720,31 +720,31 @@ Foam::radiation::solarLoad::solarLoad
         )
     ),
     coarseMesh_(),
-    Qr_
+    qr_
     (
         IOobject
         (
-            "Qr",
+            "qr",
             mesh_.time().timeName(),
             mesh_,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    QsecondRad_
+    qsecondRad_
     (
         IOobject
         (
-            "QsecondRad",
+            "qsecondRad",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("QsecondRad", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qsecondRad", dimMass/pow3(dimTime), 0.0)
     ),
     hitFaces_(),
     Ru_
@@ -812,7 +812,7 @@ Foam::radiation::solarLoad::solarLoad
         )
     ),
     coarseMesh_(),
-    Qr_
+    qr_
     (
         IOobject
         (
@@ -823,20 +823,20 @@ Foam::radiation::solarLoad::solarLoad
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    QsecondRad_
+    qsecondRad_
     (
         IOobject
         (
-            "QsecondRad",
+            "qsecondRad",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("QsecondRad", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qsecondRad", dimMass/pow3(dimTime), 0.0)
     ),
     hitFaces_(),
     Ru_
@@ -942,13 +942,13 @@ void Foam::radiation::solarLoad::calculate()
     }
 
     bool facesChanged = updateHitFaces();
-    volScalarField::Boundary& QrBf = Qr_.boundaryFieldRef();
+    volScalarField::Boundary& qrBf = qr_.boundaryFieldRef();
 
     if (facesChanged)
     {
-        // Reset Ru and Qr
+        // Reset Ru and qr
         Ru_ = dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), 0.0);
-        QrBf = 0.0;
+        qrBf = 0.0;
 
         // Add direct hit radation
         const labelList& hitFacesId = hitFaces_->rayStartFaces();
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.H b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.H
index 824f97e23242d6bce8f66955af15a83f2d25a5da..52e2e9a36e2b6ab6f44f77b38cd0f17d6758d3da 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.H
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.H
@@ -105,10 +105,10 @@ private:
         autoPtr<singleCellFvMesh> coarseMesh_;
 
         //- Net radiative heat flux [W/m2]
-        volScalarField Qr_;
+        volScalarField qr_;
 
         //- Secondary solar radiative heat flux [W/m2]
-        volScalarField QsecondRad_;
+        volScalarField qsecondRad_;
 
         //- Direct hit faces Ids
         autoPtr<faceShading> hitFaces_;
@@ -143,10 +143,10 @@ private:
         //- Face-compact map
         labelListIOList visibleFaceFaces_;
 
-        //- Couple solids through mapped boundary patch using Qr (default:true)
+        //- Couple solids through mapped boundary patch using qr (default:true)
         bool solidCoupled_;
 
-        //- Couple wall patches using Qr (default:false)
+        //- Couple wall patches using qr (default:false)
         bool wallCoupled_;
 
         //- Absorptivity list
diff --git a/src/thermophysicalModels/thermophysicalPropertiesFvPatchFields/liquidProperties/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C b/src/thermophysicalModels/thermophysicalPropertiesFvPatchFields/liquidProperties/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C
index f0f9d054700fb335a7ba3fb5840c2d32204bb58e..ca59d0e9c57bf333f7d3ead5729f0103bd4f7352 100644
--- a/src/thermophysicalModels/thermophysicalPropertiesFvPatchFields/liquidProperties/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/thermophysicalPropertiesFvPatchFields/liquidProperties/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C
@@ -149,8 +149,8 @@ humidityTemperatureCoupledMixedFvPatchScalarField
     rhoName_("rho"),
     muName_("thermo:mu"),
     TnbrName_("T"),
-    QrNbrName_("none"),
-    QrName_("none"),
+    qrNbrName_("none"),
+    qrName_("none"),
     specieName_("none"),
     liquid_(nullptr),
     liquidDict_(nullptr),
@@ -189,8 +189,8 @@ humidityTemperatureCoupledMixedFvPatchScalarField
     rhoName_(psf.rhoName_),
     muName_(psf.muName_),
     TnbrName_(psf.TnbrName_),
-    QrNbrName_(psf.QrNbrName_),
-    QrName_(psf.QrName_),
+    qrNbrName_(psf.qrNbrName_),
+    qrName_(psf.qrName_),
     specieName_(psf.specieName_),
     liquid_(psf.liquid_, false),
     liquidDict_(psf.liquidDict_),
@@ -224,8 +224,8 @@ humidityTemperatureCoupledMixedFvPatchScalarField
     rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
     muName_(dict.lookupOrDefault<word>("mu", "thermo:mu")),
     TnbrName_(dict.lookupOrDefault<word>("Tnbr", "T")),
-    QrNbrName_(dict.lookupOrDefault<word>("QrNbr", "none")),
-    QrName_(dict.lookupOrDefault<word>("Qr", "none")),
+    qrNbrName_(dict.lookupOrDefault<word>("qrNbr", "none")),
+    qrName_(dict.lookupOrDefault<word>("qr", "none")),
     specieName_(dict.lookupOrDefault<word>("specie", "none")),
     liquid_(nullptr),
     liquidDict_(),
@@ -348,8 +348,8 @@ humidityTemperatureCoupledMixedFvPatchScalarField
     rhoName_(psf.rhoName_),
     muName_(psf.muName_),
     TnbrName_(psf.TnbrName_),
-    QrNbrName_(psf.QrNbrName_),
-    QrName_(psf.QrName_),
+    qrNbrName_(psf.qrNbrName_),
+    qrName_(psf.qrName_),
     specieName_(psf.specieName_),
     liquid_(psf.liquid_, false),
     liquidDict_(psf.liquidDict_),
@@ -675,26 +675,26 @@ void Foam::humidityTemperatureCoupledMixedFvPatchScalarField::updateCoeffs()
         mpp.distribute(dmHfgNbr);
     }
 
-    // Obtain Rad heat (Qr)
-    scalarField Qr(Tp.size(), 0.0);
-    if (QrName_ != "none")
+    // Obtain Rad heat (qr)
+    scalarField qr(Tp.size(), 0.0);
+    if (qrName_ != "none")
     {
-        Qr = patch().lookupPatchField<volScalarField, scalar>(QrName_);
+        qr = patch().lookupPatchField<volScalarField, scalar>(qrName_);
     }
 
-    scalarField QrNbr(Tp.size(), 0.0);
-    if (QrNbrName_ != "none")
+    scalarField qrNbr(Tp.size(), 0.0);
+    if (qrNbrName_ != "none")
     {
-        QrNbr = nbrPatch.lookupPatchField<volScalarField, scalar>(QrNbrName_);
-        mpp.distribute(QrNbr);
+        qrNbr = nbrPatch.lookupPatchField<volScalarField, scalar>(qrNbrName_);
+        mpp.distribute(qrNbr);
     }
 
     const scalarField dmHfg(dmHfgNbr + dmHfg_);
 
     const scalarField mpCpdt(mpCpTpNbr + mpCpTp_);
 
-    // Qr > 0 (heat up the wall)
-    scalarField alpha(KDeltaNbr + mpCpdt - (Qr + QrNbr)/Tp);
+    // qr > 0 (heat up the wall)
+    scalarField alpha(KDeltaNbr + mpCpdt - (qr + qrNbr)/Tp);
 
     valueFraction() = alpha/(alpha + myKDelta_);
 
@@ -742,8 +742,8 @@ void Foam::humidityTemperatureCoupledMixedFvPatchScalarField::write
     writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
     writeEntryIfDifferent<word>(os, "mu", "thermo:mu", muName_);
     writeEntryIfDifferent<word>(os, "Tnbr", "T", TnbrName_);
-    writeEntryIfDifferent<word>(os, "QrNbr", "none", QrNbrName_);
-    writeEntryIfDifferent<word>(os, "Qr", "none", QrName_);
+    writeEntryIfDifferent<word>(os, "qrNbr", "none", qrNbrName_);
+    writeEntryIfDifferent<word>(os, "qr", "none", qrName_);
 
     if (fluid_)
     {
diff --git a/src/thermophysicalModels/thermophysicalPropertiesFvPatchFields/liquidProperties/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.H b/src/thermophysicalModels/thermophysicalPropertiesFvPatchFields/liquidProperties/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.H
index 9eea39ef7783c8c78f79f5c6a437433d2d127543..3f5b1c0eead451bad20441aae85071336c87b63e 100644
--- a/src/thermophysicalModels/thermophysicalPropertiesFvPatchFields/liquidProperties/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.H
+++ b/src/thermophysicalModels/thermophysicalPropertiesFvPatchFields/liquidProperties/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.H
@@ -212,10 +212,10 @@ private:
             const word TnbrName_;
 
             //- Name of the radiative heat flux in the neighbout region
-            const word QrNbrName_;
+            const word qrNbrName_;
 
             //- Name of the radiative heat flux field
-            const word QrName_;
+            const word qrName_;
 
             //- Name of the species on which the mass transfered (default H2O)
             const word specieName_;
diff --git a/tutorials/combustion/fireFoam/LES/compartmentFire/0/T b/tutorials/combustion/fireFoam/LES/compartmentFire/0/T
index 3242ebc6bdad37284b61480540151cf1f8f18f61..090543962072f5932bc88774b8589eca0c2739c3 100644
--- a/tutorials/combustion/fireFoam/LES/compartmentFire/0/T
+++ b/tutorials/combustion/fireFoam/LES/compartmentFire/0/T
@@ -43,8 +43,8 @@ boundaryField
         type            compressible::turbulentTemperatureRadCoupledMixed;
         value           uniform 294.75;
         Tnbr            T;
-        QrNbr           none;
-        Qr              Qr;
+        qrNbr           none;
+        qr              qr;
         kappaMethod     fluidThermo;
     }
     region0_to_panelRegion_internalWallPanel_top
@@ -52,8 +52,8 @@ boundaryField
         type            compressible::turbulentTemperatureRadCoupledMixed;
         value           uniform 294.75;
         Tnbr            T;
-        QrNbr           none;
-        Qr              Qr;
+        qrNbr           none;
+        qr              qr;
         kappaMethod     fluidThermo;
     }
     region0_to_panelRegion_internalWallPanel_bottom
@@ -61,8 +61,8 @@ boundaryField
         type            compressible::turbulentTemperatureRadCoupledMixed;
         value           uniform 294.75;
         Tnbr            T;
-        QrNbr           none;
-        Qr              Qr;
+        qrNbr           none;
+        qr              qr;
         kappaMethod     fluidThermo;
     }
 }
diff --git a/tutorials/combustion/fireFoam/LES/compartmentFire/0/panelRegion/T b/tutorials/combustion/fireFoam/LES/compartmentFire/0/panelRegion/T
index 8c508d723893c3be092e84476f73be34043518cf..e787443c2000ef4940d228d2bf34e00a256b367c 100644
--- a/tutorials/combustion/fireFoam/LES/compartmentFire/0/panelRegion/T
+++ b/tutorials/combustion/fireFoam/LES/compartmentFire/0/panelRegion/T
@@ -47,8 +47,8 @@ boundaryField
         type                compressible::turbulentTemperatureRadCoupledMixed;
         neighbourField      T;
         kappaMethod         solidThermo;
-        QrNbr               Qr;
-        Qr                  none;
+        qrNbr               qr;
+        qr                  none;
         value               $internalField;
     }
 }
diff --git a/tutorials/combustion/fireFoam/LES/compartmentFire/system/controlDict b/tutorials/combustion/fireFoam/LES/compartmentFire/system/controlDict
index a3010f4c5504e25bc2e317cc1eabd7eff2b74d60..0c0df8816b6ddffd133d6954b0aba4f4686d8088 100644
--- a/tutorials/combustion/fireFoam/LES/compartmentFire/system/controlDict
+++ b/tutorials/combustion/fireFoam/LES/compartmentFire/system/controlDict
@@ -70,7 +70,7 @@ functions
         fields              (phi);
     }
 
-    wallPanel_Qin
+    wallPanel_qin
     {
         type                patchProbes;
         libs                ("libsampling.so");
@@ -84,10 +84,10 @@ functions
             (0.2  0.2 0.01)  // HF3
             (0.0  0.4 0.01)  // HF4
         );
-        fields              (Qin);
+        fields              (qin);
     }
 
-    inletQr_Qin
+    inletqr_qin
     {
         type                patchProbes;
         libs                ("libsampling.so");
@@ -103,7 +103,7 @@ functions
             (-0.02 0.0  0.02) // HF4
             (-0.02 0.0 -0.02) // HF5
         );
-        fields    (Qr Qin);
+        fields    (qr qin);
     }
 
     thermoCouple
diff --git a/tutorials/combustion/fireFoam/LES/simplePMMApanel/0/panelRegion/T b/tutorials/combustion/fireFoam/LES/simplePMMApanel/0/panelRegion/T
index 86c1260be083230c14b4050f88ffb8f3a6fc90d6..07e537337fe1a5b7935059f690c9499d3f143228 100644
--- a/tutorials/combustion/fireFoam/LES/simplePMMApanel/0/panelRegion/T
+++ b/tutorials/combustion/fireFoam/LES/simplePMMApanel/0/panelRegion/T
@@ -38,7 +38,7 @@ boundaryField
         type                fixedIncidentRadiation;
         kappaMethod         solidThermo;
         kappa               none;
-        QrIncident          uniform 60000.0;   //W
+        qrIncident          uniform 60000.0;   //W
         value               uniform 298.15;
     }
 
diff --git a/tutorials/combustion/fireFoam/LES/simplePMMApanel/constant/pyrolysisZones b/tutorials/combustion/fireFoam/LES/simplePMMApanel/constant/pyrolysisZones
index f5fd1efacc503a9dee9d6980ecc49005c85c105d..64b0ab71653bbc6a6cada3d248d3f6af3d1fb66d 100644
--- a/tutorials/combustion/fireFoam/LES/simplePMMApanel/constant/pyrolysisZones
+++ b/tutorials/combustion/fireFoam/LES/simplePMMApanel/constant/pyrolysisZones
@@ -25,7 +25,7 @@ FoamFile
 
         reactingOneDimCoeffs
         {
-            QrHSource       no; //Energy source term due in depht radiation
+            qrHSource       no; //Energy source term due in depht radiation
 
             filmCoupled     false;
 
diff --git a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict
index 6098ecaa4945e022aa9a1769980d29541438d08d..eb59adc50e9bfeff4cdd04800ddda58d6596ccab 100644
--- a/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict
+++ b/tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict
@@ -7,7 +7,7 @@
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
-    version         2;
+    version         2.0;
     format          ascii;
     class           dictionary;
     location        "system";
@@ -23,7 +23,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         7000;
+endTime         6000;
 
 deltaT          1;
 
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
index 5dc858a4d40859793b8ee4a2c409e71c0379f93b..828b16db4fc689ea030fa89f9f96831967670f43 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
@@ -54,6 +54,7 @@ PIMPLE
     nCorrectors         1;
     nNonOrthogonalCorrectors 0;
     consistent          yes;
+    SIMPLErho           yes;
 
     pMaxFactor          1.5;
     pMinFactor          0.9;
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/T b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/T
new file mode 100644
index 0000000000000000000000000000000000000000..2df93ba37eda80c3a2c2c94e2a1a55659c5aa4f6
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/T
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            fixedValue;
+        value           uniform 350;
+    }
+
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/U b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..03d1d6a8d856e001773f425f707425c765f90b21
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/U
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            noSlip;
+    }
+
+    inlet
+    {
+        type            flowRateInletVelocity;
+        massFlowRate    constant 5;
+        rhoInlet        1000;    // Guess for rho
+    }
+
+    outlet
+    {
+        type            pressureInletOutletVelocity;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/alphat b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/alphat
new file mode 100644
index 0000000000000000000000000000000000000000..61043c15b176c51a24630b5f6f8f22b9921eda30
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/alphat
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/epsilon b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..20d9e01d3fd428b66e421cc471ac3f25a6dfe970
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/epsilon
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 200;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 200;
+    }
+
+    inlet
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.005;
+        value           uniform 200;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 200;
+        value           uniform 200;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/k b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/k
new file mode 100644
index 0000000000000000000000000000000000000000..4bddfd6596532883ae70970018ab4c36036bc17d
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/k
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+
+    inlet
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.05;
+        value           uniform 1;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/nut b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/nut
new file mode 100644
index 0000000000000000000000000000000000000000..a7ffddb820aee299c106ee28be41ab84c21c18c9
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/nut
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/p b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..b4c015f22217d0f122bd94673345805f084a6201
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/p
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/thermophysicalProperties b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/thermophysicalProperties
new file mode 100644
index 0000000000000000000000000000000000000000..426ab721160866a4dbd6eda166dfb8a5e083df74
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/thermophysicalProperties
@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    properties      liquid;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    H2O;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/turbulenceProperties b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..cd2daf8229ba0b2be3dca97cab3a5c08f20b0e8a
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/turbulenceProperties
@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RAS;
+
+RAS
+{
+    RASModel        kEpsilon;
+
+    turbulence      on;
+
+    printCoeffs     on;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..9bec6a9547f35be865bbf60a0a4d8c18b0258980
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict
@@ -0,0 +1,127 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.001;
+
+vertices
+(
+    // front-plane: z = +25mm
+    // inlet region
+    (   -50  25   25)           // pt 0
+    (     0  25   25)           // pt 1
+    (   -50  75   25)           // pt 2
+    (     0  75   25)           // pt 3
+    // outlet region
+    (  -500 -75   25)           // pt 4
+    (     0 -75   25)           // pt 5
+    (  -500 -25   25)           // pt 6
+    (     0 -25   25)           // pt 7
+    // bend mid-points
+    (    25   0   25)           // pt 8
+    (    75   0   25)           // pt 9
+    // back-plane: z = -25mm
+    // inlet region
+    (   -50  25   -25)          // pt 0 + 10
+    (     0  25   -25)          // pt 1 + 10
+    (   -50  75   -25)          // pt 2 + 10
+    (     0  75   -25)          // pt 3 + 10
+    // outlet region
+    (  -500 -75   -25)          // pt 4 + 10
+    (     0 -75   -25)          // pt 5 + 10
+    (  -500 -25   -25)          // pt 7 + 10
+    (     0 -25   -25)          // pt 8 + 10
+    // bend mid-points
+    (    25   0   -25)          // pt 8 + 10
+    (    75   0   -25)          // pt 9 + 10
+);
+
+blocks
+(
+    hex (0 1 11 10  2 3 13 12) inlet  ( 20 20 20)  simpleGrading (1 1 1)
+    hex (4 5 15 14  6 7 17 16) outlet (200 20 20)  simpleGrading (1 1 1)
+
+    hex (1 8 18 11  3 9 19 13) bend1  ( 30 20 20)  simpleGrading (1 1 1)
+    hex (5 9 19 15  7 8 18 17) bend2  ( 30 20 20)  simpleGrading (1 1 1)
+);
+
+edges
+(
+   // block 2
+   arc  1  8  ( 17.678  17.678  25)
+   arc 11 18  ( 17.678  17.678 -25)
+   arc  3  9  ( 53.033  53.033  25)
+   arc 13 19  ( 53.033  53.033 -25)
+   // block 3
+   arc  7  8  ( 17.678  -17.678  25)
+   arc 17 18  ( 17.678  -17.678 -25)
+   arc  5  9  ( 53.033  -53.033  25)
+   arc 15 19  ( 53.033  -53.033 -25)
+);
+
+boundary
+(
+        // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    Default_Boundary_Region
+    {
+        type wall;
+        faces
+        (
+            // block0
+            ( 0 1 3 2 )
+            ( 11 10 12 13 )
+            ( 0 10 11 1 )
+            ( 2 3 13 12 )
+            // block1
+            ( 4 5 7 6 )
+            ( 15 14 16 17 )
+            ( 4 14 15 5 )
+            ( 6 7 17 16 )
+            // block2
+            ( 1 8 9 3 )
+            ( 18 11 13 19 )
+            ( 3 9 19 13 )
+            ( 1 11 18 8 )
+            // block3
+            ( 5 9 8 7 )
+            ( 19 15 17 18 )
+            ( 5 15 19 9 )
+            ( 7 8 18 17 )
+        );
+    }
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 2 12 10)
+        );
+    }
+
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (4 6 16 14)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/controlDict b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..07384c85ac8ed84b73dc54c7c8848c466b5037f9
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/controlDict
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     rhoPimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.5;
+
+deltaT          2e-3;
+
+writeControl    timeStep;
+
+writeInterval   10;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+graphFormat     raw;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/decomposeParDict b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..9e844a62ba7901b16f2d0c8c470ae7b5acb2504a
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method          hierarchical;
+
+simpleCoeffs
+{
+    n               (8 1 1);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (4 2 1);
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..8c9c2b631b2a4c7272f8a23eae7edb61bc216d4f
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes
@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default             Euler;
+}
+
+gradSchemes
+{
+    default             Gauss linear;
+
+    limited             cellLimited Gauss linear 1;
+}
+
+divSchemes
+{
+    default             none;
+
+    div(phi,U)          Gauss linearUpwind limited;
+    div(phi,e)          Gauss linearUpwind limited;
+    div(phi,epsilon)    Gauss linearUpwind limited;
+    div(phi,k)          Gauss linearUpwind limited;
+    div(phi,K)          Gauss linearUpwind limited;
+    div(phiv,p)         Gauss linearUpwind limited;
+
+    div(((rho*nuEff)*dev2(T(grad(U)))))      Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..173adf6c59d86acbb1d11971622a1487082f3aec
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSolution
@@ -0,0 +1,63 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "rho.*"
+    {
+        solver nthn;
+    }
+
+    "p.*"
+    {
+        solver          GAMG;
+        smoother        GaussSeidel;
+
+        tolerance       1e-7;
+        relTol          0;
+    }
+
+    "(U|e|k|epsilon).*"
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+
+        tolerance       1e-7;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nCorrectors     3;
+    nNonOrthogonalCorrectors 0;
+    pMinFactor      0.1;
+    pMaxFactor      1.5;
+
+    transonic       no;
+    consistent      no;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        ".*"            1;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/include/1DBaffle/1DTemperatureMasterBafflePatches b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/include/1DBaffle/1DTemperatureMasterBafflePatches
index 88eac454cadf670ed9c1eaa981473d4f7a20c92a..d2bc0fbe59817a0559c8b810435ad39ba816c3c3 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/include/1DBaffle/1DTemperatureMasterBafflePatches
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/include/1DBaffle/1DTemperatureMasterBafflePatches
@@ -11,7 +11,7 @@ T
     type   compressible::thermalBaffle1D<hConstSolidThermoPhysics>;
 
     thickness       uniform 0.005;  // thickness [m]
-    Qs              uniform 100;    // heat flux [W/m2]
+    qs              uniform 100;    // heat flux [W/m2]
 
     # include "1DbaffleSolidThermo"
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/T b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/T
index a7d921215a3c3703368d3346e3d966429a973fb2..97e9af0ab3fa4140c71163a52aeb515d956da1ee 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/T
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/T
@@ -58,8 +58,8 @@ boundaryField
         inletValue      uniform 300;
         Tnbr            T;
         kappaMethod     fluidThermo;
-        QrNbr           none;
-        Qr              Qr;
+        qrNbr           none;
+        qr              qr;
         kappa           none;
     }
     air_to_solid
@@ -69,8 +69,8 @@ boundaryField
         inletValue      uniform 300;
         Tnbr            T;
         kappaMethod     fluidThermo;
-        QrNbr           none;
-        Qr              Qr;
+        qrNbr           none;
+        qr              qr;
         kappa           none;
     }
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/Qr b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/qr
similarity index 98%
rename from tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/Qr
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/qr
index e86e40d3851468980199638c28c4fbd6c5a405dc..dd7559ff857282d4ff4f26c10c52b2f53baa4276 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/Qr
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/qr
@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0/air";
-    object      Qr;
+    object      qr;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/floor/T b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/floor/T
index 35ceb480f022b89ccaebd32024835753920bb2c7..209671145c8790aec4daebde7da9dbc1f56bad2a 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/floor/T
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/floor/T
@@ -64,8 +64,8 @@ boundaryField
         value           uniform 300;
         Tnbr            T;
         kappaMethod     solidThermo;
-        QrNbr           Qr;
-        Qr              none;
+        qrNbr           qr;
+        qr              none;
         kappa           none;
     }
     floor_to_solid
@@ -74,8 +74,8 @@ boundaryField
         value           uniform 300;
         Tnbr            T;
         kappaMethod     solidThermo;
-        QrNbr           none;
-        Qr              none;
+        qrNbr           none;
+        qr              none;
         kappa           none;
     }
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/solid/T b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/solid/T
index 4e1f26a8c48a2fe050afe249eb90e5f745e78282..196552bde6470c131271205930b6ef4bd8c58a48 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/solid/T
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/solid/T
@@ -32,8 +32,8 @@ boundaryField
         value           uniform 300;
         Tnbr            T;
         kappaMethod     solidThermo;
-        QrNbr           Qr;
-        Qr              none;
+        qrNbr           qr;
+        qr              none;
         kappa           none;
     }
     solid_to_floor
@@ -42,8 +42,8 @@ boundaryField
         value           uniform 300;
         Tnbr            T;
         kappaMethod     solidThermo;
-        QrNbr           none;
-        Qr              none;
+        qrNbr           none;
+        qr              none;
         kappa           none;
     }
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun.pre b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun.pre
index 37a630197599f4e696c7b42ec2acf6e8625d9ead..78752106d3d6bd53217090e120b38a004a5bbb6d 100755
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun.pre
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun.pre
@@ -16,7 +16,7 @@ rm -r system/domain3
 # remove fluid fields from solid regions (important for post-processing)
 for i in solid floor
 do
-   rm -f 0*/$i/{rho,mut,alphat,epsilon,k,U,p_rgh,Qr,G,IDefault}
+   rm -f 0*/$i/{rho,mut,alphat,epsilon,k,U,p_rgh,qr,G,IDefault}
 done
 
 for i in air solid floor
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/constant/air/radiationProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/constant/air/radiationProperties
index dd4e6f8da4a19a7c08072d511acac88629670bc1..fe50c86eeacdd19c3f8e816f1001eb9a6fce7ee5 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/constant/air/radiationProperties
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/constant/air/radiationProperties
@@ -72,8 +72,8 @@ solarLoadCoeffs
         C   0.058;  // Model constant
 
     // Radiative flux coupling flags
-    solidCoupled    true;  //Couple through Qr the solid regions (default true)
-    wallCoupled     false; //Couple through Qr wall patches (default false)
+    solidCoupled    true;  //Couple through qr the solid regions (default true)
+    wallCoupled     false; //Couple through qr wall patches (default false)
 }
 
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/air/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/air/changeDictionaryDict
index 33ce84e5598af295e95415c33fab6d2a861601b9..46208629875248ea9a0723a09c2d3026cd632907 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/air/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/air/changeDictionaryDict
@@ -69,8 +69,8 @@ dictionaryReplacement
                 type            compressible::turbulentTemperatureRadCoupledMixed;
                 Tnbr            T;
                 kappaMethod     fluidThermo;
-                QrNbr           none;
-                Qr              Qr;
+                qrNbr           none;
+                qr              qr;
                 kappa           none;
                 value           uniform 300;
             }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/floor/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/floor/changeDictionaryDict
index 6ae161489a78d4bb9c31f0a2871d7720ce3ca02b..3860f7dc3fcc24b619a2f3a19ac13d94d2d8e0ef 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/floor/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/floor/changeDictionaryDict
@@ -40,8 +40,8 @@ dictionaryReplacement
                 type            compressible::turbulentTemperatureRadCoupledMixed;
                 Tnbr            T;
                 kappaMethod     solidThermo;
-                QrNbr           none;
-                Qr              none;
+                qrNbr           none;
+                qr              none;
                 kappa           none;
                 value           uniform 300;
             }
@@ -51,8 +51,8 @@ dictionaryReplacement
                 type            compressible::turbulentTemperatureRadCoupledMixed;
                 Tnbr            T;
                 kappaMethod     solidThermo;
-                QrNbr           Qr;
-                Qr              none;
+                qrNbr           qr;
+                qr              none;
                 kappa           none;
                 value           uniform 300;
             }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/solid/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/solid/changeDictionaryDict
index 1a3e5fe65c845c29b3d24f04f445a2eadb64a591..080b5fba86def0d2b8101ace168d0f0e7217dfbb 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/solid/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/solid/changeDictionaryDict
@@ -34,8 +34,8 @@ dictionaryReplacement
                 type            compressible::turbulentTemperatureRadCoupledMixed;
                 Tnbr            T;
                 kappaMethod     solidThermo;
-                QrNbr           Qr;
-                Qr              none;
+                qrNbr           qr;
+                qr              none;
                 kappa           none;
                 value           uniform 300;
             }
@@ -45,8 +45,8 @@ dictionaryReplacement
                 type            compressible::turbulentTemperatureRadCoupledMixed;
                 Tnbr            T;
                 kappaMethod     solidThermo;
-                QrNbr           none;
-                Qr              none;
+                qrNbr           none;
+                qr              none;
                 kappa           none;
                 value           uniform 300;
             }