diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H
index 1576f6ba787b037d97c84368836912f0c7d129dd..138e58fc7f72b23abdcba18deef3da53e3b6b951 100644
--- a/applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H
@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
+                  - fvc::snGrad(p)
                 ) * mesh.magSf()
             )
         );
diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H b/applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H
index 4152105784639530f725c53b7b561841b798c624..d82a03edb59a2d419420700638404e6f90b197c8 100644
--- a/applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H
@@ -12,7 +12,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 
diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H b/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H
index dc04fb454680bbeadec55c8946227270601651c7..6b8e67cc4920d3da95a06c22b81db6dee407bc84 100644
--- a/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H
@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -88,24 +88,6 @@
 
     dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        )
-    );
-
     volScalarField rho1 = rho10 + psi1*p;
     volScalarField rho2 = rho20 + psi2*p;
 
@@ -152,11 +134,11 @@
     );
 
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H
index 7e4b37061fdede3fccc625875cc38c447657a198..e6004eb9de9f6d157dd5a42637540fa3d70c6ddb 100644
--- a/applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H
@@ -2,17 +2,17 @@
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rUAf = fvc::interpolate(rUA);
 
-    tmp<fvScalarMatrix> pdEqnComp;
+    tmp<fvScalarMatrix> pEqnComp;
 
     if (transonic)
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
     }
     else
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
     }
 
 
@@ -26,16 +26,16 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqnIncomp
+        fvScalarMatrix pEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, pd)
+          - fvm::laplacian(rUAf, p)
         );
 
         if
@@ -51,9 +51,9 @@
                     max(alpha1, scalar(0))*(psi1/rho1)
                   + max(alpha2, scalar(0))*(psi2/rho2)
                 )
-               *pdEqnComp()
-              + pdEqnIncomp,
-                mesh.solver(pd.name() + "Final")
+               *pEqnComp()
+              + pEqnIncomp,
+                mesh.solver(p.name() + "Final")
             );
         }
         else
@@ -64,8 +64,8 @@
                     max(alpha1, scalar(0))*(psi1/rho1)
                   + max(alpha2, scalar(0))*(psi2/rho2)
                 )
-               *pdEqnComp()
-              + pdEqnIncomp
+               *pEqnComp()
+              + pEqnIncomp
             );
         }
 
@@ -73,26 +73,21 @@
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(pdEqnComp & pd);
-            phi += pdEqnIncomp.flux();
+               *(pEqnComp & p);
+            phi += pEqnIncomp.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 
-    p = max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        );
+    p.max(pMin);
 
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(pd) " << min(pd).value() << endl;
+    Info<< "min(p) " << min(p).value() << endl;
 
     // Make the fluxes relative to the mesh motion
     fvc::makeRelative(phi, U);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
index 528e0aaafd892c2995fa94848c0d59b66af41942..0b1a9ac029d776dbd57da20f91c14da8ee88f5a7 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
+                  - fvc::snGrad(p)
                 ) * mesh.magSf()
             )
         );
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 1f579d245bdcac21f647ac2a90a516b5f612b3c9..3e6904d383e677c976ac8daacd110b31af414f75 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -46,11 +46,6 @@
     #include "createPhi.H"
 
 
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi);
 
@@ -88,24 +83,6 @@
 
     dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        )
-    );
-
     volScalarField rho1 = rho10 + psi1*p;
     volScalarField rho2 = rho20 + psi2*p;
 
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index ebf24498ade0bd1d1572d6055da8b5bccd369ffa..9d2dc23916babb9c6c94acd2b53125c8cfee1dad 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -2,17 +2,17 @@
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rUAf = fvc::interpolate(rUA);
 
-    tmp<fvScalarMatrix> pdEqnComp;
+    tmp<fvScalarMatrix> pEqnComp;
 
     if (transonic)
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
     }
     else
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
     }
 
 
@@ -26,16 +26,16 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf;
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqnIncomp
+        fvScalarMatrix pEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, pd)
+          - fvm::laplacian(rUAf, p)
         );
 
         solve
