diff --git a/applications/solvers/multiphase/bubbleFoam/pEqn.H b/applications/solvers/multiphase/bubbleFoam/pEqn.H
index 631c9cbba6..af3f6e2822 100644
--- a/applications/solvers/multiphase/bubbleFoam/pEqn.H
+++ b/applications/solvers/multiphase/bubbleFoam/pEqn.H
@@ -14,6 +14,12 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
+    U.boundaryField() =
+        alpha1.boundaryField()*U1.boundaryField()
+      + alpha2.boundaryField()*U2.boundaryField();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
+
     surfaceScalarField phiDrag1
     (
         fvc::interpolate(alpha2/rho1*dragCoef*rAU1)*phi2
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 9e67a47c4f..e0c13b16e0 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -29,6 +29,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
index 709d30c398..915e714082 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
@@ -22,6 +22,12 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
+    U.boundaryField() =
+        alpha1.boundaryField()*U1.boundaryField()
+      + alpha2.boundaryField()*U2.boundaryField();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
+
     surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
 
     if (g0.value() > 0.0)
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
index 9c56c4d26f..28de546ce0 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
@@ -5,6 +5,11 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
+    mrfZones.absoluteFlux(phi);
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
+    mrfZones.relativeFlux(phi);
+
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H
index d1bc7bf1ed..a21d5b7bc8 100644
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@@ -4,6 +4,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index e3ef0dc9ef..aa3ed1fe85 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -4,6 +4,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
index 7bc1eeb12c..5e87fa0319 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
@@ -5,6 +5,11 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
+    mrfZones.absoluteFlux(phi);
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
+    mrfZones.relativeFlux(phi);
+
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
index c1e346ce43..2e26fb538a 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@@ -4,6 +4,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H
index 9f27e992a8..99c15656dc 100644
--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@@ -5,6 +5,9 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(rho.boundaryField()*U.boundaryField())
+      & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
index ddfca4e3ea..810a753c9d 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@@ -4,6 +4,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 9b543bcd25..d6bde3812b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -18,6 +18,13 @@
     mrfZones.absoluteFlux(phi1);
     mrfZones.absoluteFlux(phi2.oldTime());
     mrfZones.absoluteFlux(phi2);
+    mrfZones.absoluteFlux(phi);
+
+    U.boundaryField() =
+        alpha1.boundaryField()*U1.boundaryField()
+      + alpha2.boundaryField()*U2.boundaryField();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
 
@@ -56,6 +63,7 @@
     mrfZones.relativeFlux(phi1);
     mrfZones.relativeFlux(phi2.oldTime());
     mrfZones.relativeFlux(phi2);
+    mrfZones.relativeFlux(phi);
 
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H.old b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H.old
deleted file mode 100644
index 348cc847d7..0000000000
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H.old
+++ /dev/null
@@ -1,118 +0,0 @@
-{
-    surfaceScalarField alpha1f(fvc::interpolate(alpha1));
-    surfaceScalarField alpha2f(scalar(1) - alpha1f);
-
-    volScalarField rAU1(1.0/U1Eqn.A());
-    volScalarField rAU2(1.0/U2Eqn.A());
-
-    rAU1f = fvc::interpolate(rAU1);
-    surfaceScalarField rAU2f(fvc::interpolate(rAU2));
-
-    volVectorField HbyA1("HbyA1", U1);
-    HbyA1 = rAU1*U1Eqn.H();
-
-    volVectorField HbyA2("HbyA2", U2);
-    HbyA2 = rAU2*U2Eqn.H();
-
-    mrfZones.absoluteFlux(phi1.oldTime());
-    mrfZones.absoluteFlux(phi1);
-
-    mrfZones.absoluteFlux(phi2.oldTime());
-    mrfZones.absoluteFlux(phi2);
-
-    surfaceScalarField phiDrag1
-    (
-        fvc::interpolate(alpha2/rho1*K*rAU1)*phi2 + rAU1f*(g & mesh.Sf())
-    );
-
-    if (g0.value() > 0.0)
-    {
-        phiDrag1 -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
-    }
-
-    if (kineticTheory.on())
-    {
-        phiDrag1 -= rAU1f*fvc::snGrad(kineticTheory.pa()/rho1)*mesh.magSf();
-    }
-
-
-    surfaceScalarField phiDrag2
-    (
-        fvc::interpolate(alpha1/rho2*K*rAU2)*phi1 + rAU2f*(g & mesh.Sf())
-    );
-
-    // Fix for gravity on outlet boundary.
-    forAll(p.boundaryField(), patchi)
-    {
-        if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
-        {
-            phiDrag1.boundaryField()[patchi] = 0.0;
-            phiDrag2.boundaryField()[patchi] = 0.0;
-        }
-    }
-
-    surfaceScalarField phiHbyA1
-    (
-        "phiHbyA1",
-        (fvc::interpolate(HbyA1) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU1, U1, phi1)
-      + phiDrag1
-    );
-    mrfZones.relativeFlux(phiHbyA1);
-
-    surfaceScalarField phiHbyA2
-    (
-        "phiHbyA2",
-        (fvc::interpolate(HbyA2) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU2, U2, phi2)
-      + phiDrag2
-    );
-    mrfZones.relativeFlux(phiHbyA2);
-
-    mrfZones.relativeFlux(phi1.oldTime());
-    mrfZones.relativeFlux(phi1);
-    mrfZones.relativeFlux(phi2.oldTime());
-    mrfZones.relativeFlux(phi2);
-
-    surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
-
-    surfaceScalarField Dp
-    (
-        "Dp",
-        alpha1f*rAU1f/rho1 + alpha2f*rAU2f/rho2
-    );
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
-        );
-
-        pEqn.setReference(pRefCell, pRefValue);
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            surfaceScalarField SfGradp(pEqn.flux()/Dp);
-
-            phi1 = phiHbyA1 - rAU1f*SfGradp/rho1;
-            phi2 = phiHbyA2 - rAU2f*SfGradp/rho2;
-            phi = alpha1f*phi1 + alpha2f*phi2;
-
-            p.relax();
-            SfGradp = pEqn.flux()/Dp;
-
-            U1 = HbyA1 + fvc::reconstruct(phiDrag1 - rAU1f*SfGradp/rho1);
-            U1.correctBoundaryConditions();
-
-            U2 = HbyA2 + fvc::reconstruct(phiDrag2 - rAU2f*SfGradp/rho2);
-            U2.correctBoundaryConditions();
-
-            U = alpha1*U1 + alpha2*U2;
-        }
-    }
-}
-
-#include "continuityErrs.H"
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
index f4f4f54f65..1cb9326e32 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
@@ -333,7 +333,7 @@ tmp<FieldField<fvsPatchField, Type> > interpolate
     const tmp<FieldField<fvPatchField, Type> >& tfvpff
 )
 {
-    tmp<FieldField<fvPatchField, Type> > tfvspff = interpolate(tfvpff());
+    tmp<FieldField<fvsPatchField, Type> > tfvspff = interpolate(tfvpff());
     tfvpff.clear();
     return tfvspff;
 }