diff --git a/applications/solvers/multiphase/interFoam/alphaCourantNo.H b/applications/solvers/multiphase/VoF/alphaCourantNo.H
similarity index 96%
rename from applications/solvers/multiphase/interFoam/alphaCourantNo.H
rename to applications/solvers/multiphase/VoF/alphaCourantNo.H
index a084155c8ae70e1e0b4a6cc9ae9b6f5fba828235..24f6e48a22cb58366d5cfec831d0720ded233f86 100644
--- a/applications/solvers/multiphase/interFoam/alphaCourantNo.H
+++ b/applications/solvers/multiphase/VoF/alphaCourantNo.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
diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/VoF/alphaEqn.H
similarity index 86%
rename from applications/solvers/multiphase/interFoam/alphaEqn.H
rename to applications/solvers/multiphase/VoF/alphaEqn.H
index 312c1d6723df78c8853f033e7d2d341efc605ab3..570442a304914253f103a22e9740298073298879 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/VoF/alphaEqn.H
@@ -79,6 +79,8 @@
 
     if (MULESCorr)
     {
+        #include "alphaSuSp.H"
+
         fvScalarMatrix alpha1Eqn
         (
             (
@@ -92,6 +94,8 @@
                 phiCN,
                 upwind<scalar>(mesh, phiCN)
             ).fvmDiv(phiCN, alpha1)
+         ==
+            Su + fvm::Sp(Sp + divU, alpha1)
         );
 
         alpha1Eqn.solve();
@@ -124,6 +128,8 @@
 
     for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
     {
+        #include "alphaSuSp.H"
+
         surfaceScalarField phir(phic*mixture.nHatf());
 
         tmp<surfaceScalarField> talphaPhiUn
@@ -154,7 +160,17 @@
             tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
             volScalarField alpha10("alpha10", alpha1);
 
-            MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr.ref(), 1, 0);
+            MULES::correct
+            (
+                geometricOneField(),
+                alpha1,
+                talphaPhiUn(),
+                talphaPhiCorr.ref(),
+                Sp,
+                (-Sp*alpha1)(),
+                1,
+                0
+            );
 
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
@@ -171,7 +187,17 @@
         {
             alphaPhi = talphaPhiUn;
 
-            MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0);
+            MULES::explicitSolve
+            (
+                geometricOneField(),
+                alpha1,
+                phiCN,
+                alphaPhi,
+                Sp,
+                (Su + divU*min(alpha1(), scalar(1)))(),
+                1,
+                0
+            );
         }
 
         alpha2 = 1.0 - alpha1;
@@ -195,7 +221,8 @@
      == fv::EulerDdtScheme<vector>::typeName
     )
     {
-        rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2;
+        #include "rhofs.H"
+        rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f;
     }
     else
     {
@@ -206,7 +233,8 @@
         }
 
         // Calculate the end-of-time-step mass flux
-        rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
+        #include "rhofs.H"
+        rhoPhi = alphaPhi*(rho1f - rho2f) + phi*rho2f;
     }
 
     Info<< "Phase-1 volume fraction = "
diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/VoF/alphaEqnSubCycle.H
similarity index 100%
rename from applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
rename to applications/solvers/multiphase/VoF/alphaEqnSubCycle.H
diff --git a/applications/solvers/multiphase/interFoam/setDeltaT.H b/applications/solvers/multiphase/VoF/setDeltaT.H
similarity index 95%
rename from applications/solvers/multiphase/interFoam/setDeltaT.H
rename to applications/solvers/multiphase/VoF/setDeltaT.H
index 9cc860c032eed2e9a46138518589b83f1e8aed19..924d24c8ea728f2543f192ee5dfd9e5425ffa7b1 100644
--- a/applications/solvers/multiphase/interFoam/setDeltaT.H
+++ b/applications/solvers/multiphase/VoF/setDeltaT.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/applications/solvers/multiphase/interFoam/setRDeltaT.H b/applications/solvers/multiphase/VoF/setRDeltaT.H
similarity index 100%
rename from applications/solvers/multiphase/interFoam/setRDeltaT.H
rename to applications/solvers/multiphase/VoF/setRDeltaT.H
diff --git a/applications/solvers/multiphase/compressibleInterFoam/Make/options b/applications/solvers/multiphase/compressibleInterFoam/Make/options
index 27e73ef53aa6b642f5e16f8aa0a7929b383a1102..b947836f65eda540a6d54e4770c67f263be75458 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/Make/options
+++ b/applications/solvers/multiphase/compressibleInterFoam/Make/options
@@ -1,6 +1,6 @@
 EXE_INC = \
     -I. \
