diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index c7289b23f9..c4edc961d8 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 602affd4e0..12dc93795c 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 317a18120c..e9aa9fb9b9 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 0fa1a50ee2..56527781a2 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 b1cfcac9bd..8811a53df6 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 b0dd8ebef2..f6e57df771 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 81c65c1caf..5db9f16546 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 f0d1809848..3039171f77 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 4a673c6a3f..d56e300581 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();
-
-        mesh.update();
-
-        if (mesh.changing())
-        {
-            Info<< "Execution time for mesh.update() = "
-                << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
-                << " s" << endl;
-
-            gh = g & mesh.C();
-            ghf = g & mesh.Cf();
-        }
-
-        if (mesh.changing() && correctPhi)
-        {
-            // Calculate absolute flux from the mapped surface velocity
-            phi = mesh.Sf() & Uf;
-
-            #include "correctPhi.H"
-
-            // Make the flux relative to the mesh motion
-            fvc::makeRelative(phi, U);
-
-            interface.correct();
-        }
-
-        if (mesh.changing() && checkMeshCourantNo)
-        {
-            #include "meshCourantNo.H"
-        }
-
-        twoPhaseProperties.correct();
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
+        tmp<surfaceScalarField> tphiAlpha;
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    Info<< "Execution time for mesh.update() = "
+                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                        << " s" << endl;
+
+                    gh = g & mesh.C();
+                    ghf = g & mesh.Cf();
+                }
+
+                if (mesh.changing() && correctPhi)
+                {
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & Uf;
+
+                    #include "correctPhi.H"
+
+                    // Make the flux relative to the mesh motion
+                    fvc::makeRelative(phi, U);
+
+                    interface.correct();
+                }
+
+                if (mesh.changing() && checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
+            }
+
+            #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/interDyMFoam/readControls.H b/applications/solvers/multiphase/interFoam/interDyMFoam/readControls.H
index d4e332ff38..bf8b38f426 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 4ce08d27db..0ddfed5039 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 196c82194c..3fd5874d5b 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 12bc57d785..3321da57f9 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 3423a4654a..15b4439220 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 a4338b907f..077738ccb2 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 31fa619754..0bfd8ad56e 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();
-
-        mesh.update();
-
-        if (mesh.changing())
-        {
-            Info<< "Execution time for mesh.update() = "
-                << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
-                << " s" << endl;
-
-            gh = g & mesh.C();
-            ghf = g & mesh.Cf();
-        }
-
-        if (mesh.changing() && correctPhi)
-        {
-            // Calculate absolute flux from the mapped surface velocity
-            phi = mesh.Sf() & Uf;
-
-            #include "../interFoam/interDyMFoam/correctPhi.H"
-
-            // Make the flux relative to the mesh motion
-            fvc::makeRelative(phi, U);
-        }
-
-        if (mesh.changing() && checkMeshCourantNo)
-        {
-            #include "meshCourantNo.H"
-        }
-
-        twoPhaseProperties->correct();
-
-        #include "alphaEqnSubCycle.H"
-        interface.correct();
-
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    Info<< "Execution time for mesh.update() = "
+                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                        << " s" << endl;
+
+                    gh = g & mesh.C();
+                    ghf = g & mesh.Cf();
+                }
+
+                if (mesh.changing() && correctPhi)
+                {
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & Uf;
+
+                    #include "../interFoam/interDyMFoam/correctPhi.H"
+
+                    // Make the flux relative to the mesh motion
+                    fvc::makeRelative(phi, U);
+                }
+
+                if (mesh.changing() && checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
+            }
+
+            #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/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
index 28ef7f0a69..2a653dcf37 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 20afd1f72f..e41b012860 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 18b686f4f8..31e1e61d5d 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/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 65bc15d852..c349740ea7 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -1253,6 +1253,9 @@ int main(int argc, char *argv[])
     // Snap parameters
     const snapParameters snapParams(snapDict);
 
+    // Layer addition parameters
+    const layerParameters layerParams(layerDict, mesh.boundaryMesh());
+
 
     if (wantRefine)
     {
@@ -1346,9 +1349,6 @@ int main(int argc, char *argv[])
             globalToSlavePatch
         );
 
-        // Layer addition parameters
-        layerParameters layerParams(layerDict, mesh.boundaryMesh());
-
         // Use the maxLocalCells from the refinement parameters
         bool preBalance = returnReduce
         (
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 9db10ea85e..b850033a44 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -304,11 +304,13 @@ addLayersControls
     // size of the refined cell outside layer (true) or absolute sizes (false).
     relativeSizes true;
 
-    // Layer thickness specification. This can be specified in one of four ways
+    // Layer thickness specification. This can be specified in one of following
+    // ways:
     // - expansionRatio and finalLayerThickness (cell nearest internal mesh)
     // - expansionRatio and firstLayerThickness (cell on surface)
     // - overall thickness and firstLayerThickness
     // - overall thickness and finalLayerThickness
+    // - overall thickness and expansionRatio
 
         // Expansion factor for layer mesh
         expansionRatio 1.0;
diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
index 3ad100d4f2..aa7fa5554b 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
@@ -720,6 +720,20 @@ int main(int argc, char *argv[])
                 PtrList<faceSet> faceSets(fSetNames.size());
                 PtrList<pointSet> pointSets(pSetNames.size());
 
+                Info<< "Reconstructing sets:" << endl;
+                if (cSetNames.size())
+                {
+                    Info<< "    cellSets " << cSetNames.sortedToc() << endl;
+                }
+                if (fSetNames.size())
+                {
+                    Info<< "    faceSets " << fSetNames.sortedToc() << endl;
+                }
+                if (pSetNames.size())
+                {
+                    Info<< "    pointSets " << pSetNames.sortedToc() << endl;
+                }
+
                 // Load sets
                 forAll(procMeshes.meshes(), procI)
                 {
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
index 7c94b9b138..96004c3e26 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
@@ -542,8 +542,11 @@ public:
             //- Compact maps. Gets per field a bool whether it is used (locally)
             //  and works out itself what this side and sender side can remove
             //  from maps.
-            void compact(const boolList& elemIsUsed, const int tag);
-
+            void compact
+            (
+                const boolList& elemIsUsed,
+                const int tag = UPstream::msgType()
+            );
 
             //- Distribute data. Note:schedule only used for Pstream::scheduled
             //  for now, all others just use send-to-all, receive-from-all.
diff --git a/src/finiteVolume/cfdTools/general/include/alphaControls.H b/src/finiteVolume/cfdTools/general/include/alphaControls.H
index e7b3cf7794..e57869309a 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 2c89cc4cae..8329580bff 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 6cb9d3a21f..510e308188 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.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
@@ -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 1d8917de89..8ac86bdfa8 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/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 2c6add9fc9..7bb72786b8 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -52,6 +52,7 @@ Description
 #include "fixedValuePointPatchFields.H"
 #include "calculatedPointPatchFields.H"
 #include "cyclicSlipPointPatchFields.H"
+#include "fixedValueFvPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -110,6 +111,31 @@ void Foam::autoLayerDriver::dumpDisplacement
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::autoLayerDriver::avgPointData
+(
+    const indirectPrimitivePatch& pp,
+    const scalarField& pointFld
+)
+{
+    tmp<scalarField> tfaceFld(new scalarField(pp.size(), 0.0));
+    scalarField& faceFld = tfaceFld();
+
+    forAll(pp.localFaces(), faceI)
+    {
+        const face& f = pp.localFaces()[faceI];
+        if (f.size())
+        {
+            forAll(f, fp)
+            {
+                faceFld[faceI] += pointFld[f[fp]];
+            }
+            faceFld[faceI] /= f.size();
+        }
+    }
+    return tfaceFld;
+}
+
+
 // Check that primitivePatch is not multiply connected. Collect non-manifold
 // points in pointSet.
 void Foam::autoLayerDriver::checkManifold
@@ -2391,19 +2417,21 @@ Foam::label Foam::autoLayerDriver::countExtrusion
 }
 
 
-// Collect layer faces and layer cells into bools for ease of handling
+// Collect layer faces and layer cells into mesh fields for ease of handling
 void Foam::autoLayerDriver::getLayerCellsFaces
 (
     const polyMesh& mesh,
     const addPatchCellLayer& addLayer,
-    boolList& flaggedCells,
-    boolList& flaggedFaces
+    const scalarField& oldRealThickness,
+
+    labelList& cellNLayers,
+    scalarField& faceRealThickness
 )
 {
-    flaggedCells.setSize(mesh.nCells());
-    flaggedCells = false;
-    flaggedFaces.setSize(mesh.nFaces());
-    flaggedFaces = false;
+    cellNLayers.setSize(mesh.nCells());
+    cellNLayers = 0;
+    faceRealThickness.setSize(mesh.nFaces());
+    faceRealThickness = 0;
 
     // Mark all faces in the layer
     const labelListList& layerFaces = addLayer.layerFaces();
@@ -2415,27 +2443,195 @@ void Foam::autoLayerDriver::getLayerCellsFaces
     {
         const labelList& added = addedCells[oldPatchFaceI];
 
-        forAll(added, i)
+        const labelList& layer = layerFaces[oldPatchFaceI];
+
+        if (layer.size())
         {
-            flaggedCells[added[i]] = true;
+            forAll(added, i)
+            {
+                cellNLayers[added[i]] = layer.size()-1;
+            }
         }
     }
 
     forAll(layerFaces, oldPatchFaceI)
     {
         const labelList& layer = layerFaces[oldPatchFaceI];
+        const scalar realThickness = oldRealThickness[oldPatchFaceI];
 
         if (layer.size())
         {
-            for (label i = 1; i < layer.size()-1; i++)
+            // Layer contains both original boundary face and new boundary
+            // face so is nLayers+1
+            forAll(layer, i)
             {
-                flaggedFaces[layer[i]] = true;
+                faceRealThickness[layer[i]] = realThickness;
             }
         }
     }
 }
 
 
+bool Foam::autoLayerDriver::writeLayerData
+(
+    const fvMesh& mesh,
+    const labelList& patchIDs,
+    const labelList& cellNLayers,
+    const scalarField& faceWantedThickness,
+    const scalarField& faceRealThickness
+) const
+{
+    bool allOk = true;
+
+    {
+        label nAdded = 0;
+        forAll(cellNLayers, cellI)
+        {
+            if (cellNLayers[cellI] > 0)
+            {
+                nAdded++;
+            }
+        }
+        cellSet addedCellSet(mesh, "addedCells", nAdded);
+        forAll(cellNLayers, cellI)
+        {
+            if (cellNLayers[cellI] > 0)
+            {
+                addedCellSet.insert(cellI);
+            }
+        }
+        addedCellSet.instance() = meshRefiner_.timeName();
+        Info<< "Writing "
+            << returnReduce(addedCellSet.size(), sumOp<label>())
+            << " added cells to cellSet "
+            << addedCellSet.name() << endl;
+        bool ok = addedCellSet.write();
+        allOk = allOk & ok;
+    }
+    {
+        label nAdded = 0;
+        forAll(faceRealThickness, faceI)
+        {
+            if (faceRealThickness[faceI] > 0)
+            {
+                nAdded++;
+            }
+        }
+
+        faceSet layerFacesSet(mesh, "layerFaces", nAdded);
+        forAll(faceRealThickness, faceI)
+        {
+            if (faceRealThickness[faceI] > 0)
+            {
+                layerFacesSet.insert(faceI);
+            }
+        }
+        layerFacesSet.instance() = meshRefiner_.timeName();
+        Info<< "Writing "
+            << returnReduce(layerFacesSet.size(), sumOp<label>())
+            << " faces inside added layer to faceSet "
+            << layerFacesSet.name() << endl;
+        bool ok = layerFacesSet.write();
+        allOk = allOk & ok;
+    }
+    {
+        volScalarField fld
+        (
+            IOobject
+            (
+                "nSurfaceLayers",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0),
+            fixedValueFvPatchScalarField::typeName
+        );
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
+            const polyPatch& pp = pbm[patchI];
+            const labelList& faceCells = pp.faceCells();
+            scalarField pfld(faceCells.size());
+            forAll(faceCells, i)
+            {
+                pfld[i] = cellNLayers[faceCells[i]];
+            }
+            fld.boundaryField()[patchI] == pfld;
+        }
+        Info<< "Writing volScalarField " << fld.name()
+            << " with actual number of layers" << endl;
+        bool ok = fld.write();
+        allOk = allOk & ok;
+    }
+    {
+        volScalarField fld
+        (
+            IOobject
+            (
+                "wantedThickness",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0),
+            fixedValueFvPatchScalarField::typeName
+        );
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
+            fld.boundaryField()[patchI] == pbm[patchI].patchSlice
+            (
+                faceWantedThickness
+            );
+        }
+        Info<< "Writing volScalarField " << fld.name()
+            << " with wanted thickness" << endl;
+        bool ok = fld.write();
+        allOk = allOk & ok;
+    }
+    {
+        volScalarField fld
+        (
+            IOobject
+            (
+                "thickness",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0),
+            fixedValueFvPatchScalarField::typeName
+        );
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
+            fld.boundaryField()[patchI] == pbm[patchI].patchSlice
+            (
+                faceRealThickness
+            );
+        }
+        Info<< "Writing volScalarField " << fld.name()
+            << " with layer thickness" << endl;
+        bool ok = fld.write();
+        allOk = allOk & ok;
+    }
+    return allOk;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::autoLayerDriver::autoLayerDriver
@@ -2618,7 +2814,7 @@ void Foam::autoLayerDriver::addLayers
     labelList patchNLayers(pp().nPoints(), 0);
 
     // Ideal number of cells added
-    label nIdealAddedCells = 0;
+    label nIdealTotAddedCells = 0;
 
     // Whether to add edge for all pp.localPoints.
     List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
@@ -2637,7 +2833,7 @@ void Foam::autoLayerDriver::addLayers
             patchDisp,
             patchNLayers,
             extrudeStatus,
-            nIdealAddedCells
+            nIdealTotAddedCells
         );
 
         // Precalculate mesh edge labels for patch edges
@@ -2835,11 +3031,20 @@ void Foam::autoLayerDriver::addLayers
     // Saved old points
     pointField oldPoints(mesh.points());
 
-    // Last set of topology changes. (changing mesh clears out polyTopoChange)
+    // Current set of topology changes. (changing mesh clears out
+    // polyTopoChange)
     polyTopoChange savedMeshMod(mesh.boundaryMesh().size());
+    // Per cell 0 or number of layers in the cell column it is part of
+    labelList cellNLayers;
+    // Per face actual overall layer thickness
+    scalarField faceRealThickness;
+    // Per face wanted overall layer thickness
+    scalarField faceWantedThickness(mesh.nFaces(), 0.0);
+    {
+        UIndirectList<scalar>(faceWantedThickness, pp().addressing()) =
+            avgPointData(pp, thickness);
+    }
 
-    boolList flaggedCells;
-    boolList flaggedFaces;
 
     for (label iteration = 0; iteration < layerParams.nLayerIter(); iteration++)
     {
@@ -3088,35 +3293,52 @@ void Foam::autoLayerDriver::addLayers
         (
             newMesh,
             addLayer,
-            flaggedCells,
-            flaggedFaces
+            avgPointData(pp, mag(patchDisp))(), // current thickness
+
+            cellNLayers,
+            faceRealThickness
         );
 
 
+        // Count number of added cells
+        label nAddedCells = 0;
+        forAll(cellNLayers, cellI)
+        {
+            if (cellNLayers[cellI] > 0)
+            {
+                nAddedCells++;
+            }
+        }
+
+
         if (debug&meshRefinement::MESH)
         {
             Info<< "Writing layer mesh to time " << meshRefiner_.timeName()
                 << endl;
             newMesh.write();
-            cellSet addedCellSet
-            (
-                newMesh,
-                "addedCells",
-                findIndices(flaggedCells, true)
-            );
+
+            cellSet addedCellSet(newMesh," addedCells", nAddedCells);
+            forAll(cellNLayers, cellI)
+            {
+                if (cellNLayers[cellI] > 0)
+                {
+                    addedCellSet.insert(cellI);
+                }
+            }
             addedCellSet.instance() = meshRefiner_.timeName();
             Info<< "Writing "
                 << returnReduce(addedCellSet.size(), sumOp<label>())
-                << " added cells to cellSet "
-                << addedCellSet.name() << endl;
+                << " added cells to cellSet " << addedCellSet.name() << endl;
             addedCellSet.write();
 
-            faceSet layerFacesSet
-            (
-                newMesh,
-                "layerFaces",
-                findIndices(flaggedCells, true)
-            );
+            faceSet layerFacesSet(newMesh, "layerFaces", newMesh.nFaces()/100);
+            forAll(faceRealThickness, faceI)
+            {
+                if (faceRealThickness[faceI] > 0)
+                {
+                    layerFacesSet.insert(faceI);
+                }
+            }
             layerFacesSet.instance() = meshRefiner_.timeName();
             Info<< "Writing "
                 << returnReduce(layerFacesSet.size(), sumOp<label>())
@@ -3140,26 +3362,17 @@ void Foam::autoLayerDriver::addLayers
             extrudeStatus
         );
 
-        label nExtruded = countExtrusion(pp, extrudeStatus);
+        label nTotExtruded = countExtrusion(pp, extrudeStatus);
         label nTotFaces = returnReduce(pp().size(), sumOp<label>());
-        label nAddedCells = 0;
-        {
-            forAll(flaggedCells, cellI)
-            {
-                if (flaggedCells[cellI])
-                {
-                    nAddedCells++;
-                }
-            }
-            reduce(nAddedCells, sumOp<label>());
-        }
-        Info<< "Extruding " << nExtruded
+        label nTotAddedCells = returnReduce(nAddedCells, sumOp<label>());
+
+        Info<< "Extruding " << nTotExtruded
             << " out of " << nTotFaces
-            << " faces (" << 100.0*nExtruded/nTotFaces << "%)."
+            << " faces (" << 100.0*nTotExtruded/nTotFaces << "%)."
             << " Removed extrusion at " << nTotChanged << " faces."
             << endl
-            << "Added " << nAddedCells << " out of " << nIdealAddedCells
-            << " cells (" << 100.0*nAddedCells/nIdealAddedCells << "%)."
+            << "Added " << nTotAddedCells << " out of " << nIdealTotAddedCells
+            << " cells (" << 100.0*nTotAddedCells/nIdealTotAddedCells << "%)."
             << endl;
 
         if (nTotChanged == 0)
@@ -3217,6 +3430,8 @@ void Foam::autoLayerDriver::addLayers
 
     meshRefiner_.updateMesh(map, labelList(0));
 
+    // Update numbering of faceWantedThickness
+    meshRefinement::updateList(map().faceMap(), 0.0, faceWantedThickness);
 
     // Update numbering on baffles
     forAll(baffles, i)
@@ -3237,8 +3452,9 @@ void Foam::autoLayerDriver::addLayers
 
         autoPtr<mapPolyMesh> map = meshRefiner_.mergeBaffles(baffles);
 
-        inplaceReorder(map().reverseCellMap(), flaggedCells);
-        inplaceReorder(map().reverseFaceMap(), flaggedFaces);
+        inplaceReorder(map().reverseCellMap(), cellNLayers);
+        inplaceReorder(map().reverseFaceMap(), faceWantedThickness);
+        inplaceReorder(map().reverseFaceMap(), faceRealThickness);
 
         Info<< "Converted baffles in = "
             << meshRefiner_.mesh().time().cpuTimeIncrement()
@@ -3272,29 +3488,23 @@ void Foam::autoLayerDriver::addLayers
         );
 
         // Re-distribute flag of layer faces and cells
-        map().distributeCellData(flaggedCells);
-        map().distributeFaceData(flaggedFaces);
+        map().distributeCellData(cellNLayers);
+        map().distributeFaceData(faceWantedThickness);
+        map().distributeFaceData(faceRealThickness);
     }
 
 
-    // Write mesh
-    // ~~~~~~~~~~
+    // Write mesh data
+    // ~~~~~~~~~~~~~~~
 
-    cellSet addedCellSet(mesh, "addedCells", findIndices(flaggedCells, true));
-    addedCellSet.instance() = meshRefiner_.timeName();
-    Info<< "Writing "
-        << returnReduce(addedCellSet.size(), sumOp<label>())
-        << " added cells to cellSet "
-        << addedCellSet.name() << endl;
-    addedCellSet.write();
-
-    faceSet layerFacesSet(mesh, "layerFaces", findIndices(flaggedFaces, true));
-    layerFacesSet.instance() = meshRefiner_.timeName();
-    Info<< "Writing "
-        << returnReduce(layerFacesSet.size(), sumOp<label>())
-        << " faces inside added layer to faceSet "
-        << layerFacesSet.name() << endl;
-    layerFacesSet.write();
+    writeLayerData
+    (
+        mesh,
+        patchIDs,
+        cellNLayers,
+        faceWantedThickness,
+        faceRealThickness
+    );
 }
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index df7ea815a0..4534ebe00d 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -120,6 +120,13 @@ class autoLayerDriver
                 const List<extrudeMode>&
             );
 
+            //- Average point wise data to face wise
+            static tmp<scalarField> avgPointData
+            (
+                const indirectPrimitivePatch&,
+                const scalarField& pointFld
+            );
+
             //- Check that primitivePatch is not multiply connected.
             //  Collect non-manifold points in pointSet.
             static void checkManifold
@@ -367,10 +374,23 @@ class autoLayerDriver
                 (
                     const polyMesh&,
                     const addPatchCellLayer&,
-                    boolList&,
-                    boolList&
+                    const scalarField& oldRealThickness,
+
+                    labelList& cellStatus,
+                    scalarField& faceRealThickness
                 );
 
+                //- Write cellSet,faceSet for layers
+                bool writeLayerData
+                (
+                    const fvMesh& mesh,
+                    const labelList& patchIDs,
+                    const labelList& cellNLayers,
+                    const scalarField& faceWantedThickness,
+                    const scalarField& faceRealThickness
+                ) const;
+
+
             // Mesh shrinking (to create space for layers)
 
                 //- Average field (over all subset of mesh points) by
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
index 8444b6632c..fd64485b62 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
@@ -235,6 +235,12 @@ Foam::layerParameters::layerParameters
         Info<< "Layer thickness specified as final layer and expansion ratio."
             << endl;
     }
+    else if (haveTotal && haveExp)
+    {
+        layerSpec_ = TOTAL_AND_EXPANSION;
+        Info<< "Layer thickness specified as overall thickness"
+            << " and expansion ratio." << endl;
+    }
 
 
     if (layerSpec_ == ILLEGAL || nSpec != 2)
@@ -253,7 +259,9 @@ Foam::layerParameters::layerParameters
             << "    final layer thickness ('finalLayerThickness')"
             << " and expansion ratio ('expansionRatio') or" << nl
             << "    final layer thickness ('finalLayerThickness')"
-            << " and overall thickness ('thickness')"
+            << " and overall thickness ('thickness') or" << nl
+            << "    overall thickness ('thickness')"
+            << " and expansion ratio ('expansionRatio'"
             << exit(FatalIOError);
     }
 
@@ -354,6 +362,19 @@ Foam::layerParameters::layerParameters
                             );
                         break;
 
+                        case TOTAL_AND_EXPANSION:
+                            layerDict.readIfPresent
+                            (
+                                "thickness",
+                                thickness_[patchI]
+                            );
+                            layerDict.readIfPresent
+                            (
+                                "expansionRatio",
+                                expansionRatio_[patchI]
+                            );
+                        break;
+
                         default:
                             FatalIOErrorIn
                             (
@@ -390,6 +411,7 @@ Foam::scalar Foam::layerParameters::layerThickness
     {
         case FIRST_AND_TOTAL:
         case FINAL_AND_TOTAL:
+        case TOTAL_AND_EXPANSION:
         {
             return totalThickness;
         }
@@ -450,6 +472,7 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
     {
         case FIRST_AND_EXPANSION:
         case FINAL_AND_EXPANSION:
+        case TOTAL_AND_EXPANSION:
         {
             return expansionRatio;
         }
@@ -524,6 +547,18 @@ Foam::scalar Foam::layerParameters::firstLayerThickness
         }
         break;
 
+        case TOTAL_AND_EXPANSION:
+        {
+            scalar r = finalLayerThicknessRatio
+            (
+                nLayers,
+                expansionRatio
+            );
+            scalar finalThickness = r*totalThickness;
+            return finalThickness/pow(expansionRatio, nLayers-1);
+        }
+        break;
+
         default:
         {
             FatalErrorIn("layerParameters::layerThickness(..)")
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
index 4005dd5df8..5eaa529647 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
@@ -64,13 +64,15 @@ public:
         //  - first and expansion ratio specified
         //  - final and total thickness specified
         //  - final and expansion ratio specified
+        //  - total thickness and expansion ratio specified
         enum layerSpecification
         {
             ILLEGAL,
             FIRST_AND_TOTAL,
             FIRST_AND_EXPANSION,
             FINAL_AND_TOTAL,
-            FINAL_AND_EXPANSION
+            FINAL_AND_EXPANSION,
+            TOTAL_AND_EXPANSION
         };
 
 
diff --git a/src/parallel/reconstruct/reconstruct/processorMeshes.C b/src/parallel/reconstruct/reconstruct/processorMeshes.C
index b3729c8b7e..875e96c38e 100644
--- a/src/parallel/reconstruct/reconstruct/processorMeshes.C
+++ b/src/parallel/reconstruct/reconstruct/processorMeshes.C
@@ -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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,15 @@ License
 
 void Foam::processorMeshes::read()
 {
+    forAll(databases_, procI)
+    {
+        meshes_.set(procI, NULL);
+        pointProcAddressing_.set(procI, NULL);
+        faceProcAddressing_.set(procI, NULL);
+        cellProcAddressing_.set(procI, NULL);
+        boundaryProcAddressing_.set(procI, NULL);
+    }
+
     forAll(databases_, procI)
     {
         meshes_.set
diff --git a/tutorials/multiphase/interFoam/laminar/damBreak/system/controlDict b/tutorials/multiphase/interFoam/laminar/damBreak/system/controlDict
index bae61318bb..77ec4eba87 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 36156e8614..9e473a7865 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 8b9ef8e645..17c34d0093 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;
     }
 )