@@ -44,31 +44,27 @@
                 max(alpha1, scalar(0))*(psi1/rho1)
               + max(alpha2, scalar(0))*(psi2/rho2)
             )
-           *pdEqnComp()
-          + pdEqnIncomp
+           *pEqnComp()
+          + pEqnIncomp
         );
 
         if (nonOrth == nNonOrthCorr)
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(pdEqnComp & pd);
-            phi += pdEqnIncomp.flux();
+               *(pEqnComp & p);
+            phi += pEqnIncomp.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 
-    p = max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        );
+    p.max(pMin);
 
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(pd) " << min(pd).value() << endl;
+    Info<< "min(p) " << min(p).value() << endl;
 }
diff --git a/applications/solvers/multiphase/interDyMFoam/correctPhi.H b/applications/solvers/multiphase/interDyMFoam/correctPhi.H
index c975c9b37416e10a066bb578ac05b531a4217473..1f7845f347d25a2e9e45cdbd6539e1998022f5e4 100644
--- a/applications/solvers/multiphase/interDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interDyMFoam/correctPhi.H
@@ -12,7 +12,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 
@@ -27,7 +27,7 @@
             fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
         );
 
-        pcorrEqn.setReference(pdRefCell, pdRefValue);
+        pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
         if (nonOrth == nNonOrthCorr)
diff --git a/applications/solvers/multiphase/interDyMFoam/createFields.H b/applications/solvers/multiphase/interDyMFoam/createFields.H
index bcceb9d748193843a3124b585b538497b61e9d22..26644e17f472fe4dd8605280c0ed8e5a1ee0a5dd 100644
--- a/applications/solvers/multiphase/interDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/interDyMFoam/createFields.H
@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -92,47 +92,17 @@
         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
     }
 
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*(g & mesh.C())
-    );
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
-
+    label pRefCell = 0;
     scalar pRefValue = 0.0;
-
-    if (pd.needReference())
-    {
-        pRefValue = readScalar
-        (
-            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
-        );
-
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pdRefCell)
-        );
-    }
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
diff --git a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
index 046503a8c459abc2f1c6caeb958c39a2ee985668..cf8ec8eb323bdc42352669aae04249baa07426f1 100644
--- a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
@@ -114,18 +114,6 @@ int main(int argc, char *argv[])
             #include "pEqn.H"
         }
 
-        p = pd + rho*gh;
-
-        if (pd.needReference())
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pdRefCell)
-            );
-        }
-
         turbulence->correct();
 
         runTime.write();
diff --git a/applications/solvers/multiphase/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interDyMFoam/pEqn.H
index f8ff18f60631c3cd66f02a8b9b4e8adc5ceb19e9..8825661c8307d70f2c8cb8bf6f18675579aaf6cb 100644
--- a/applications/solvers/multiphase/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interDyMFoam/pEqn.H
@@ -7,38 +7,38 @@
 
     phi = phiU +
     (
-        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-      - ghf*fvc::snGrad(rho)
-    )*rAUf*mesh.magSf();
+        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+      + fvc::interpolate(rho)*(g & mesh.Sf())
+    )*rAUf;
 
-    if (pd.needReference())
+    if (p.needReference())
     {
         fvc::makeRelative(phi, U);
-        adjustPhi(phi, U, pd);
+        adjustPhi(phi, U, p);
         fvc::makeAbsolute(phi, U);
     }
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAUf, pd) == fvc::div(phi)
+            fvm::laplacian(rAUf, p) == fvc::div(phi)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi -= pdEqn.flux();
+            phi -= pEqn.flux();
         }
     }
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
index c59137c7b87a9be3c799e03567dc8209341798f5..757e37a2cfdc7d9dbdacd0621120452219ac2a2c 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
@@ -25,10 +25,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
+                  - fvc::snGrad(p)
                 ) * mesh.magSf()
             )
         );
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H b/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H
index 171e1670f47dbdf29d3904c47af7d196b706296e..642aa1c5d2a88c157d1d704f555c9b87ddefcc48 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H
@@ -1,11 +1,11 @@
 {
 #   include "continuityErrs.H"
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
@@ -22,7 +22,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 
@@ -37,7 +37,7 @@
             fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
         );
 