-    -I../interFoam \
+    -I../VoF \
     -ItwoPhaseMixtureThermo \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
diff --git a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
index 2eaeaa25f4b5dab15477d015be5458bd26a58d18..b5e66ec33fb5194fa8d154aa79d8351008b2f501 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
@@ -1,7 +1,7 @@
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
-      + fvm::div(rhoPhi, U)
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         fvOptions(rho, U)
diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaControls.H b/applications/solvers/multiphase/compressibleInterFoam/alphaControls.H
deleted file mode 100644
index a99a0b39e289f60589a61d3c43b1fb94fcff81dc..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/compressibleInterFoam/alphaControls.H
+++ /dev/null
@@ -1,5 +0,0 @@
-const dictionary& alphaControls = mesh.solverDict(alpha1.name());
-
-label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
-
-label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaEqn.H b/applications/solvers/multiphase/compressibleInterFoam/alphaEqn.H
deleted file mode 100644
index 8105650296f4ca0221dc19bcce249cb4a7228a13..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/compressibleInterFoam/alphaEqn.H
+++ /dev/null
@@ -1,51 +0,0 @@
-{
-    word alphaScheme("div(phi,alpha)");
-    word alpharScheme("div(phirb,alpha)");
-
-    surfaceScalarField phir(phic*mixture.nHatf());
-
-    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
-    {
-        #include "alphaSuSp.H"
-
-        surfaceScalarField alphaPhi1
-        (
-            fvc::flux
-            (
-                phi,
-                alpha1,
-                alphaScheme
-            )
-          + fvc::flux
-            (
-                -fvc::flux(-phir, alpha2, alpharScheme),
-                alpha1,
-                alpharScheme
-            )
-        );
-
-        MULES::explicitSolve
-        (
-            geometricOneField(),
-            alpha1,
-            phi,
-            alphaPhi1,
-            Sp,
-            Su,
-            1,
-            0
-        );
-
-        surfaceScalarField rho1f(fvc::interpolate(rho1));
-        surfaceScalarField rho2f(fvc::interpolate(rho2));
-        rhoPhi = alphaPhi1*(rho1f - rho2f) + phi*rho2f;
-
-        alpha2 = scalar(1) - alpha1;
-    }
-
-    Info<< "Liquid phase volume fraction = "
-        << alpha1.weightedAverage(mesh.V()).value()
-        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
-        << "  Min(" << alpha2.name() << ") = " << min(alpha2).value()
-        << endl;
-}
diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/compressibleInterFoam/alphaEqnSubCycle.H
deleted file mode 100644
index 54141d34a7875bc10865929a033fbb8188f10f86..0000000000000000000000000000000000000000
--- a/applications/solvers/multiphase/compressibleInterFoam/alphaEqnSubCycle.H
+++ /dev/null
@@ -1,38 +0,0 @@
-{
-    surfaceScalarField phic(mag(phi/mesh.magSf()));
-    phic = min(mixture.cAlpha()*phic, max(phic));
-
-    volScalarField divU(fvc::div(fvc::absolute(phi, U)));
-
-    if (nAlphaSubCycles > 1)
-    {
-        dimensionedScalar totalDeltaT = runTime.deltaT();
-        surfaceScalarField rhoPhiSum
-        (
-            IOobject
-            (
-                "rhoPhiSum",
-                runTime.timeName(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar("0", rhoPhi.dimensions(), 0)
-        );
-
-        for
-        (
-            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
-            !(++alphaSubCycle).end();
-        )
-        {
-            #include "alphaEqn.H"
-            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
-        }
-
-        rhoPhi = rhoPhiSum;
-    }
-    else
-    {
-        #include "alphaEqn.H"
-    }
-}
diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H b/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H
index a1b383e89e20032ac01a0e13a12463d3cbf389bf..81309cd091ffbed46ba2ce456c9afaf110394f22 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H
@@ -1,37 +1,43 @@
-        volScalarField::Internal Sp
-        (
-            IOobject
-            (
-                "Sp",
-                runTime.timeName(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
-        );
+volScalarField::Internal Sp
+(
+    IOobject
+    (
+        "Sp",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("Sp", dgdt.dimensions(), 0)
+);
 
-        volScalarField::Internal Su
-        (
-            IOobject
-            (
-                "Su",
-                runTime.timeName(),
-                mesh
-            ),
-            // Divergence term is handled explicitly to be
-            // consistent with the explicit transport solution
-            divU*min(alpha1, scalar(1))
-        );
+volScalarField::Internal Su
+(
+    IOobject
+    (
+        "Su",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("Su", dgdt.dimensions(), 0)
+);
 
-        forAll(dgdt, celli)
-        {
-            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
-            {
-                Sp[celli] -= dgdt[celli]*alpha1[celli];
-                Su[celli] += dgdt[celli]*alpha1[celli];
-            }
-            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
-            {
-                Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
-            }
-        }
+forAll(dgdt, celli)
+{
+    if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
+    {
+        Sp[celli] -= dgdt[celli]*alpha1[celli];
+        Su[celli] += dgdt[celli]*alpha1[celli];
+    }
+    else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
+    {
+        Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
+    }
+}
+
+volScalarField::Internal divU
+(
+    mesh.moving()
+  ? fvc::div(phiCN() + mesh.phi())
+  : fvc::div(phiCN())
+);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H
new file mode 100644
index 0000000000000000000000000000000000000000..f8d0a15836077f9505224bc5acf8b5966dddf2e2
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H
@@ -0,0 +1,3 @@
+{
+    #include "alphaEqnSubCycle.H"
+}
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options
index 73e50e5ef644884216dbf43426820f34007b3dec..b44af44a7231a401a9421a4c641e382e69b5bec7 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options
@@ -1,6 +1,7 @@
 EXE_INC = \
+    -I. \
     -I.. \
-    -I../../interFoam \
+    -I../../VoF \
     -I../twoPhaseMixtureThermo \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index 9f09b35184bd885602e8979c1d0a1666d2c38dfb..f82665a49e44292e193317e6d7d4804ceafc6b8d 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -39,7 +39,10 @@ Description
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
-#include "MULES.H"
+#include "CMULES.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
 #include "subCycle.H"
 #include "twoPhaseMixture.H"
 #include "twoPhaseMixtureThermo.H"
@@ -47,7 +50,6 @@ Description
 #include "pimpleControl.H"
 #include "fvOptions.H"
 #include "CorrectPhi.H"
-#include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -141,7 +143,7 @@ int main(int argc, char *argv[])
             }
 
             #include "alphaControls.H"
-            #include "alphaEqnSubCycle.H"
+            #include "compressibleAlphaEqnSubCycle.H"
 
             solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index d0dbfada5634f7fd997b8cc0ff6b1a3eb3cfad73..7733f6a28880acddb92c7b572955f5e87ae56b7f 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -8,6 +8,7 @@
         fvc::flux(HbyA)
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
     );
+    MRF.makeRelative(phiHbyA);
 
     surfaceScalarField phig
     (
@@ -20,7 +21,7 @@
     phiHbyA += phig;
 
     // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, U, phiHbyA, rAUf);
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
 
     // Make the fluxes relative to the mesh motion
     fvc::makeRelative(phiHbyA, U);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index bcc28de1bdcb0996359ad6a1b69ea5c505d94f67..f6ec5d20cd2bca42b0bc60b6a1447dbab94fa08c 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -36,7 +36,10 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "MULES.H"
+#include "CMULES.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
 #include "subCycle.H"
 #include "rhoThermo.H"
 #include "twoPhaseMixture.H"
@@ -44,7 +47,6 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 7af8f15fe4fd6ccb7e8dfdcd852f7799361044ab..9bac4d825922be2cae8a78f1f656001b49136c99 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -63,6 +63,7 @@ dimensionedScalar pMin
 );
 
 mesh.setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(alpha1.name());
 
 
 #include "readGravitationalAcceleration.H"
@@ -99,3 +100,22 @@ autoPtr<compressible::turbulenceModel> turbulence
 
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
+
+// MULES flux from previous time-step
+surfaceScalarField alphaPhi
+(
+    IOobject
+    (
+        "alphaPhi",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    phi*fvc::interpolate(alpha1)
+);
+
+// MULES Correction
+tmp<surfaceScalarField> talphaPhiCorr0;
+
+#include "createMRF.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index a40a98fd0ace0b284c7a420804e051c2313af3b5..51c58ba53689d0b3582c9c0407c81d0136b90f01 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -8,6 +8,7 @@
         fvc::flux(HbyA)
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
+    MRF.makeRelative(phiHbyA);
 
     surfaceScalarField phig
     (
@@ -20,7 +21,7 @@
     phiHbyA += phig;
 
     // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, U, phiHbyA, rAUf);
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
 
     tmp<fvScalarMatrix> p_rghEqnComp1;
     tmp<fvScalarMatrix> p_rghEqnComp2;
diff --git a/applications/solvers/multiphase/compressibleInterFoam/rhofs.H b/applications/solvers/multiphase/compressibleInterFoam/rhofs.H
new file mode 100644
index 0000000000000000000000000000000000000000..67dc10f378441716045250b114dca8ab6db486fe
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleInterFoam/rhofs.H
@@ -0,0 +1,2 @@
+surfaceScalarField rho1f(fvc::interpolate(rho1));
+surfaceScalarField rho2f(fvc::interpolate(rho2));
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
index f3010e0e57fbc734a79ed4618c39b70a0dd28775..9df3add9e22497c1bdb24393900f4d56004d9997 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I. \
+    -I../VoF \
     -I../interFoam \
     -ImultiphaseMixtureThermo/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
diff --git a/applications/solvers/multiphase/interFoam/Make/options b/applications/solvers/multiphase/interFoam/Make/options
index cfe4a920729ed8bd1497f89de334caf51b85a312..7f2fed1d162beb222aaf0653c179cd1b3e843b2c 100644
--- a/applications/solvers/multiphase/interFoam/Make/options
+++ b/applications/solvers/multiphase/interFoam/Make/options
@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I../VoF \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
diff --git a/applications/solvers/multiphase/interFoam/alphaSuSp.H b/applications/solvers/multiphase/interFoam/alphaSuSp.H
new file mode 100644
index 0000000000000000000000000000000000000000..0b8e05e18710e49d98c8942cb771ee274ae1786f
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/alphaSuSp.H
@@ -0,0 +1,3 @@
+zeroField Su;
+zeroField Sp;
+zeroField divU;
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/Make/options b/applications/solvers/multiphase/interFoam/interDyMFoam/Make/options
index bceea20c8147c9bf4b7ff720f6d00f6ba21b65bd..3d028feaf7ee7b395b33738f6763bcf407bd1b7d 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/Make/options
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/Make/options
@@ -1,6 +1,7 @@
 EXE_INC = \
     -I. \
     -I.. \
+    -I../../VoF \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index 2c7df14ae7864b32f598331db27d1b3cacefaef9..aafc10bcc457beba7bd6ab0dccfaa84cea662bd0 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@@ -28,7 +28,7 @@
     phiHbyA += phig;
 
     // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, U, phiHbyA, rAUf);
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
 
     while (pimple.correctNonOrthogonal())
     {
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/Make/options b/applications/solvers/multiphase/interFoam/interMixingFoam/Make/options
index ce7222ff1bcf43af0ef9ab7a62e7d43b469072b3..a07545b21637aa42b030b81a719722aae244a5d9 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/Make/options
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/Make/options
@@ -1,6 +1,7 @@
 EXE_INC = \
     -I. \
     -I.. \
+    -I../../VoF \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -IimmiscibleIncompressibleThreePhaseMixture \
     -IincompressibleThreePhaseMixture \
diff --git a/applications/solvers/multiphase/interFoam/rhofs.H b/applications/solvers/multiphase/interFoam/rhofs.H
new file mode 100644
index 0000000000000000000000000000000000000000..5949bf1e0a60d21934f09e1fbe8780c5cf7e43d5
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/rhofs.H
@@ -0,0 +1,2 @@
+const dimensionedScalar& rho1f(rho1);
+const dimensionedScalar& rho2f(rho2);
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/Make/options b/applications/solvers/multiphase/multiphaseInterFoam/Make/options
index 99850d69be78f2de6c487d44459a19e5d43aee51..576e7aa224c16f8d2c0eda66dc38e8f5ed744523 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/Make/options
+++ b/applications/solvers/multiphase/multiphaseInterFoam/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I. \
+    -I../VoF \
     -I../interFoam \
     -ImultiphaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels \
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterDyMFoam/Make/options b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterDyMFoam/Make/options
index 1a210c20e682b6da1647db40cbb548907985dd8b..b7d051f1ca8d83e173ce73644787401c3038b239 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterDyMFoam/Make/options
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterDyMFoam/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I.. \
+    -I../../VoF \
     -I../../interFoam/interDyMFoam \
     -I../../interFoam \
     -I../multiphaseMixture/lnInclude \
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/Make/options b/applications/solvers/multiphase/twoLiquidMixingFoam/Make/options
index 4dfc63419d1f824c5509433bbf7f45cdbefd25f8..e9549744cea30e721d7f48650648dfb2e3ae70c5 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/Make/options
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/Make/options
@@ -1,6 +1,6 @@
 EXE_INC = \
     -I. \
-    -I../interFoam \
+    -I../VoF \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
diff --git a/src/OpenFOAM/fields/Fields/zeroField/zeroField.H b/src/OpenFOAM/fields/Fields/zeroField/zeroField.H
index f8e53959e631fa94c4fb34170122e8f579b3c5f0..15f25f504db759b30dedd4ff1c2815a86925f4a0 100644
--- a/src/OpenFOAM/fields/Fields/zeroField/zeroField.H
+++ b/src/OpenFOAM/fields/Fields/zeroField/zeroField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,6 +66,10 @@ public:
         inline scalar operator[](const label) const;
 
         inline zeroField field() const;
+
+        inline zeroField operator()() const;
+
+        inline zeroField operator-() const;
 };
 
 
diff --git a/src/OpenFOAM/fields/Fields/zeroField/zeroFieldI.H b/src/OpenFOAM/fields/Fields/zeroField/zeroFieldI.H
index fb63a48bf9e13e7764ee7d4fadf6fe650d724824..59287815fd55259380c13407362806c15e5f1808 100644
--- a/src/OpenFOAM/fields/Fields/zeroField/zeroFieldI.H
+++ b/src/OpenFOAM/fields/Fields/zeroField/zeroFieldI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,4 +39,80 @@ inline Foam::zeroField Foam::zeroField::field() const
 }
 
 
+inline Foam::zeroField Foam::zeroField::operator()() const
+{
+    return zeroField();
+}
+
+
+inline Foam::zeroField Foam::zeroField::operator-() const
+{
+    return zeroField();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline const zeroField operator+(const zeroField&, const zeroField&)
+{
+    return zeroField();
+}
+
+template<class Type>
+inline const Type& operator+(const Type& t, const zeroField&)
+{
+    return t;
+}
+
+template<class Type>
+inline const Type& operator+(const zeroField&, const Type& t)
+{
+    return t;
+}
+
+inline const zeroField operator-(const zeroField&, const zeroField&)
+{
+    return zeroField();
+}
+
+template<class Type>
+inline const Type& operator-(const Type& t, const zeroField&)
+{
+    return t;
+}
+
+template<class Type>
+inline Type operator-(const zeroField&, const Type& t)
+{
+    return -t;
+}
+
+template<class Type>
+inline zeroField operator*(const Type& t, const zeroField&)
+{
+    return zeroField();
+}
+
+template<class Type>
+inline zeroField operator*(const zeroField&, const Type& t)
+{
+    return zeroField();
+}
+
+template<class Type>
+inline zeroField operator/(const zeroField&, const Type& t)
+{
+    return zeroField();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroField.H b/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroField.H
index 0eb7c87235d0cfc6fb79e7d10599d8383ba18c30..2a9b4fe591a84ff83c9dea70acd8e6d47e101222 100644
--- a/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroField.H
+++ b/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroField.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
@@ -72,6 +72,8 @@ public:
 
         inline zeroField field() const;
 
+        inline zeroField operator()() const;
+
         inline zeroField oldTime() const;
 
         inline zeroFieldField boundaryField() const;
@@ -84,7 +86,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    #include "geometricZeroFieldI.H"
+#include "geometricZeroFieldI.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroFieldI.H b/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroFieldI.H
index d6188702d58937a8d260261dd045e66395021dd6..ca5d89f3f05b5a2ece4f3326acb712632ef94c2e 100644
--- a/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroFieldI.H
+++ b/src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroFieldI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,6 +42,11 @@ inline Foam::zeroField Foam::geometricZeroField::field() const
     return zeroField();
 }
 
+inline Foam::zeroField Foam::geometricZeroField::operator()() const
+{
+    return zeroField();
+}
+
 inline Foam::zeroField Foam::geometricZeroField::oldTime() const
 {
     return zeroField();
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/controlDict b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/controlDict
index c1d09089d0eedca9370e03b2ce22eb91da531058..dc0b7630a109f1aa5961af24bf96b2b614198fff 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/controlDict
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/controlDict
@@ -29,15 +29,15 @@ deltaT          0.0001;
 
 writeControl    adjustableRunTime;
 
-writeInterval   0.005;
+writeInterval   0.01;
 
 purgeWrite      0;
 
-writeFormat     ascii;
+writeFormat     binary;
 
 writePrecision  8;
 
-writeCompression compressed;
+writeCompression off;
 
 timeFormat      general;
 
@@ -47,10 +47,9 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.25;
-
+maxCo           0.5;
+maxAlphaCo      0.5;
 maxDeltaT       1;
-maxAlphaCo      1;
 
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution
index 7c97abb05a8db73aeb72559796c9ad08574c353b..041112222d93c29f80ee09327873fb576b3cf3d7 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSolution
@@ -17,11 +17,19 @@ FoamFile
 
 solvers
 {
-    alpha.water
+    "alpha.water.*"
     {
         nAlphaCorr      1;
-        nAlphaSubCycles 1;
+        nAlphaSubCycles 2;
         cAlpha          1;
+
+        MULESCorr       no;
+        nLimiterIter    5;
+
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-8;
+        relTol          0;
     }
 
     ".*(rho|rhoFinal)"
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/controlDict b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/controlDict
index bb61973589ff9314fc924f3b85d2f1511fe52e7a..a049876aa2b077e88dc943a40ad3319a4adbc890 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/controlDict
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/controlDict
@@ -29,7 +29,7 @@ deltaT          0.0001;
 
 writeControl    adjustableRunTime;
 
-writeInterval   0.005;
+writeInterval   0.01;
 
 purgeWrite      0;
 
@@ -47,10 +47,9 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.25;
-
+maxCo           0.5;
+maxAlphaCo      0.5;
 maxDeltaT       1;
-maxAlphaCo      1;
 
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSolution b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSolution
index f75504c67558f70d71463f857efb2ee41520e563..dac0efd3b8bcc8be1822aeca219f8520e7cde83d 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSolution
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSolution
@@ -17,11 +17,19 @@ FoamFile
 
 solvers
 {
-    alpha.water
+    "alpha.water.*"
     {
         nAlphaCorr      1;
-        nAlphaSubCycles 1;
+        nAlphaSubCycles 2;
         cAlpha          1;
+
+        MULESCorr       no;
+        nLimiterIter    5;
+
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-8;
+        relTol          0;
     }
 
     pcorr