diff --git a/applications/solvers/heatTransfer/buoyantFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantFoam/UEqn.H
index 988268c4ef..24cfbf4f68 100644
--- a/applications/solvers/heatTransfer/buoyantFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantFoam/UEqn.H
@@ -9,4 +9,18 @@
 
     UEqn.relax();
 
-    solve(UEqn == -fvc::grad(pd) - fvc::grad(rho)*gh);
+    if (momentumPredictor)
+    {
+        solve
+        (
+            UEqn
+         ==
+           -fvc::reconstruct
+            (
+                (
+                    fvc::snGrad(pd)
+                  + ghf*fvc::snGrad(rho)
+                ) * mesh.magSf()
+            )
+        );
+    }
diff --git a/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C b/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C
index cafb677900..8b3cca8d21 100644
--- a/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C
+++ b/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C
@@ -41,13 +41,12 @@ Description
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -55,23 +54,24 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
-//#       include "setDeltaT.H"
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
 
         for (int corr=0; corr<nCorr; corr++)
         {
-#           include "hEqn.H"
-#           include "pEqn.H"
+            #include "hEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();
diff --git a/applications/solvers/heatTransfer/buoyantFoam/createFields.H b/applications/solvers/heatTransfer/buoyantFoam/createFields.H
index 3d604d4daf..aa32325c34 100644
--- a/applications/solvers/heatTransfer/buoyantFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantFoam/createFields.H
@@ -57,6 +57,7 @@
 
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("gh", g & mesh.Cf());
 
     dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
 
diff --git a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H
index 36378f7f6f..9feededcd3 100644
--- a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H
@@ -1,60 +1,67 @@
-bool closedVolume = pd.needReference();
-
-rho = thermo->rho();
-
-volScalarField rUA = 1.0/UEqn.A();
-U = rUA*UEqn.H();
-
-phi =
-    fvc::interpolate(rho)
-   *(
-        (fvc::interpolate(U) & mesh.Sf())
-      + fvc::ddtPhiCorr(rUA, rho, U, phi)
-    )
-  - fvc::interpolate(rho*rUA*gh)*fvc::snGrad(rho)*mesh.magSf();
-
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
-    fvScalarMatrix pdEqn
+    bool closedVolume = pd.needReference();
+
+    rho = thermo->rho();
+
+    volScalarField rUA = 1.0/UEqn.A();
+    surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA);
+
+    U = rUA*UEqn.H();
+
+    surfaceScalarField phiU
     (
-        fvm::ddt(psi, pd)
-      + fvc::ddt(psi)*pRef
-      + fvc::ddt(psi, rho)*gh
-      + fvc::div(phi)
-      - fvm::laplacian(rho*rUA, pd)
+        fvc::interpolate(rho)
+       *(
+           (fvc::interpolate(U) & mesh.Sf())
+         + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
     );
 
-    if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+    phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        pdEqn.solve(mesh.solver(pd.name() + "Final"));
-    }
-    else
-    {
-        pdEqn.solve(mesh.solver(pd.name()));
+        fvScalarMatrix pdEqn
+        (
+            fvm::ddt(psi, pd)
+          + fvc::ddt(psi)*pRef
+          + fvc::ddt(psi, rho)*gh
+          + fvc::div(phi)
+          - fvm::laplacian(rhorUAf, pd)
+        );
+
+        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        {
+            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+        }
+        else
+        {
+            pdEqn.solve(mesh.solver(pd.name()));
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pdEqn.flux();
+        }
     }
 
-    if (nonOrth == nNonOrthCorr)
+    U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
+    U.correctBoundaryConditions();
+
+    p == pd + rho*gh + pRef;
+    dpdt = fvc::ddt(p);
+
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
+
+    // For closed-volume cases adjust the pressure and density levels
+    // to obey overall mass continuity
+    if (closedVolume)
     {
-        phi += pdEqn.flux();
+        p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
+            /fvc::domainIntegrate(thermo->psi());
+        rho = thermo->rho();
     }
-}
 
-p == pd + rho*gh + pRef;
-dpdt = fvc::ddt(p);
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
-
-U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
-U.correctBoundaryConditions();
-
-
-// For closed-volume cases adjust the pressure and density levels
-// to obey overall mass continuity
-if (closedVolume)
-{
-    p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
-        /fvc::domainIntegrate(thermo->psi());
     pd == p - (rho*gh + pRef);
-    rho = thermo->rho();
 }
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
index 5bc4c4738c..71a57cd2a9 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
@@ -11,7 +11,15 @@
 
     eqnResidual = solve
     (
-        UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh
+        UEqn()
+     ==
+       -fvc::reconstruct
+        (
+            (
+                fvc::snGrad(pd)
+              + ghf*fvc::snGrad(rho)
+            ) * mesh.magSf()
+        )
     ).initialResidual();
 
     maxResidual = max(eqnResidual, maxResidual);
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index 184242be81..d26da5cb32 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -53,6 +53,7 @@
 
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("gh", g & mesh.Cf());
 
     dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index 8c3621205b..3c257d249d 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -1,55 +1,65 @@
-volScalarField rUA = 1.0/UEqn().A();
-U = rUA*UEqn().H();
-UEqn.clear();
-
-phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
-bool closedVolume = adjustPhi(phi, U, p);
-phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf();
-
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
-    fvScalarMatrix pdEqn
-    (
-        fvm::laplacian(rho*rUA, pd) == fvc::div(phi)
-    );
+    volScalarField rUA = 1.0/UEqn().A();
+    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
 
-    pdEqn.setReference(pdRefCell, pdRefValue);
-    // retain the residual from the first iteration
-    if (nonOrth == 0)
+    U = rUA*UEqn().H();
+    UEqn.clear();
+
+    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
+    bool closedVolume = adjustPhi(phi, U, p);
+    surfaceScalarField buoyancyPhi = ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
+    phi -= buoyancyPhi;
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        eqnResidual = pdEqn.solve().initialResidual();
-        maxResidual = max(eqnResidual, maxResidual);
-    }
-    else
-    {
-        pdEqn.solve();
+        fvScalarMatrix pdEqn
+        (
+            fvm::laplacian(rhorUAf, pd) == fvc::div(phi)
+        );
+
+        pdEqn.setReference(pdRefCell, pdRefValue);
+
+        // retain the residual from the first iteration
+        if (nonOrth == 0)
+        {
+            eqnResidual = pdEqn.solve().initialResidual();
+            maxResidual = max(eqnResidual, maxResidual);
+        }
+        else
+        {
+            pdEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            // Calculate the conservative fluxes
+            phi -= pdEqn.flux();
+
+            // Explicitly relax pressure for momentum corrector
+            pd.relax();
+
+            // Correct the momentum source with the pressure gradient flux
+            // calculated from the relaxed pressure
+            U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rhorUAf);
+            U.correctBoundaryConditions();
+        }
     }
 
-    if (nonOrth == nNonOrthCorr)
+    #include "continuityErrs.H"
+
+    p == pd + rho*gh + pRef;
+
+    // For closed-volume cases adjust the pressure and density levels
+    // to obey overall mass continuity
+    if (closedVolume)
     {
-        phi -= pdEqn.flux();
+        p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
+            /fvc::domainIntegrate(thermo->psi());
     }
-}
 
-#include "continuityErrs.H"
+    rho = thermo->rho();
+    rho.relax();
+    Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
 
-// Explicitly relax pressure for momentum corrector
-pd.relax();
-
-p = pd + rho*gh + pRef;
-
-U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
-U.correctBoundaryConditions();
-
-// For closed-volume cases adjust the pressure and density levels
-// to obey overall mass continuity
-if (closedVolume)
-{
-    p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
-        /fvc::domainIntegrate(thermo->psi());
     pd == p - (rho*gh + pRef);
 }
-
-rho = thermo->rho();
-rho.relax();
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options
index 22d41fdb80..89d7787af6 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options
@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I../buoyantSimpleFoam \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H
deleted file mode 100644
index 5bc4c4738c..0000000000
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H
+++ /dev/null
@@ -1,18 +0,0 @@
-    // Solve the Momentum equation
-
-    tmp<fvVectorMatrix> UEqn
-    (
-        fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
-      + turbulence->divDevRhoReff(U)
-    );
-
-    UEqn().relax();
-
-    eqnResidual = solve
-    (
-        UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh
-    ).initialResidual();
-
-    maxResidual = max(eqnResidual, maxResidual);
-
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H
index 76e3805bca..62c06ec38d 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H
@@ -54,6 +54,7 @@
 
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("gh", g & mesh.Cf());
 
     dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H
deleted file mode 100644
index 2a713bac62..0000000000
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H
+++ /dev/null
@@ -1,54 +0,0 @@
-volScalarField rUA = 1.0/UEqn().A();
-U = rUA*UEqn().H();
-UEqn.clear();
-phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
-bool closedVolume = adjustPhi(phi, U, p);
-phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf();
-
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-{
-    fvScalarMatrix pdEqn
-    (
-        fvm::laplacian(rho*rUA, pd) == fvc::div(phi)
-    );
-
-    pdEqn.setReference(pdRefCell, pdRefValue);
-    // retain the residual from the first iteration
-    if (nonOrth == 0)
-    {
-        eqnResidual = pdEqn.solve().initialResidual();
-        maxResidual = max(eqnResidual, maxResidual);
-    }
-    else
-    {
-        pdEqn.solve();
-    }
-
-    if (nonOrth == nNonOrthCorr)
-    {
-        phi -= pdEqn.flux();
-    }
-}
-
-#include "continuityErrs.H"
-
-// Explicitly relax pressure for momentum corrector
-pd.relax();
-
-p = pd + rho*gh + pRef;
-
-U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
-U.correctBoundaryConditions();
-
-// For closed-volume cases adjust the pressure and density levels
-// to obey overall mass continuity
-if (closedVolume)
-{
-    p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
-        /fvc::domainIntegrate(thermo->psi());
-    pd == p - (rho*gh + pRef);
-}
-
-rho = thermo->rho();
-rho.relax();
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;