-        pcorrEqn.setReference(pdRefCell, pdRefValue);
+        pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
         if (nonOrth == nNonOrthCorr)
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
index 5e727dea3eca4828580cadaf9a10531a7719fd62..98d166214b2ca5d34fa2aa554dfab679e16df7ee 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -66,26 +66,9 @@
     rho.oldTime();
 
 
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
-
-    Info<< "Calculating field g.h" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*gh
-    );
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
 
 
     // Construct interface from alpha1 distribution
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index 0037d71cf52ce5dd50dd298767b27987d4e03e38..4f290157f64bbf3de676acc26c9f5e708af91ba6 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -13,11 +13,11 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf;
 
-    adjustPhi(phi, U, pd);
+    adjustPhi(phi, U, p);
 
     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
     const volScalarField& vDotcP = vDotP[0]();
@@ -25,31 +25,29 @@
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvc::div(phi) - fvm::laplacian(rUAf, pd)
-          + (vDotvP - vDotcP)*(rho*gh - pSat) + fvm::Sp(vDotvP - vDotcP, pd)
+            fvc::div(phi) - fvm::laplacian(rUAf, p)
+          - (vDotvP - vDotcP)*pSat + fvm::Sp(vDotvP - vDotcP, p)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi += pdEqn.flux();
+            phi += pEqn.flux();
         }
     }
 
-    p = pd + rho*gh;
-
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 }
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
index 6efcdbc8b1fe550c9b9ee5d2f24b470751ca347e..9d7e36dc95a0e695946f20a2e0fb4d9f80512b56 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@@ -12,9 +12,9 @@
 
     phi = phiU +
         (
-            mixture.surfaceTensionForce()
+            mixture.surfaceTensionForce()*mesh.magSf()
           + fvc::interpolate(rho)*(g & mesh.Sf())
-        )*rUAf*mesh.magSf();
+        )*rUAf;
 
     adjustPhi(phi, U, p);
 
diff --git a/applications/solvers/multiphase/settlingFoam/UEqn.H b/applications/solvers/multiphase/settlingFoam/UEqn.H
index ac4c18ab19dc4b70fe0770f23e0a24bc4af3eebd..04e9194363d9ee4edc6aaa65e032fcebcaeddfca 100644
--- a/applications/solvers/multiphase/settlingFoam/UEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/UEqn.H
@@ -22,10 +22,8 @@
           ==
             fvc::reconstruct
             (
-                (
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p)
-                )*mesh.magSf()
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              - fvc::snGrad(p)*mesh.magSf()
             )
         );
     }
diff --git a/applications/solvers/multiphase/settlingFoam/createFields.H b/applications/solvers/multiphase/settlingFoam/createFields.H
index b13649c3124ddaf5ba06c5e13e3b090b746798ec..7c8445a71d4465bae5514c86afbf83f6611ce057 100644
--- a/applications/solvers/multiphase/settlingFoam/createFields.H
+++ b/applications/solvers/multiphase/settlingFoam/createFields.H
@@ -337,6 +337,3 @@
         ),
         mut + mul
     );
-
-    Info<< "Calculating field (g.h)f\n" << endl;
-    surfaceScalarField ghf = surfaceScalarField("ghf", g & mesh.Cf());
diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H
index 65b954c6da84e1fe6b98ee5a99f295f7a14ff324..b4fc05b32ab56f04646427310478e9cc92429b11 100644
--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@@ -15,7 +15,7 @@ phi =
     );
 
 surfaceScalarField phiU("phiU", phi);
-phi -= ghf*fvc::snGrad(rho)*rUAf*mesh.magSf();
+phi += fvc::interpolate(rho)*(g & mesh.Sf())*rUAf;
 
 for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {