diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 2f6333cc9761b89834245073ac3edfad744a318f..bd12a82495b7b8bb8ee00733d5cc374c6b44a0d6 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@@ -91,10 +91,7 @@ K = 0.5*magSqr(U);
 {
     rhoUf = fvc::interpolate(rho*U);
     surfaceVectorField n(mesh.Sf()/mesh.magSf());
-    rhoUf +=
-        mesh.Sf()
-       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
-       /sqr(mesh.magSf());
+    rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
 }
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index fd892bb9327e766841defe9190cfedaa49d382b1..1c58a02c33950e10054b8a9748c25278d939f114 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -107,10 +107,7 @@ K = 0.5*magSqr(U);
 {
     rhoUf = fvc::interpolate(rho*U);
     surfaceVectorField n(mesh.Sf()/mesh.magSf());
-    rhoUf +=
-        mesh.Sf()
-       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
-       /sqr(mesh.magSf());
+    rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
 }
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
index 65afd1aeb88642cae2b880ac21ba6951bab763f0..4799f6456728144ee8a1327dc2fbb064b2cd8df0 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
@@ -42,8 +42,5 @@ U.correctBoundaryConditions();
 {
     rhoUf = fvc::interpolate(rho*U);
     surfaceVectorField n(mesh.Sf()/mesh.magSf());
-    rhoUf +=
-        mesh.Sf()
-       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
-       /sqr(mesh.magSf());
+    rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
 }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index 7ecfe1f66236e7458835eb35dde5a76a7038b573..e2938bdb949776932b8e15d11e405fc52cb26553 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
@@ -51,7 +51,7 @@ fvOptions.correct(U);
 {
     Uf = fvc::interpolate(U);
     surfaceVectorField n(mesh.Sf()/mesh.magSf());
-    Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+    Uf += n*(phi/mesh.magSf() - (n & Uf));
 }
 
 // Make the fluxes relative to the mesh motion
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
index b10c24aa0b7c86eb281cea29e612929b31e1ad50..d23ec0f95058a2ea46ee547ba9833b2e4851e9e9 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
@@ -93,10 +93,7 @@ K = 0.5*magSqr(U);
 {
     rhoUf = fvc::interpolate(rho*U);
     surfaceVectorField n(mesh.Sf()/mesh.magSf());
-    rhoUf +=
-        mesh.Sf()
-       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
-       /sqr(mesh.magSf());
+    rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
 }
 
 if (thermo.dpdt())
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index 0f0d5768414323fb1df135a80b7bd1b5e27d2832..a0e0e6c74dc740175847af7d89fd41388a3dbb38 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@@ -85,6 +85,6 @@
     {
         Uf = fvc::interpolate(U);
         surfaceVectorField n(mesh.Sf()/mesh.magSf());
-        Uf += mesh.Sf()*(phiv - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+        Uf += n*(phiv/mesh.magSf() - (n & Uf));
     }
 }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index c86921371f52ec96f0cff265257110280f460546..a014a5ab4edcc2ca517b2552ed655bca42297b6c 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -114,7 +114,7 @@
     {
         Uf = fvc::interpolate(U);
         surfaceVectorField n(mesh.Sf()/mesh.magSf());
-        Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+        Uf += n*(phi/mesh.magSf() - (n & Uf));
     }
 
     // Make the fluxes relative to the mesh motion
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index 0a3d7ed11e97ed60f03d44a5b8e61705d303c4f5..ba7b4efa9ae1094b0e14f3c5a43d7a76868a3102 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -55,10 +55,10 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
-    #include "createUf.H"
     #include "readTimeControls.H"
     #include "createPrghCorrTypes.H"
     #include "correctPhi.H"
+    #include "createUf.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index c30b361e72f9badfc98172483c455fc469552275..76441207919fe4427f6734750ff4de1a26c544f4 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@@ -67,7 +67,7 @@
     {
         Uf = fvc::interpolate(U);
         surfaceVectorField n(mesh.Sf()/mesh.magSf());
-        Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+        Uf += n*(phi/mesh.magSf() - (n & Uf));
     }
 
     // Make the fluxes relative to the mesh motion
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
index 0bfd8ad56eab6c4ca652b238ec5f4e8c36a3847b..810edeff308fb8ee57549e504f71aeb2a53b7ac4 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
@@ -65,10 +65,10 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
-    #include "createUf.H"
     #include "readTimeControls.H"
     #include "createPcorrTypes.H"
     #include "../interFoam/interDyMFoam/correctPhi.H"
+    #include "createUf.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
index 7863f5d74851ebba806897821174ff43b6484bea..6fac1df95c2236b94da9914f01a8de5416ed230f 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
@@ -68,7 +68,7 @@
     {
         Uf = fvc::interpolate(U);
         surfaceVectorField n(mesh.Sf()/mesh.magSf());
-        Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+        Uf += n*(phi/mesh.magSf() - (n & Uf));
     }
 
     // Make the fluxes relative to the mesh motion
diff --git a/src/finiteVolume/cfdTools/compressible/createRhoUf.H b/src/finiteVolume/cfdTools/compressible/createRhoUf.H
index e91115ae75b2a5b9a6e6f932352c8129bc4b693f..312aec60dca9cdaffbd448bd17b3053e9731a3fe 100644
--- a/src/finiteVolume/cfdTools/compressible/createRhoUf.H
+++ b/src/finiteVolume/cfdTools/compressible/createRhoUf.H
@@ -46,9 +46,14 @@ surfaceVectorField rhoUf
         IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
-    linearInterpolate(rho*U)
+    fvc::interpolate(rho*U)
 );
 
+{
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    rhoUf += n*(phi/mesh.magSf() - (n & rhoUf));
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/src/finiteVolume/cfdTools/incompressible/createUf.H b/src/finiteVolume/cfdTools/incompressible/createUf.H
index aab92ce0e2aec37304bd634f72bcbee68ed80189..0b65f2daed23a820fd4850cea3852e61d4ddf21d 100644
--- a/src/finiteVolume/cfdTools/incompressible/createUf.H
+++ b/src/finiteVolume/cfdTools/incompressible/createUf.H
@@ -46,9 +46,15 @@ surfaceVectorField Uf
         IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
-    linearInterpolate(U)
+    fvc::interpolate(U)
 );
 
+{
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    Uf += n*(phi/mesh.magSf() - (n & Uf));
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif