diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index c7289b23f901906c4fe7c96366c041ac8f0a3e9d..c4edc961d887231dfa7be0afe4f105f39379862c 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -70,7 +70,7 @@
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
index 602affd4e00fbc1644b5decb6faf3f3f3451c448..12dc93795c250ca173b90710b80e8c9bc84dd4f3 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
@@ -83,14 +83,14 @@ Foam::multiphaseMixtureThermo::multiphaseMixtureThermo
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("rho*phi", dimMass/dimTime, 0.0)
+        dimensionedScalar("rhoPhi", dimMass/dimTime, 0.0)
     ),
 
     alphas_
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
index 317a18120c18394070602caed5845f7b186d0999..e9aa9fb9b9efa54295d25c629fa901ca72bc9073 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
@@ -77,19 +77,26 @@ int main(int argc, char *argv[])
 
         #include "setrDeltaT.H"
 
-        twoPhaseProperties.correct();
-
-        #define LTSSOLVE
-        #include "alphaEqnSubCycle.H"
-        #undef LTSSOLVE
-
-        interface.correct();
-
-        turbulence->correct();
+        tmp<surfaceScalarField> tphiAlpha;
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties.correct();
+
+                #define LTSSOLVE
+                #include "alphaEqnSubCycle.H"
+                #undef LTSSOLVE
+
+                interface.correct();
+            }
+
+            turbulence->correct();
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
index 0fa1a50ee2f6613bc8dbd9e4b3b4800b5e0ca9c4..56527781a25b4f3b5093f04ccd5c370438e48635 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
@@ -81,15 +81,22 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        twoPhaseProperties.correct();
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
-        #include "zonePhaseVolumes.H"
+        tmp<surfaceScalarField> tphiAlpha;
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties.correct();
+
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+                #include "zonePhaseVolumes.H"
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop
diff --git a/applications/solvers/multiphase/interFoam/Make/options b/applications/solvers/multiphase/interFoam/Make/options
index b1cfcac9bd812c22bd62a4560d5c5b2d42ac20e6..8811a53df650d1ae1646edae2769d78846fbe4fa 100644
--- a/applications/solvers/multiphase/interFoam/Make/options
+++ b/applications/solvers/multiphase/interFoam/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -ggdb3 \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H
index b0dd8ebef2abd5bd3e59dc21ce4cc18b8267821c..f6e57df771d897333bcc161ab381f66f4767f4b5 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqn.H
@@ -6,9 +6,7 @@
     phic = min(interface.cAlpha()*phic, max(phic));
     surfaceScalarField phir(phic*interface.nHatf());
 
-    tmp<surfaceScalarField> tphiAlpha;
-
-    if (MULESCorr)
+    if (pimple.firstIter() && MULESCorr)
     {
         fvScalarMatrix alpha1Eqn
         (
diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
index 81c65c1caf737480956f72ff2d3952f7311436c5..5db9f16546180ca8b428eb2527dd5798c4acddba 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
@@ -1,5 +1,3 @@
-#include "alphaControls.H"
-
 if (nAlphaSubCycles > 1)
 {
     dimensionedScalar totalDeltaT = runTime.deltaT();
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index f0d18098483634633a10f2c7ffb9068a06f91cb6..3039171f77a74c38d46bb5770ed85d8357f3d597 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -62,13 +62,14 @@
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        rho1*phi
+        mesh,
+        dimensionedScalar("0", dimMass/dimTime, 0)
     );
 
 
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index 4a673c6a3f3770f07b7bfcd5937c4d48254ca9ad..d56e300581d9d9e4cbd46dc79270392f4d51f887 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -77,46 +77,56 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+        tmp<surfaceScalarField> tphiAlpha;
 
-        mesh.update();
-
-        if (mesh.changing())
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
         {
-            Info<< "Execution time for mesh.update() = "
-                << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
-                << " s" << endl;
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
-            gh = g & mesh.C();
-            ghf = g & mesh.Cf();
-        }
+                mesh.update();
 
-        if (mesh.changing() && correctPhi)
-        {
-            // Calculate absolute flux from the mapped surface velocity
-            phi = mesh.Sf() & Uf;
+                if (mesh.changing())
+                {
+                    Info<< "Execution time for mesh.update() = "
+                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                        << " s" << endl;
 
-            #include "correctPhi.H"
+                    gh = g & mesh.C();
+                    ghf = g & mesh.Cf();
+                }
 
-            // Make the flux relative to the mesh motion
-            fvc::makeRelative(phi, U);
+                if (mesh.changing() && correctPhi)
+                {
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & Uf;
 
-            interface.correct();
-        }
+                    #include "correctPhi.H"
 
-        if (mesh.changing() && checkMeshCourantNo)
-        {
-            #include "meshCourantNo.H"
-        }
+                    // Make the flux relative to the mesh motion
+                    fvc::makeRelative(phi, U);
 
-        twoPhaseProperties.correct();
+                    interface.correct();
+                }
 
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
+                if (mesh.changing() && checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
+            }
+
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties.correct();
+
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+            }
 
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
             #include "UEqn.H"
 
             // --- Pressure corrector loop
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/readControls.H b/applications/solvers/multiphase/interFoam/interDyMFoam/readControls.H
index d4e332ff38e9d2699c0919ef45a1a9edd3769fc1..bf8b38f42616a82904a1aa8ff4dd5eba9cffb629 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/readControls.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/readControls.H
@@ -1,6 +1,16 @@
-    #include "readTimeControls.H"
+#include "readTimeControls.H"
 
-    bool correctPhi =
-        pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
-    bool checkMeshCourantNo =
-        pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);
+bool correctPhi
+(
+    pimple.dict().lookupOrDefault<Switch>("correctPhi", true)
+);
+
+bool checkMeshCourantNo
+(
+    pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false)
+);
+
+bool moveMeshOuterCorrectors
+(
+    pimple.dict().lookupOrDefault<Switch>("moveMeshOuterCorrectors", false)
+);
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index 4ce08d27db070538ee81555549c40272704a79a4..0ddfed50396eefbe37953f16d608fa60b69d83f1 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -80,14 +80,21 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        twoPhaseProperties.correct();
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
+        tmp<surfaceScalarField> tphiAlpha;
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties.correct();
+
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H b/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
index 196c82194c9fdac753607d426491c67f5dbaf972..3fd5874d5b479efc9288e9bf143970c2227a481e 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
@@ -63,7 +63,7 @@
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
index 12bc57d785f7fc5af819bde912225b69c7a3edcf..3321da57f9dffa699ab0a136187f94c34810f49f 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
@@ -74,16 +74,23 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        threePhaseProperties.correct();
-
-        #include "alphaEqnsSubCycle.H"
-        interface.correct();
-
-        #define twoPhaseProperties threePhaseProperties
+        tmp<surfaceScalarField> tphiAlpha;
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                threePhaseProperties.correct();
+
+                #include "alphaEqnsSubCycle.H"
+                interface.correct();
+
+                #define twoPhaseProperties threePhaseProperties
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop
diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
index 3423a4654affdbe01bea8824f4c693ec0a2b0be9..15b4439220846d0079623443d49ef38b0530ac86 100644
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
@@ -83,14 +83,21 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        twoPhaseProperties.correct();
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
+        tmp<surfaceScalarField> tphiAlpha;
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties.correct();
+
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
index a4338b907f8c0f633d89919d4896d276d21cc16e..077738ccb2945dc7b3609742c86a0974434bb2e1 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
@@ -1,18 +1,4 @@
-surfaceScalarField rhoPhi
-(
-    IOobject
-    (
-        "rhoPhi",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("0", dimMass/dimTime, 0)
-);
-
 {
-    #include "alphaControls.H"
-
     surfaceScalarField phic(mag(phi/mesh.magSf()));
     phic = min(interface.cAlpha()*phic, max(phic));
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
index 31fa619754e1537253ec4ee1783f8b5d92929730..0bfd8ad56eab6c4ca652b238ec5f4e8c36a3847b 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
@@ -86,44 +86,64 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
-        mesh.update();
+                mesh.update();
 
-        if (mesh.changing())
-        {
-            Info<< "Execution time for mesh.update() = "
-                << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
-                << " s" << endl;
+                if (mesh.changing())
+                {
+                    Info<< "Execution time for mesh.update() = "
+                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                        << " s" << endl;
 
-            gh = g & mesh.C();
-            ghf = g & mesh.Cf();
-        }
+                    gh = g & mesh.C();
+                    ghf = g & mesh.Cf();
+                }
 
-        if (mesh.changing() && correctPhi)
-        {
-            // Calculate absolute flux from the mapped surface velocity
-            phi = mesh.Sf() & Uf;
+                if (mesh.changing() && correctPhi)
+                {
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & Uf;
 
-            #include "../interFoam/interDyMFoam/correctPhi.H"
+                    #include "../interFoam/interDyMFoam/correctPhi.H"
 
-            // Make the flux relative to the mesh motion
-            fvc::makeRelative(phi, U);
-        }
+                    // Make the flux relative to the mesh motion
+                    fvc::makeRelative(phi, U);
+                }
 
-        if (mesh.changing() && checkMeshCourantNo)
-        {
-            #include "meshCourantNo.H"
-        }
+                if (mesh.changing() && checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
+            }
 
-        twoPhaseProperties->correct();
+            #include "alphaControls.H"
+
+            surfaceScalarField rhoPhi
+            (
+                IOobject
+                (
+                    "rhoPhi",
+                    runTime.timeName(),
+                    mesh
+                ),
+                mesh,
+                dimensionedScalar("0", dimMass/dimTime, 0)
+            );
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties->correct();
 
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+            }
 
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
             #include "UEqn.H"
 
             // --- Pressure corrector loop
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
index 28ef7f0a69166534f9860cdd8e55f7ea4fa7c768..2a653dcf37fb115d9ad1dc5c4de7552f168bea0f 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
@@ -83,14 +83,31 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        twoPhaseProperties->correct();
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
-
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            #include "alphaControls.H"
+
+            surfaceScalarField rhoPhi
+            (
+                IOobject
+                (
+                    "rhoPhi",
+                    runTime.timeName(),
+                    mesh
+                ),
+                mesh,
+                dimensionedScalar("0", dimMass/dimTime, 0)
+            );
+
+            if (pimple.firstIter() || alphaOuterCorrectors)
+            {
+                twoPhaseProperties->correct();
+
+                #include "alphaEqnSubCycle.H"
+                interface.correct();
+            }
+
             #include "UEqn.H"
 
             // --- Pressure corrector loop
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index 20afd1f72ffabb375ea915744af2b7c8d0da73f3..e41b0128603ef02eb028310bfa1576869ab7b54e 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -87,14 +87,14 @@ Foam::multiphaseMixture::multiphaseMixture
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("rho*phi", dimMass/dimTime, 0.0)
+        dimensionedScalar("rhoPhi", dimMass/dimTime, 0.0)
     ),
 
     alphas_
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
index 18b686f4f8afcfdcbf92f61f2e4f7c3a54d26893..31e1e61d5d3a6688952ee5b543c2efe96968c390 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
@@ -54,7 +54,7 @@
     (
         IOobject
         (
-            "rho*phi",
+            "rhoPhi",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
diff --git a/src/finiteVolume/cfdTools/general/include/alphaControls.H b/src/finiteVolume/cfdTools/general/include/alphaControls.H
index e7b3cf7794c14d1dc34d0c5d42ec2873ef73a333..e57869309aa3da416416158d413863ac6919b016 100644
--- a/src/finiteVolume/cfdTools/general/include/alphaControls.H
+++ b/src/finiteVolume/cfdTools/general/include/alphaControls.H
@@ -1,4 +1,12 @@
 const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+
 label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
+
 label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
-Switch MULESCorr(alphaControls.lookupOrDefault<Switch>("MULESCorr", false));
+
+bool MULESCorr(alphaControls.lookupOrDefault<Switch>("MULESCorr", false));
+
+bool alphaOuterCorrectors
+(
+    alphaControls.lookupOrDefault<Switch>("alphaOuterCorrectors", false)
+);
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
index 2c89cc4cae0a0eb08e80800de12ca2a12b3d1a43..8329580bffb30afa36a6cd77e9692b0703564deb 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
index 6cb9d3a21f652df26ddf6251473c7281de0e7611..43ac8d94f7de8416d2f1a30fba4412758524dc8f 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
@@ -129,6 +129,9 @@ public:
             //- Helper function to identify when to store the intial residuals
             inline bool storeInitialResiduals() const;
 
+            //- Helper function to identify first PIMPLE (outer) iteration
+            inline bool firstIter() const;
+
             //- Helper function to identify final PIMPLE (outer) iteration
             inline bool finalIter() const;
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
index 1d8917de89b34dfa38f65ef411f555de35ffe303..8ac86bdfa886f1c37929c6c2586931eea2b3951d 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,14 +66,20 @@ inline bool Foam::pimpleControl::correct()
 
 inline bool Foam::pimpleControl::storeInitialResiduals() const
 {
-    // start from second PIMPLE iteration
+    // Start from second PIMPLE iteration
     return (corr_ == 2) && (corrPISO_ == 0) && (corrNonOrtho_ == 0);
 }
 
 
+inline bool Foam::pimpleControl::firstIter() const
+{
+    return corr_ == 1;
+}
+
+
 inline bool Foam::pimpleControl::finalIter() const
 {
-    return converged_ || (corr_ == nCorrPIMPLE_);
+    return converged_ || (corr_ == corrPISO_);
 }
 
 
diff --git a/tutorials/multiphase/interFoam/laminar/damBreak/system/controlDict b/tutorials/multiphase/interFoam/laminar/damBreak/system/controlDict
index bae61318bbdfcedfb88a544d74c1a9e069d543bd..77ec4eba8758af0c0513ea7ca1fdcf130678df6d 100644
--- a/tutorials/multiphase/interFoam/laminar/damBreak/system/controlDict
+++ b/tutorials/multiphase/interFoam/laminar/damBreak/system/controlDict
@@ -47,8 +47,8 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.5;
-maxAlphaCo      0.5;
+maxCo           1;
+maxAlphaCo      1;
 
 maxDeltaT       1;
 
diff --git a/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution b/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution
index 36156e86145c247b6b99b8a34bd10a54fbbc1434..9e473a7865023654a2e5865a837a2989d6071492 100644
--- a/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution
+++ b/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSolution
@@ -19,16 +19,25 @@ solvers
 {
     alpha.water
     {
-        nAlphaCorr      1;
-        nAlphaSubCycles 2;
+        nAlphaCorr      2;
+        nAlphaSubCycles 1;
+        alphaOuterCorrectors yes;
         cAlpha          1;
+
+        MULESCorr       yes;
+        nLimiterIter    3;
+
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-8;
+        relTol          0;
     }
 
     pcorr
     {
         solver          PCG;
         preconditioner  DIC;
-        tolerance       1e-10;
+        tolerance       1e-5;
         relTol          0;
     }
 
@@ -43,14 +52,13 @@ solvers
     p_rghFinal
     {
         $p_rgh;
-        tolerance       1e-07;
         relTol          0;
     }
 
     U
     {
-        solver          PBiCG;
-        preconditioner  DILU;
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
         tolerance       1e-06;
         relTol          0;
     }
@@ -58,10 +66,22 @@ solvers
 
 PIMPLE
 {
-    momentumPredictor no;
-    nCorrectors     3;
+    momentumPredictor   no;
+    nOuterCorrectors    1;
+    nCorrectors         3;
     nNonOrthogonalCorrectors 0;
 }
 
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        ".*" 1;
+    }
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/polyMesh/boundary b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/polyMesh/boundary
index 8b9ef8e645041c16a8f156769cc5d6f6d74a0456..17c34d0093bec0469e5824177f10948c1f2b86ba 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/polyMesh/boundary
+++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/polyMesh/boundary
@@ -39,7 +39,7 @@ FoamFile
     bullet
     {
         type            wall;
-        nFaces          37752;
+        nFaces          37743;
         startFace       1133431;
     }
 )