diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index 3f4f485cf6..eb263b8a50 100644
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@@ -4,6 +4,9 @@ volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
+phi.boundaryField() =
+    fvc::interpolate(rho.boundaryField())
+   *(fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField());
 
 surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/UEqn.H b/applications/solvers/multiphase/cavitatingFoam/UEqn.H
index 9a5761b59f..1ca0b9f0a5 100644
--- a/applications/solvers/multiphase/cavitatingFoam/UEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/UEqn.H
@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
-      - fvm::laplacian(muEff, U)
-    //- (fvc::grad(U) & fvc::grad(muf))
-      - fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
index 257f6d48b5..7cc250a66a 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
index 67717b0738..709d30c398 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
@@ -22,27 +22,29 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
+    surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
+
+    if (g0.value() > 0.0)
+    {
+        ppDrag -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
+    }
+
+    if (kineticTheory.on())
+    {
+        ppDrag -=
+            fvc::interpolate(1.0/rho1)*rAlphaAU1f
+           *fvc::snGrad(kineticTheory.pa())*mesh.magSf();
+    }
+
     surfaceScalarField phiHbyA1
     (
         "phiHbyA1",
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
       + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
+      + ppDrag
       + rAlphaAU1f*(g & mesh.Sf())
     );
-
-    if (g0.value() > 0.0)
-    {
-        phiHbyA1 -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
-    }
-
-    if (kineticTheory.on())
-    {
-        phiHbyA1 -=
-            fvc::interpolate((1.0/rho1)*rAU1)
-           *fvc::snGrad(kineticTheory.pa())*mesh.magSf();
-    }
-
     mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2
@@ -135,7 +137,8 @@
             U1 = HbyA1
                + fvc::reconstruct
                  (
-                     rAlphaAU1f
+                     ppDrag
+                   + rAlphaAU1f
                     *(
                          (g & mesh.Sf())
                        + mSfGradp/fvc::interpolate(rho1)
diff --git a/applications/solvers/multiphase/interFoam/UEqn.H b/applications/solvers/multiphase/interFoam/UEqn.H
index 257f6d48b5..7cc250a66a 100644
--- a/applications/solvers/multiphase/interFoam/UEqn.H
+++ b/applications/solvers/multiphase/interFoam/UEqn.H
@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index d942f2c501..42721c0018 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@@ -4,6 +4,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
index b8fe82ff14..062c5523c9 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
@@ -1,18 +1,9 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties->muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
       - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
-      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
index 43ee1e3eb9..fbcba7db72 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        mixture.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(mixture.rhoPhi(), U)
-      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/settlingFoam/UEqn.H b/applications/solvers/multiphase/settlingFoam/UEqn.H
index d6232da309..ff6323bca2 100644
--- a/applications/solvers/multiphase/settlingFoam/UEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/UEqn.H
@@ -11,7 +11,6 @@
         )
       - fvm::laplacian(muEff, U, "laplacian(muEff,U)")
       - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*dev2(T(fvc::grad(U))))
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H
index 3dbb374774..7afd323dd0 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H
@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+      + turbulence->divDevRhoReff(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatchDict b/applications/utilities/mesh/manipulation/createPatch/createPatchDict
index 35e0b86926..f731f78df9 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatchDict
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatchDict
@@ -20,6 +20,7 @@ FoamFile
 // - always: order faces on coupled patches such that they are opposite. This
 //   is done for all coupled faces, not just for any patches created.
 // - optional: synchronise points on coupled patches.
+// - always: remove zero-sized (non-coupled) patches (that were not added)
 
 // 1. Create cyclic:
 // - specify where the faces should come from
diff --git a/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C b/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C
index 2518dc0450..0969b1e560 100644
--- a/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C
+++ b/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,6 @@ Usage
 
 #include "argList.H"
 #include "Time.H"
-#include "timeSelector.H"
 #include "fvMesh.H"
 #include "unitConversion.H"
 #include "polyTopoChange.H"
@@ -354,7 +353,6 @@ int main(int argc, char *argv[])
 {
 #   include "addOverwriteOption.H"
     argList::noParallel();
-    timeSelector::addOptions(true, false);
 
     argList::validArgs.append("featureAngle [0-180]");
     argList::addBoolOption
@@ -376,9 +374,6 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
-
-    instantList timeDirs = timeSelector::select0(runTime, args);
-
 #   include "createMesh.H"
 
     const word oldInstance = mesh.pointsInstance();
diff --git a/applications/utilities/preProcessing/setFields/setFields.C b/applications/utilities/preProcessing/setFields/setFields.C
index 41c09b6053..ebc43d84f2 100644
--- a/applications/utilities/preProcessing/setFields/setFields.C
+++ b/applications/utilities/preProcessing/setFields/setFields.C
@@ -27,7 +27,6 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 #include "fvMesh.H"
 #include "topoSetSource.H"
@@ -360,14 +359,8 @@ public:
 
 int main(int argc, char *argv[])
 {
-    timeSelector::addOptions();
-
 #   include "setRootCase.H"
 #   include "createTime.H"
-
-    // Get times list
-    instantList timeDirs = timeSelector::select0(runTime, args);
-
 #   include "createMesh.H"
 
     Info<< "Reading setFieldsDict\n" << endl;
diff --git a/applications/utilities/surface/surfaceMeshConvert/surfaceMeshConvert.C b/applications/utilities/surface/surfaceMeshConvert/surfaceMeshConvert.C
index 357b2bd198..305786f963 100644
--- a/applications/utilities/surface/surfaceMeshConvert/surfaceMeshConvert.C
+++ b/applications/utilities/surface/surfaceMeshConvert/surfaceMeshConvert.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,7 +58,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "MeshedSurfaces.H"
diff --git a/applications/utilities/surface/surfaceMeshConvertTesting/surfaceMeshConvertTesting.C b/applications/utilities/surface/surfaceMeshConvertTesting/surfaceMeshConvertTesting.C
index 92f0ffbd3d..cbc5059f10 100644
--- a/applications/utilities/surface/surfaceMeshConvertTesting/surfaceMeshConvertTesting.C
+++ b/applications/utilities/surface/surfaceMeshConvertTesting/surfaceMeshConvertTesting.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,7 +52,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 #include "polyMesh.H"
 #include "triSurface.H"
diff --git a/applications/utilities/surface/surfaceMeshExport/surfaceMeshExport.C b/applications/utilities/surface/surfaceMeshExport/surfaceMeshExport.C
index c9b7d62113..7c9b3be904 100644
--- a/applications/utilities/surface/surfaceMeshExport/surfaceMeshExport.C
+++ b/applications/utilities/surface/surfaceMeshExport/surfaceMeshExport.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "MeshedSurfaces.H"
diff --git a/applications/utilities/surface/surfaceMeshImport/surfaceMeshImport.C b/applications/utilities/surface/surfaceMeshImport/surfaceMeshImport.C
index e9f7f63710..3a14ac63c4 100644
--- a/applications/utilities/surface/surfaceMeshImport/surfaceMeshImport.C
+++ b/applications/utilities/surface/surfaceMeshImport/surfaceMeshImport.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "MeshedSurfaces.H"
diff --git a/applications/utilities/surface/surfaceMeshInfo/surfaceMeshInfo.C b/applications/utilities/surface/surfaceMeshInfo/surfaceMeshInfo.C
index 08e6a99b38..80638bf4db 100644
--- a/applications/utilities/surface/surfaceMeshInfo/surfaceMeshInfo.C
+++ b/applications/utilities/surface/surfaceMeshInfo/surfaceMeshInfo.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "UnsortedMeshedSurfaces.H"
diff --git a/src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionaryIO.C b/src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionaryIO.C
index 8585392e27..881250f3eb 100644
--- a/src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionaryIO.C
+++ b/src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionaryIO.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,7 +102,7 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
                 myComm.above(),
                 0,
                 Pstream::msgType(),
-                IOstream::ASCII
+                IOstream::BINARY
             );
             IOdictionary::readData(fromAbove);
         }
@@ -121,7 +121,7 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
                 myComm.below()[belowI],
                 0,
                 Pstream::msgType(),
-                IOstream::ASCII
+                IOstream::BINARY
             );
             IOdictionary::writeData(toBelow);
         }
diff --git a/src/OpenFOAM/db/Time/timeSelector.C b/src/OpenFOAM/db/Time/timeSelector.C
index 865f3ee911..60b589a347 100644
--- a/src/OpenFOAM/db/Time/timeSelector.C
+++ b/src/OpenFOAM/db/Time/timeSelector.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -147,7 +147,7 @@ void Foam::timeSelector::addOptions
     (
         "time",
         "ranges",
-        "comma-separated time ranges - eg, ':10,20,40-70,1000:'"
+        "comma-separated time ranges - eg, ':10,20,40:70,1000:'"
     );
 }
 
diff --git a/src/OpenFOAM/db/dictionary/dictionaryIO.C b/src/OpenFOAM/db/dictionary/dictionaryIO.C
index 03ce2b76a7..aa1321d914 100644
--- a/src/OpenFOAM/db/dictionary/dictionaryIO.C
+++ b/src/OpenFOAM/db/dictionary/dictionaryIO.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,6 +80,12 @@ Foam::autoPtr<Foam::dictionary> Foam::dictionary::New(Istream& is)
 
 bool Foam::dictionary::read(Istream& is, const bool keepHeader)
 {
+    // Check for empty dictionary
+    if (is.eof())
+    {
+        return true;
+    }
+
     if (!is.good())
     {
         FatalIOErrorIn("dictionary::read(Istream&, bool)", is)
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
index 3a3c234821..1f9c51ac9f 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
@@ -323,8 +323,6 @@ class globalMeshData
             // its own master. Maybe store as well?
 
             void calcGlobalCoPointSlaves() const;
-            const labelListList& globalCoPointSlaves() const;
-            const mapDistribute& globalCoPointSlavesMap() const;
 
 
         //- Disallow default bitwise copy construct
@@ -547,6 +545,11 @@ public:
                 //- Is my edge same orientation as master edge
                 const PackedBoolList& globalEdgeOrientation() const;
 
+            // Collocated point to collocated point
+
+                const labelListList& globalCoPointSlaves() const;
+                const mapDistribute& globalCoPointSlavesMap() const;
+
             // Coupled point to boundary faces. These are uncoupled boundary
             // faces only but include empty patches.
 
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C
index 7ecd1b84a1..b4333fc3a7 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -339,6 +339,31 @@ bool Foam::globalPoints::storeInitialInfo
 }
 
 
+void Foam::globalPoints::printProcPoint
+(
+    const labelList& patchToMeshPoint,
+    const labelPair& pointInfo
+) const
+{
+    label procI = globalIndexAndTransform::processor(pointInfo);
+    label index = globalIndexAndTransform::index(pointInfo);
+    label trafoI = globalIndexAndTransform::transformIndex(pointInfo);
+
+    Pout<< "    proc:" << procI;
+    Pout<< " localpoint:";
+    Pout<< index;
+    Pout<< " through transform:"
+        << trafoI << " bits:"
+        << globalTransforms_.decodeTransformIndex(trafoI);
+
+    if (procI == Pstream::myProcNo())
+    {
+        label meshPointI = localToMeshPoint(patchToMeshPoint, index);
+        Pout<< " at:" <<  mesh_.points()[meshPointI];
+    }
+}
+
+
 void Foam::globalPoints::printProcPoints
 (
     const labelList& patchToMeshPoint,
@@ -347,23 +372,7 @@ void Foam::globalPoints::printProcPoints
 {
     forAll(pointInfo, i)
     {
-        label procI = globalIndexAndTransform::processor(pointInfo[i]);
-        label index = globalIndexAndTransform::index(pointInfo[i]);
-        label trafoI = globalIndexAndTransform::transformIndex(pointInfo[i]);
-
-        Pout<< "    proc:" << procI;
-        Pout<< " localpoint:";
-        Pout<< index;
-        Pout<< " through transform:"
-            << trafoI << " bits:"
-            << globalTransforms_.decodeTransformIndex(trafoI);
-
-        if (procI == Pstream::myProcNo())
-        {
-            label meshPointI = localToMeshPoint(patchToMeshPoint, index);
-            Pout<< " at:" <<  mesh_.points()[meshPointI];
-        }
-
+        printProcPoint(patchToMeshPoint, pointInfo[i]);
         Pout<< endl;
     }
 }
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.H
index daef012fe7..de7a6b88d0 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -206,6 +206,12 @@ class globalPoints
         );
 
         //- Debug printing
+        void printProcPoint
+        (
+            const labelList& patchToMeshPoint,
+            const labelPair& pointInfo
+        ) const;
+
         void printProcPoints
         (
             const labelList& patchToMeshPoint,
diff --git a/src/Pstream/mpi/UIPread.C b/src/Pstream/mpi/UIPread.C
index bcc6f4c412..19d6b16112 100644
--- a/src/Pstream/mpi/UIPread.C
+++ b/src/Pstream/mpi/UIPread.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,12 +107,7 @@ Foam::UIPstream::UIPstream
 
         if (!messageSize_)
         {
-            FatalErrorIn
-            (
-                "UIPstream::UIPstream(const commsTypes, const int, "
-                "DynamicList<char>&, streamFormat, versionNumber)"
-            )   << "read failed"
-                << Foam::abort(FatalError);
+            setEof();
         }
     }
 }
@@ -199,11 +194,7 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
 
         if (!messageSize_)
         {
-            FatalErrorIn
-            (
-                "UIPstream::UIPstream(const int, PstreamBuffers&)"
-            )   << "read failed"
-                << Foam::abort(FatalError);
+            setEof();
         }
     }
 }
diff --git a/src/finiteVolume/fields/fvPatchFields/basic/calculated/calculatedFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/basic/calculated/calculatedFvPatchField.C
index 1f179134fe..79b8c3750e 100644
--- a/src/finiteVolume/fields/fvPatchFields/basic/calculated/calculatedFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/basic/calculated/calculatedFvPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -98,20 +98,19 @@ calculatedFvPatchField<Type>::calculatedFvPatchField
 
 
 template<class Type>
-template<class Type2>
 tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
 (
-    const fvPatchField<Type2>& pf
+    const fvPatch& p
 )
 {
     typename patchConstructorTable::iterator patchTypeCstrIter =
-        patchConstructorTablePtr_->find(pf.patch().type());
+        patchConstructorTablePtr_->find(p.type());
 
     if (patchTypeCstrIter != patchConstructorTablePtr_->end())
     {
         return patchTypeCstrIter()
         (
-            pf.patch(),
+            p,
             DimensionedField<Type, volMesh>::null()
         );
     }
@@ -121,7 +120,7 @@ tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
         (
             new calculatedFvPatchField<Type>
             (
-                pf.patch(),
+                p,
                 DimensionedField<Type, volMesh>::null()
             )
         );
@@ -129,6 +128,17 @@ tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
 }
 
 
+template<class Type>
+template<class Type2>
+tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
+(
+    const fvPatchField<Type2>& pf
+)
+{
+    return NewCalculatedType(pf.patch());
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
diff --git a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H
index 4dc1041ffc..ecace2a051 100644
--- a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H
@@ -257,6 +257,13 @@ public:
             const dictionary&
         );
 
+        //- Return a pointer to a new calculatedFvPatchField created on
+        //  freestore without setting patchField values
+        static tmp<fvPatchField<Type> > NewCalculatedType
+        (
+            const fvPatch&
+        );
+
         //- Return a pointer to a new calculatedFvPatchField created on
         //  freestore without setting patchField values
         template<class Type2>
diff --git a/src/finiteVolume/fields/fvsPatchFields/basic/calculated/calculatedFvsPatchField.C b/src/finiteVolume/fields/fvsPatchFields/basic/calculated/calculatedFvsPatchField.C
index a22514c857..efe80d2f73 100644
--- a/src/finiteVolume/fields/fvsPatchFields/basic/calculated/calculatedFvsPatchField.C
+++ b/src/finiteVolume/fields/fvsPatchFields/basic/calculated/calculatedFvsPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -97,21 +97,20 @@ calculatedFvsPatchField<Type>::calculatedFvsPatchField
 
 
 template<class Type>
-template<class Type2>
 tmp<fvsPatchField<Type> > fvsPatchField<Type>::NewCalculatedType
 (
-    const fvsPatchField<Type2>& pf
+    const fvPatch& p
 )
 {
     typename patchConstructorTable::iterator patchTypeCstrIter =
-        patchConstructorTablePtr_->find(pf.patch().type());
+        patchConstructorTablePtr_->find(p.type());
 
     if (patchTypeCstrIter != patchConstructorTablePtr_->end())
     {
         return patchTypeCstrIter()
         (
-            pf.patch(),
-            Field<Type>::null()
+            p,
+            DimensionedField<Type, surfaceMesh>::null()
         );
     }
     else
@@ -120,14 +119,25 @@ tmp<fvsPatchField<Type> > fvsPatchField<Type>::NewCalculatedType
         (
             new calculatedFvsPatchField<Type>
             (
-                pf.patch(),
-                Field<Type>::null()
+                p,
+                DimensionedField<Type, surfaceMesh>::null()
             )
         );
     }
 }
 
 
+template<class Type>
+template<class Type2>
+tmp<fvsPatchField<Type> > fvsPatchField<Type>::NewCalculatedType
+(
+    const fvsPatchField<Type2>& pf
+)
+{
+    return NewCalculatedType(pf.patch());
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H b/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H
index 7da763f09f..88fb9a7d58 100644
--- a/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H
+++ b/src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -247,6 +247,13 @@ public:
             const dictionary&
         );
 
+        //- Return a pointer to a new calculatedFvsPatchField created on
+        //  freestore without setting patchField values
+        static tmp<fvsPatchField<Type> > NewCalculatedType
+        (
+            const fvPatch&
+        );
+
         //- Return a pointer to a new calculatedFvsPatchField created on
         //  freestore without setting patchField values
         template<class Type2>
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
index 84f04e3a94..f4f4f54f65 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -302,6 +302,43 @@ interpolate
 }
 
 
+template<class Type>
+tmp<FieldField<fvsPatchField, Type> > interpolate
+(
+    const FieldField<fvPatchField, Type>& fvpff
+)
+{
+    FieldField<fvsPatchField, Type>* fvspffPtr
+    (
+        new FieldField<fvsPatchField, Type>(fvpff.size())
+    );
+
+    forAll(*fvspffPtr, patchi)
+    {
+        fvspffPtr->set
+        (
+            patchi,
+            fvsPatchField<Type>::NewCalculatedType(fvpff[patchi].patch()).ptr()
+        );
+        (*fvspffPtr)[patchi] = fvpff[patchi];
+    }
+
+    return tmp<FieldField<fvsPatchField, Type> >(fvspffPtr);
+}
+
+
+template<class Type>
+tmp<FieldField<fvsPatchField, Type> > interpolate
+(
+    const tmp<FieldField<fvPatchField, Type> >& tfvpff
+)
+{
+    tmp<FieldField<fvPatchField, Type> > tfvspff = interpolate(tfvpff());
+    tfvpff.clear();
+    return tfvspff;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fvc
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.H b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.H
index 00c0da4f65..f6f0cf04b6 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -163,12 +163,26 @@ namespace fvc
         const tmp<GeometricField<Type, fvPatchField, volMesh> >& tvf
     );
 
-    //- Interpolate tmp field onto faces using 'interpolate(\<name\>)'
+    //- Interpolate field onto faces using 'interpolate(\<name\>)'
     template<class Type>
     static tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
     (
         const GeometricField<Type, fvPatchField, volMesh>& tvf
     );
+
+    //- Interpolate boundary field onto faces (simply a type conversion)
+    template<class Type>
+    static tmp<FieldField<fvsPatchField, Type> > interpolate
+    (
+        const FieldField<fvPatchField, Type>& fvpff
+    );
+
+    //- Interpolate boundary field onto faces (simply a type conversion)
+    template<class Type>
+    static tmp<FieldField<fvsPatchField, Type> > interpolate
+    (
+        const tmp<FieldField<fvPatchField, Type> >& tfvpff
+    );
 }
 
 
diff --git a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
index 1e0aa90e41..0f2b751bbd 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,8 +51,8 @@ void volPointInterpolation::syncUntransformedData
     const indirectPrimitivePatch& cpp = gmd.coupledPatch();
     const labelList& meshPoints = cpp.meshPoints();
 
-    const mapDistribute& slavesMap = gmd.globalPointSlavesMap();
-    const labelListList& slaves = gmd.globalPointSlaves();
+    const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
+    const labelListList& slaves = gmd.globalCoPointSlaves();
 
     List<Type> elems(slavesMap.constructSize());
     forAll(meshPoints, i)
@@ -105,8 +105,8 @@ void volPointInterpolation::pushUntransformedData
     const indirectPrimitivePatch& cpp = gmd.coupledPatch();
     const labelList& meshPoints = cpp.meshPoints();
 
-    const mapDistribute& slavesMap = gmd.globalPointSlavesMap();
-    const labelListList& slaves = gmd.globalPointSlaves();
+    const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
+    const labelListList& slaves = gmd.globalCoPointSlaves();
 
     List<Type> elems(slavesMap.constructSize());
     forAll(meshPoints, i)
@@ -155,14 +155,14 @@ void volPointInterpolation::addSeparated
             refCast<coupledPointPatchField<Type> >
                 (pf.boundaryField()[patchI]).initSwapAddSeparated
                 (
-                    Pstream::blocking,  //Pstream::nonBlocking,
+                    Pstream::nonBlocking,
                     pf.internalField()
                 );
         }
     }
 
     // Block for any outstanding requests
-    //Pstream::waitRequests();
+    Pstream::waitRequests();
 
     forAll(pf.boundaryField(), patchI)
     {
@@ -171,7 +171,7 @@ void volPointInterpolation::addSeparated
             refCast<coupledPointPatchField<Type> >
                 (pf.boundaryField()[patchI]).swapAddSeparated
                 (
-                    Pstream::blocking,  //Pstream::nonBlocking,
+                    Pstream::nonBlocking,
                     pf.internalField()
                 );
         }
@@ -306,7 +306,6 @@ void volPointInterpolation::interpolateBoundaryField
     }
 
     // Sum collocated contributions
-    //mesh().globalData().syncPointData(pfi, plusEqOp<Type>());
     syncUntransformedData(pfi, plusEqOp<Type>());
 
     // And add separated contributions
@@ -314,9 +313,7 @@ void volPointInterpolation::interpolateBoundaryField
 
     // Push master data to slaves. It is possible (not sure how often) for
     // a coupled point to have its master on a different patch so
-    // to make sure just push master data to slaves. Reuse the syncPointData
-    // structure.
-    //mesh().globalData().syncPointData(pfi, nopEqOp<Type>());
+    // to make sure just push master data to slaves.
     pushUntransformedData(pfi);
 
 
diff --git a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C
index cb47e0a101..87305f2fa8 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,10 +113,8 @@ void volPointInterpolation::calcBoundaryAddressing()
     {
         boolList oldData(isPatchPoint_);
 
-        //mesh().globalData().syncPointData(isPatchPoint_, orEqOp<bool>());
         syncUntransformedData(isPatchPoint_, orEqOp<bool>());
 
-
         forAll(isPatchPoint_, pointI)
         {
             if (isPatchPoint_[pointI] != oldData[pointI])
@@ -272,7 +270,7 @@ void volPointInterpolation::makeWeights()
     makeInternalWeights(sumWeights);
 
 
-    // Create boundary weights; add to sumWeights
+    // Create boundary weights; override sumWeights
     makeBoundaryWeights(sumWeights);
 
 
@@ -292,7 +290,6 @@ void volPointInterpolation::makeWeights()
 
 
     // Sum collocated contributions
-    //mesh().globalData().syncPointData(sumWeights, plusEqOp<scalar>());
     syncUntransformedData(sumWeights, plusEqOp<scalar>());
 
     // And add separated contributions
@@ -302,7 +299,6 @@ void volPointInterpolation::makeWeights()
     // a coupled point to have its master on a different patch so
     // to make sure just push master data to slaves. Reuse the syncPointData
     // structure.
-    //mesh().globalData().syncPointData(sumWeights, nopEqOp<scalar>());
     pushUntransformedData(sumWeights);
 
 
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index 5781ef9911..0b0e3800f7 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
@@ -713,7 +713,7 @@ Foam::label Foam::AMIInterpolation<SourcePatch, TargetPatch>::findTargetFace
 
     const pointField& srcPts = srcPatch.points();
     const face& srcFace = srcPatch[srcFaceI];
-    const point& srcPt = srcFace.centre(srcPts);
+    const point srcPt = srcFace.centre(srcPts);
     const scalar srcFaceArea = srcMagSf_[srcFaceI];
 
 //    pointIndexHit sample = treePtr_->findNearest(srcPt, sqr(0.1*bb.mag()));
@@ -782,6 +782,62 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::appendNbrFaces
 }
 
 
+template<class SourcePatch, class TargetPatch>
+bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::processSourceFace
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    const label srcFaceI,
+    const label tgtStartFaceI,
+
+    // list of tgt face neighbour faces
+    DynamicList<label>& nbrFaces,
+    // list of faces currently visited for srcFaceI to avoid multiple hits
+    DynamicList<label>& visitedFaces,
+
+    // temporary storage for addressing and weights
+    List<DynamicList<label> >& srcAddr,
+    List<DynamicList<scalar> >& srcWght,
+    List<DynamicList<label> >& tgtAddr,
+    List<DynamicList<scalar> >& tgtWght
+)
+{
+    nbrFaces.clear();
+    visitedFaces.clear();
+
+    // append initial target face and neighbours
+    nbrFaces.append(tgtStartFaceI);
+    appendNbrFaces(tgtStartFaceI, tgtPatch, visitedFaces, nbrFaces);
+
+    bool faceProcessed = false;
+
+    do
+    {
+        // process new target face
+        label tgtFaceI = nbrFaces.remove();
+        visitedFaces.append(tgtFaceI);
+        scalar area = interArea(srcFaceI, tgtFaceI, srcPatch, tgtPatch);
+
+        // store when intersection area > 0
+        if (area > 0)
+        {
+            srcAddr[srcFaceI].append(tgtFaceI);
+            srcWght[srcFaceI].append(area);
+
+            tgtAddr[tgtFaceI].append(srcFaceI);
+            tgtWght[tgtFaceI].append(area);
+
+            appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
+
+            faceProcessed = true;
+        }
+
+    } while (nbrFaces.size() > 0);
+
+    return faceProcessed;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
 (
@@ -811,7 +867,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
                 label faceT = visitedFaces[j];
                 scalar area = interArea(faceS, faceT, srcPatch0, tgtPatch0);
 
-                if (area > 0)
+                // Check that faces have enough overlap for robust walking
+                if (area/srcMagSf_[srcFaceI] > faceAreaIntersect::tolerance())
                 {
                     // TODO - throwing area away - re-use in next iteration?
 
@@ -864,7 +921,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
                 << "target face" << endl;
         }
 
-//        foundNextSeed = false;
+        foundNextSeed = false;
         for (label faceI = startSeedI; faceI < mapFlag.size(); faceI++)
         {
             if (mapFlag[faceI])
@@ -887,7 +944,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
 
         FatalErrorIn
         (
-            "void Foam::cyclicAMIPolyPatch::setNextFaces"
+            "void Foam::AMIInterpolation<SourcePatch, TargetPatch>::"
+            "setNextFaces"
             "("
                 "label&, "
                 "label&, "
@@ -946,6 +1004,25 @@ Foam::scalar Foam::AMIInterpolation<SourcePatch, TargetPatch>::interArea
     {
         area = inter.calc(src, tgt, n, triMode_);
     }
+    else
+    {
+        WarningIn
+        (
+            "void Foam::AMIInterpolation<SourcePatch, TargetPatch>::"
+            "interArea"
+            "("
+                "const label, "
+                "const label, "
+                "const SourcePatch&, "
+                "const TargetPatch&"
+            ") const"
+        )   << "Invalid normal for source face " << srcFaceI
+            << " points " << UIndirectList<point>(srcPoints, src)
+            << " target face " << tgtFaceI
+            << " points " << UIndirectList<point>(tgtPoints, tgt)
+            << endl;
+    }
+
 
     if ((debug > 1) && (area > 0))
     {
@@ -956,6 +1033,105 @@ Foam::scalar Foam::AMIInterpolation<SourcePatch, TargetPatch>::interArea
 }
 
 
+template<class SourcePatch, class TargetPatch>
+void Foam::AMIInterpolation<SourcePatch, TargetPatch>::
+restartUncoveredSourceFace
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    List<DynamicList<label> >& srcAddr,
+    List<DynamicList<scalar> >& srcWght,
+    List<DynamicList<label> >& tgtAddr,
+    List<DynamicList<scalar> >& tgtWght
+)
+{
+    // Collect all src faces with a low weight
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelHashSet lowWeightFaces(100);
+    forAll(srcWght, srcFaceI)
+    {
+        scalar s = sum(srcWght[srcFaceI]);
+        scalar t = s/srcMagSf_[srcFaceI];
+
+        if (t < 0.5)
+        {
+            lowWeightFaces.insert(srcFaceI);
+        }
+    }
+
+    Info<< "AMIInterpolation : restarting search on "
+        << returnReduce(lowWeightFaces.size(), sumOp<label>())
+        << " faces since sum of weights < 0.5" << endl;
+
+    if (lowWeightFaces.size() > 0)
+    {
+        // Erase all the lowWeight source faces from the target
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        DynamicList<label> okSrcFaces(10);
+        DynamicList<scalar> okSrcWeights(10);
+        forAll(tgtAddr, tgtFaceI)
+        {
+            okSrcFaces.clear();
+            okSrcWeights.clear();
+            DynamicList<label>& srcFaces = tgtAddr[tgtFaceI];
+            DynamicList<scalar>& srcWeights = tgtWght[tgtFaceI];
+            forAll(srcFaces, i)
+            {
+                if (!lowWeightFaces.found(srcFaces[i]))
+                {
+                    okSrcFaces.append(srcFaces[i]);
+                    okSrcWeights.append(srcWeights[i]);
+                }
+            }
+            if (okSrcFaces.size() < srcFaces.size())
+            {
+                srcFaces.transfer(okSrcFaces);
+                srcWeights.transfer(okSrcWeights);
+            }
+        }
+
+
+
+        // Restart search from best hit
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        // list of tgt face neighbour faces
+        DynamicList<label> nbrFaces(10);
+
+        // list of faces currently visited for srcFaceI to avoid multiple hits
+        DynamicList<label> visitedFaces(10);
+
+        forAllConstIter(labelHashSet, lowWeightFaces, iter)
+        {
+            label srcFaceI = iter.key();
+            label tgtFaceI = findTargetFace(srcFaceI, srcPatch);
+            if (tgtFaceI != -1)
+            {
+                //bool faceProcessed =
+                processSourceFace
+                (
+                    srcPatch,
+                    tgtPatch,
+                    srcFaceI,
+                    tgtFaceI,
+
+                    nbrFaces,
+                    visitedFaces,
+
+                    srcAddr,
+                    srcWght,
+                    tgtAddr,
+                    tgtWght
+                );
+                // ? Check faceProcessed to see if restarting has worked.
+            }
+        }
+    }
+}
+
+
 template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
 (
@@ -1073,37 +1249,22 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
     DynamicList<label> nonOverlapFaces;
     do
     {
-        nbrFaces.clear();
-        visitedFaces.clear();
+        // Do advancing front starting from srcFaceI,tgtFaceI
+        bool faceProcessed = processSourceFace
+        (
+            srcPatch,
+            tgtPatch,
+            srcFaceI,
+            tgtFaceI,
 
-        // append initial target face and neighbours
-        nbrFaces.append(tgtFaceI);
-        appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
+            nbrFaces,
+            visitedFaces,
 
-        bool faceProcessed = false;
-
-        do
-        {
-            // process new target face
-            tgtFaceI = nbrFaces.remove();
-            visitedFaces.append(tgtFaceI);
-            scalar area = interArea(srcFaceI, tgtFaceI, srcPatch, tgtPatch);
-
-            // store when intersection area > 0
-            if (area > 0)
-            {
-                srcAddr[srcFaceI].append(tgtFaceI);
-                srcWght[srcFaceI].append(area);
-
-                tgtAddr[tgtFaceI].append(srcFaceI);
-                tgtWght[tgtFaceI].append(area);
-
-                appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
-
-                faceProcessed = true;
-            }
-
-        } while (nbrFaces.size() > 0);
+            srcAddr,
+            srcWght,
+            tgtAddr,
+            tgtWght
+        );
 
         mapFlag[srcFaceI] = false;
 
@@ -1140,6 +1301,22 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
         srcNonOverlap_.transfer(nonOverlapFaces);
     }
 
+
+    // Check for badly covered faces
+    if (debug)
+    {
+        restartUncoveredSourceFace
+        (
+            srcPatch,
+            tgtPatch,
+            srcAddr,
+            srcWght,
+            tgtAddr,
+            tgtWght
+        );
+    }
+
+
     // transfer data to persistent storage
     forAll(srcAddr, i)
     {
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
index b393030823..5b91afc38a 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
@@ -246,6 +246,31 @@ class AMIInterpolation
                 DynamicList<label>& faceIDs
             ) const;
 
+            bool processSourceFace
+            (
+                const SourcePatch& srcPatch,
+                const TargetPatch& tgtPatch,
+                const label srcFaceI,
+                const label tgtStartFaceI,
+
+                DynamicList<label>& nbrFaces,
+                DynamicList<label>& visitedFaces,
+                List<DynamicList<label> >& srcAddr,
+                List<DynamicList<scalar> >& srcWght,
+                List<DynamicList<label> >& tgtAddr,
+                List<DynamicList<scalar> >& tgtWght
+            );
+
+            void restartUncoveredSourceFace
+            (
+                const SourcePatch& srcPatch,
+                const TargetPatch& tgtPatch,
+                List<DynamicList<label> >& srcAddr,
+                List<DynamicList<scalar> >& srcWght,
+                List<DynamicList<label> >& tgtAddr,
+                List<DynamicList<scalar> >& tgtWght
+            );
+
             //- Set the source and target seed faces
             void setNextFaces
             (
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
index db8a0ba294..4135cfe041 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,10 @@ License
 
 #include "faceAreaIntersect.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::scalar Foam::faceAreaIntersect::tol = 1e-6;
+
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 void Foam::faceAreaIntersect::triSliceWithPlane
@@ -36,8 +40,6 @@ void Foam::faceAreaIntersect::triSliceWithPlane
     const scalar len
 )
 {
-    const scalar matchTol = 1e-6;
-
     // distance to cutting plane
     FixedList<scalar, 3> d;
 
@@ -51,7 +53,7 @@ void Foam::faceAreaIntersect::triSliceWithPlane
     {
         d[i] = ((tri[i] - p.refPoint()) & p.normal());
 
-        if (mag(d[i]) < matchTol*len)
+        if (mag(d[i]) < tol*len)
         {
             nCoPlanar++;
             copI = i;
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
index 2a27e30960..bfbe7c08a2 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,6 +78,11 @@ private:
         const bool reverseB_;
 
 
+    // Static data members
+
+        static scalar tol;
+
+
     // Private Member Functions
 
         //- Get triPoints from face
@@ -152,6 +157,9 @@ public:
 
     // Public Member Functions
 
+        //- Fraction of local length scale to use as intersection tolerance
+        inline static scalar& tolerance();
+
         //- Return area of intersection of faceA with faceB
         scalar calc
         (
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
index 0b8265aad6..87fab781e0 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,6 +23,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
 inline void Foam::faceAreaIntersect::setTriPoints
 (
     const point& a,
@@ -106,4 +108,12 @@ inline Foam::scalar Foam::faceAreaIntersect::triArea(const triPoints& t) const
 }
 
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar& Foam::faceAreaIntersect::tolerance()
+{
+    return tol;
+}
+
+
 // ************************************************************************* //
diff --git a/src/postProcessing/foamCalcFunctions/field/magSqr/magSqr.C b/src/postProcessing/foamCalcFunctions/field/magSqr/magSqr.C
index 084f87de4b..a608ade6bf 100644
--- a/src/postProcessing/foamCalcFunctions/field/magSqr/magSqr.C
+++ b/src/postProcessing/foamCalcFunctions/field/magSqr/magSqr.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ Foam::calcTypes::magSqr::~magSqr()
 
 void Foam::calcTypes::magSqr::init()
 {
-    Foam::argList::validArgs.append("magSqr");
+    argList::validArgs.append("magSqr");
     argList::validArgs.append("fieldName");
 }
 
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
index b3aaac0dbf..8dbf0ec5a8 100644
--- a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
@@ -27,8 +27,8 @@ License
 #include "volFields.H"
 #include "regionSplit.H"
 #include "fvcVolumeIntegrate.H"
-#include "Histogram.H"
 #include "mathematicalConstants.H"
+#include "stringListOps.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -56,6 +56,281 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::regionSizeDistribution::writeGraph
+(
+    const coordSet& coords,
+    const word& valueName,
+    const scalarField& values
+) const
+{
+    const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+    const wordList valNames(1, valueName);
+
+    fileName outputPath;
+    if (Pstream::parRun())
+    {
+        outputPath = mesh.time().path()/".."/name_;
+    }
+    else
+    {
+        outputPath = mesh.time().path()/name_;
+    }
+
+    if (mesh.name() != fvMesh::defaultRegion)
+    {
+        outputPath = outputPath/mesh.name();
+    }
+
+    mkDir(outputPath/mesh.time().timeName());
+    OFstream str
+    (
+        outputPath
+      / mesh.time().timeName()
+      / formatterPtr_().getFileName(coords, valNames)
+    );
+    Info<< "Writing distribution of " << valueName << " to " << str.name()
+        << endl;
+
+    List<const scalarField*> valPtrs(1);
+    valPtrs[0] = &values;
+    formatterPtr_().write(coords, valNames, valPtrs, str);
+}
+
+
+void Foam::regionSizeDistribution::writeAlphaFields
+(
+    const regionSplit& regions,
+    const Map<label>& patchRegions,
+    const Map<scalar>& regionVolume,
+    const volScalarField& alpha
+) const
+{
+    const scalar maxDropletVol = 1.0/6.0*pow(maxDiam_, 3);
+
+    // Split alpha field
+    // ~~~~~~~~~~~~~~~~~
+    // Split into
+    //  - liquidCore            : region connected to inlet patches
+    //  - per region a volume   : for all other regions
+    //  - backgroundAlpha       : remaining alpha
+
+
+    // Construct field
+    volScalarField liquidCore
+    (
+        IOobject
+        (
+            alphaName_ + "_liquidCore",
+            obr_.time().timeName(),
+            obr_,
+            IOobject::NO_READ
+        ),
+        alpha,
+        fvPatchField<scalar>::calculatedType()
+    );
+
+    volScalarField backgroundAlpha
+    (
+        IOobject
+        (
+            alphaName_ + "_background",
+            obr_.time().timeName(),
+            obr_,
+            IOobject::NO_READ
+        ),
+        alpha,
+        fvPatchField<scalar>::calculatedType()
+    );
+
+
+    // Knock out any cell not in patchRegions
+    forAll(liquidCore, cellI)
+    {
+        label regionI = regions[cellI];
+        if (patchRegions.found(regionI))
+        {
+            backgroundAlpha[cellI] = 0;
+        }
+        else
+        {
+            liquidCore[cellI] = 0;
+
+            scalar regionVol = regionVolume[regionI];
+            if (regionVol < maxDropletVol)
+            {
+                backgroundAlpha[cellI] = 0;
+            }
+        }
+    }
+    liquidCore.correctBoundaryConditions();
+    backgroundAlpha.correctBoundaryConditions();
+
+    Info<< "Volume of liquid-core = "
+        << fvc::domainIntegrate(liquidCore).value()
+        << endl;
+    Info<< "Volume of background  = "
+        << fvc::domainIntegrate(backgroundAlpha).value()
+        << endl;
+
+    Info<< "Writing liquid-core field to " << liquidCore.name() << endl;
+    liquidCore.write();
+    Info<< "Writing background field to " << backgroundAlpha.name() << endl;
+    backgroundAlpha.write();
+}
+
+
+Foam::Map<Foam::label> Foam::regionSizeDistribution::findPatchRegions
+(
+    const polyMesh& mesh,
+    const regionSplit& regions
+) const
+{
+    // Mark all regions starting at patches
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Count number of patch faces (just for initial sizing)
+    const labelHashSet patchIDs(mesh.boundaryMesh().patchSet(patchNames_));
+
+    label nPatchFaces = 0;
+    forAllConstIter(labelHashSet, patchIDs, iter)
+    {
+        nPatchFaces += mesh.boundaryMesh()[iter.key()].size();
+    }
+
+
+    Map<label> patchRegions(nPatchFaces);
+    forAllConstIter(labelHashSet, patchIDs, iter)
+    {
+        const polyPatch& pp = mesh.boundaryMesh()[iter.key()];
+
+        // Collect all regions on the patch
+        const labelList& faceCells = pp.faceCells();
+
+        forAll(faceCells, i)
+        {
+            patchRegions.insert
+            (
+                regions[faceCells[i]],
+                Pstream::myProcNo()     // dummy value
+            );
+        }
+    }
+
+
+    // Make sure all the processors have the same set of regions
+    Pstream::mapCombineGather(patchRegions, minEqOp<label>());
+    Pstream::mapCombineScatter(patchRegions);
+
+    return patchRegions;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::regionSizeDistribution::divide
+(
+    const scalarField& num,
+    const scalarField& denom
+)
+{
+    tmp<scalarField> tresult(new scalarField(num.size()));
+    scalarField& result = tresult();
+
+    forAll(denom, i)
+    {
+        if (denom[i] != 0)
+        {
+            result[i] = num[i]/denom[i];
+        }
+        else
+        {
+            result[i] = 0.0;
+        }
+    }
+    return tresult;
+}
+
+
+void Foam::regionSizeDistribution::writeGraphs
+(
+    const word& fieldName,              // name of field
+    const labelList& indices,           // index of bin for each region
+    const scalarField& sortedField,     // per region field data
+    const scalarField& binCount,        // per bin number of regions
+    const coordSet& coords              // graph data for bins
+) const
+{
+    if (Pstream::master())
+    {
+        // Calculate per-bin average
+        scalarField binSum(nBins_, 0.0);
+        forAll(sortedField, i)
+        {
+            binSum[indices[i]] += sortedField[i];
+        }
+
+        scalarField binAvg(divide(binSum, binCount));
+
+        // Per bin deviation
+        scalarField binSqrSum(nBins_, 0.0);
+        forAll(sortedField, i)
+        {
+            binSqrSum[indices[i]] += Foam::sqr(sortedField[i]);
+        }
+        scalarField binDev
+        (
+            sqrt(divide(binSqrSum, binCount) - Foam::sqr(binAvg))
+        );
+
+        // Write average
+        writeGraph(coords, fieldName + "_sum", binSum);
+        // Write average
+        writeGraph(coords, fieldName + "_avg", binAvg);
+        // Write deviation
+        writeGraph(coords, fieldName + "_dev", binDev);
+    }
+}
+
+
+void Foam::regionSizeDistribution::writeGraphs
+(
+    const word& fieldName,              // name of field
+    const scalarField& cellField,       // per cell field data
+    const regionSplit& regions,         // per cell the region(=droplet)
+    const labelList& sortedRegions,     // valid regions in sorted order
+    const scalarField& sortedNormalisation,
+
+    const labelList& indices,           // per region index of bin
+    const scalarField& binCount,        // per bin number of regions
+    const coordSet& coords              // graph data for bins
+) const
+{
+    // Sum on a per-region basis. Parallel reduced.
+    Map<scalar> regionField(regionSum(regions, cellField));
+
+    // Extract in region order
+    scalarField sortedField
+    (
+        sortedNormalisation
+      * extractData
+        (
+            sortedRegions,
+            regionField
+        )
+    );
+
+    writeGraphs
+    (
+        fieldName,      // name of field
+        indices,        // index of bin for each region
+        sortedField,    // per region field data
+        binCount,       // per bin number of regions
+        coords          // graph data for bins
+    );
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::regionSizeDistribution::regionSizeDistribution
@@ -103,8 +378,11 @@ void Foam::regionSizeDistribution::read(const dictionary& dict)
         dict.lookup("field") >> alphaName_;
         dict.lookup("patches") >> patchNames_;
         dict.lookup("threshold") >> threshold_;
-        dict.lookup("volFraction") >> volFraction_;
+        dict.lookup("maxDiameter") >> maxDiam_;
+        minDiam_ = 0.0;
+        dict.readIfPresent("minDiameter", minDiam_);
         dict.lookup("nBins") >> nBins_;
+        dict.lookup("fields") >> fields_;
 
         word format(dict.lookup("setFormat"));
         formatterPtr_ = writer<scalar>::New(format);
@@ -163,14 +441,17 @@ void Foam::regionSizeDistribution::write()
            : obr_.lookupObject<volScalarField>(alphaName_)
         );
 
-        Info<< "Volume of alpha = "
+        Info<< "Volume of alpha          = "
             << fvc::domainIntegrate(alpha).value()
             << endl;
 
         const scalar meshVol = gSum(mesh.V());
-        Info<< "Mesh volume = " << meshVol << endl;
-        Info<< "Background region volume limit = " << volFraction_*meshVol
-            << endl;
+        const scalar maxDropletVol = 1.0/6.0*pow(maxDiam_, 3);
+        const scalar delta = (maxDiam_-minDiam_)/nBins_;
+
+        Info<< "Mesh volume              = " << meshVol << endl;
+        Info<< "Maximum droplet diameter = " << maxDiam_ << endl;
+        Info<< "Maximum droplet volume   = " << maxDropletVol << endl;
 
 
         // Determine blocked faces
@@ -260,138 +541,115 @@ void Foam::regionSizeDistribution::write()
         }
 
 
+        // Determine regions connected to supplied patches
+        Map<label> patchRegions(findPatchRegions(mesh, regions));
+
+
+
         // Sum all regions
-        Map<Pair<scalar> > regionVolume(regions.nRegions()/Pstream::nProcs());
-        forAll(alpha, cellI)
-        {
-            scalar cellVol = mesh.V()[cellI];
-            scalar alphaVol = alpha[cellI]*cellVol;
-
-            label regionI = regions[cellI];
-
-            Map<Pair<scalar> >::iterator fnd = regionVolume.find(regionI);
-            if (fnd == regionVolume.end())
-            {
-                regionVolume.insert
-                (
-                    regionI,
-                    Pair<scalar>(cellVol, alphaVol)
-                );
-            }
-            else
-            {
-                fnd().first() += cellVol;
-                fnd().second() += alphaVol;
-            }
-        }
-        Pstream::mapCombineGather(regionVolume, ListPlusEqOp<scalar, 2>());
-        Pstream::mapCombineScatter(regionVolume);
-
+        const scalarField alphaVol = alpha.internalField()*mesh.V();
+        Map<scalar> allRegionVolume(regionSum(regions, mesh.V()));
+        Map<scalar> allRegionAlphaVolume(regionSum(regions, alphaVol));
+        Map<label> allRegionNumCells
+        (
+            regionSum
+            (
+                regions,
+                labelField(mesh.nCells(), 1.0)
+            )
+        );
 
         if (debug)
         {
             Info<< token::TAB << "Region"
                 << token::TAB << "Volume(mesh)"
                 << token::TAB << "Volume(" << alpha.name() << "):"
+                << token::TAB << "nCells"
                 << endl;
             scalar meshSumVol = 0.0;
             scalar alphaSumVol = 0.0;
+            label nCells = 0;
 
-            forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
+            Map<scalar>::const_iterator vIter = allRegionVolume.begin();
+            Map<scalar>::const_iterator aIter = allRegionAlphaVolume.begin();
+            Map<label>::const_iterator numIter = allRegionNumCells.begin();
+            for
+            (
+                ;
+                vIter != allRegionVolume.end()
+             && aIter != allRegionAlphaVolume.end();
+                ++vIter, ++aIter, ++numIter
+            )
             {
-                Info<< token::TAB << iter.key()
-                    << token::TAB << iter().first()
-                    << token::TAB << iter().second() << endl;
+                Info<< token::TAB << vIter.key()
+                    << token::TAB << vIter()
+                    << token::TAB << aIter()
+                    << token::TAB << numIter()
+                    << endl;
 
-                meshSumVol += iter().first();
-                alphaSumVol += iter().second();
+                meshSumVol += vIter();
+                alphaSumVol += aIter();
+                nCells += numIter();
             }
             Info<< token::TAB << "Total:"
                 << token::TAB << meshSumVol
-                << token::TAB << alphaSumVol << endl;
+                << token::TAB << alphaSumVol
+                << token::TAB << nCells
+                << endl;
             Info<< endl;
         }
 
 
 
-        // Mark all regions starting at patches
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        // Count number of patch faces (just for initial sizing)
-        label nPatchFaces = 0;
-        forAll(patchNames_, i)
         {
-            const word& pName = patchNames_[i];
-            label patchI = mesh.boundaryMesh().findPatchID(pName);
-            if (patchI == -1)
+            Info<< "Patch connected regions (liquid core):" << endl;
+            Info<< token::TAB << "Region"
+                << token::TAB << "Volume(mesh)"
+                << token::TAB << "Volume(" << alpha.name() << "):"
+                << endl;
+            forAllConstIter(Map<label>, patchRegions, iter)
             {
-                WarningIn("regionSizeDistribution::write()")
-                    << "Cannot find patch " << pName << ". Valid patches are "
-                    << mesh.boundaryMesh().names()
-                    << endl;
-            }
-            else
-            {
-                nPatchFaces += mesh.boundaryMesh()[patchI].size();
+                label regionI = iter.key();
+                Info<< token::TAB << iter.key()
+                    << token::TAB << allRegionVolume[regionI]
+                    << token::TAB << allRegionAlphaVolume[regionI] << endl;
+
             }
+            Info<< endl;
         }
 
-        Map<label> keepRegions(nPatchFaces);
-        forAll(patchNames_, i)
         {
-            const word& pName = patchNames_[i];
+            Info<< "Background regions:" << endl;
+            Info<< token::TAB << "Region"
+                << token::TAB << "Volume(mesh)"
+                << token::TAB << "Volume(" << alpha.name() << "):"
+                << endl;
+            Map<scalar>::const_iterator vIter = allRegionVolume.begin();
+            Map<scalar>::const_iterator aIter = allRegionAlphaVolume.begin();
 
-            label patchI = mesh.boundaryMesh().findPatchID(pName);
-            if (patchI != -1)
-            {
-                const polyPatch& pp = mesh.boundaryMesh()[patchI];
-
-                // Collect all regions on the patch
-                const labelList& faceCells = pp.faceCells();
-
-                forAll(faceCells, i)
-                {
-                    keepRegions.insert
-                    (
-                        regions[faceCells[i]],
-                        Pstream::myProcNo()
-                    );
-                }
-            }
-        }
-
-
-        // Make sure all the processors have the same set of regions
-        Pstream::mapCombineGather(keepRegions, minEqOp<label>());
-        Pstream::mapCombineScatter(keepRegions);
-
-        Info<< "Patch connected regions (liquid core):" << endl;
-        forAllConstIter(Map<label>, keepRegions, iter)
-        {
-            label regionI = iter.key();
-            Pair<scalar>& vols = regionVolume[regionI];
-            Info<< token::TAB << iter.key()
-                << token::TAB << vols.first()
-                << token::TAB << vols.second() << endl;
-
-        }
-        Info<< endl;
-
-        Info<< "Background regions:" << endl;
-        forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
-        {
-            if
+            for
             (
-               !keepRegions.found(iter.key())
-             && iter().first() >= volFraction_*meshVol
+                ;
+                vIter != allRegionVolume.end()
+             && aIter != allRegionAlphaVolume.end();
+                ++vIter, ++aIter
             )
             {
-                Info<< token::TAB << iter.key()
-                    << token::TAB << iter().first()
-                    << token::TAB << iter().second() << endl;
+                if
+                (
+                   !patchRegions.found(vIter.key())
+                 && vIter() >= maxDropletVol
+                )
+                {
+                    Info<< token::TAB << vIter.key()
+                        << token::TAB << vIter()
+                        << token::TAB << aIter() << endl;
+                }
             }
+            Info<< endl;
         }
-        Info<< endl;
+
 
 
         // Split alpha field
@@ -400,185 +658,197 @@ void Foam::regionSizeDistribution::write()
         //  - liquidCore            : region connected to inlet patches
         //  - per region a volume   : for all other regions
         //  - backgroundAlpha       : remaining alpha
+        writeAlphaFields(regions, patchRegions, allRegionVolume, alpha);
 
 
-        // Construct field
-        volScalarField liquidCore
-        (
-            IOobject
-            (
-                alphaName_ + "_liquidCore",
-                obr_.time().timeName(),
-                obr_,
-                IOobject::NO_READ
-            ),
-            alpha,
-            fvPatchField<scalar>::calculatedType()
-        );
-
-        volScalarField backgroundAlpha
-        (
-            IOobject
-            (
-                alphaName_ + "_background",
-                obr_.time().timeName(),
-                obr_,
-                IOobject::NO_READ
-            ),
-            alpha,
-            fvPatchField<scalar>::calculatedType()
-        );
-
-
-        // Knock out any cell not in keepRegions
-        forAll(liquidCore, cellI)
+        // Extract droplet-only allRegionVolume, i.e. delete liquid core
+        // (patchRegions) and background regions from maps.
+        // Note that we have to use mesh volume (allRegionVolume) and not
+        // allRegionAlphaVolume since background might not have alpha in it.
+        forAllIter(Map<scalar>, allRegionVolume, vIter)
         {
-            label regionI = regions[cellI];
-            if (keepRegions.found(regionI))
+            label regionI = vIter.key();
+            if
+            (
+                patchRegions.found(regionI)
+             || vIter() >= maxDropletVol
+            )
             {
-                backgroundAlpha[cellI] = 0;
-            }
-            else
-            {
-                liquidCore[cellI] = 0;
-
-                scalar regionVol = regionVolume[regionI].first();
-                if (regionVol < volFraction_*meshVol)
-                {
-                    backgroundAlpha[cellI] = 0;
-                }
+                allRegionVolume.erase(vIter);
+                allRegionAlphaVolume.erase(regionI);
+                allRegionNumCells.erase(regionI);
             }
         }
-        liquidCore.correctBoundaryConditions();
-        backgroundAlpha.correctBoundaryConditions();
 
-        Info<< "Volume of liquid-core = "
-            << fvc::domainIntegrate(liquidCore).value()
-            << endl;
-
-        Info<< "Writing liquid-core field to " << liquidCore.name() << endl;
-        liquidCore.write();
-
-        Info<< "Volume of background = "
-            << fvc::domainIntegrate(backgroundAlpha).value()
-            << endl;
-
-        Info<< "Writing background field to " << backgroundAlpha.name() << endl;
-        backgroundAlpha.write();
-
-
-
-        // Collect histogram
-        if (Pstream::master())
+        if (allRegionVolume.size())
         {
-            DynamicList<scalar> diameters(regionVolume.size());
-            forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
+            // Construct mids of bins for plotting
+            pointField xBin(nBins_);
+
+            scalar x = 0.5*delta;
+            forAll(xBin, i)
             {
-                if (!keepRegions.found(iter.key()))
-                {
-                    if (iter().first() < volFraction_*meshVol)
-                    {
-                        scalar v = iter().second();
-                      //scalar diam = Foam::cbrt(v*6/mathematicalConstant::pi);
-                        scalar diam =
-                            Foam::cbrt(v*6/constant::mathematical::pi);
-                        diameters.append(diam);
-                    }
-                }
+                xBin[i] = point(x, 0, 0);
+                x += delta;
             }
 
-            if (diameters.size())
+            const coordSet coords("diameter", "x", xBin, mag(xBin));
+
+
+            // Get in region order the alpha*volume and diameter
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            const labelList sortedRegions = allRegionAlphaVolume.sortedToc();
+
+            scalarField sortedVols
+            (
+                extractData
+                (
+                    sortedRegions,
+                    allRegionAlphaVolume
+                )
+            );
+
+            // Calculate the diameters
+            scalarField sortedDiameters(sortedVols.size());
+            forAll(sortedDiameters, i)
             {
-                scalar maxDiam = max(diameters);
-                scalar minDiam = 0.0;
-
-                Info<< "Maximum diameter:" << maxDiam << endl;
-
-                Histogram<List<scalar> > bins
+                sortedDiameters[i] = Foam::cbrt
                 (
-                    minDiam,
-                    maxDiam,
-                    nBins_,
-                    diameters
+                    sortedVols[i]
+                   *6/constant::mathematical::pi
                 );
+            }
 
-                /* 1.7.x
-                scalarField xBin(nBins_);
+            // Determine the bin index for all the diameters
+            labelList indices(sortedDiameters.size());
+            forAll(sortedDiameters, i)
+            {
+                indices[i] = (sortedDiameters[i]-minDiam_)/delta;
+            }
 
-                scalar dx = (maxDiam-minDiam)/nBins_;
-                scalar x = 0.5*dx;
-                forAll(bins.counts(), i)
+            // Calculate the counts per diameter bin
+            scalarField binCount(nBins_, 0.0);
+            forAll(sortedDiameters, i)
+            {
+                binCount[indices[i]] += 1.0;
+            }
+
+            // Write counts
+            if (Pstream::master())
+            {
+                writeGraph(coords, "count", binCount);
+            }
+
+            // Write to screen
+            {
+                Info<< "Bins:" << endl;
+                Info<< token::TAB << "Bin"
+                    << token::TAB << "Min diameter"
+                    << token::TAB << "Count:"
+                    << endl;
+
+                scalar diam = 0.0;
+                forAll(binCount, binI)
                 {
-                    xBin[i] = x;
-                    x += dx;
+                    Info<< token::TAB << binI
+                        << token::TAB << diam
+                        << token::TAB << binCount[binI] << endl;
+                    diam += delta;
                 }
+                Info<< endl;
+            }
 
-                scalarField normalisedCount(bins.counts().size());
-                forAll(bins.counts(), i)
+
+            // Write average and deviation of droplet volume.
+            writeGraphs
+            (
+                "volume",           // name of field
+                indices,            // per region the bin index
+                sortedVols,         // per region field data
+                binCount,           // per bin number of regions
+                coords              // graph data for bins
+            );
+
+            // Collect some more field
+            {
+                wordList scalarNames(obr_.names(volScalarField::typeName));
+                labelList selected = findStrings(fields_, scalarNames);
+
+                forAll(selected, i)
                 {
-                    normalisedCount[i] = 1.0*bins.counts()[i];
+                    const word& fldName = scalarNames[selected[i]];
+                    Info<< "Scalar field " << fldName << endl;
+
+                    const scalarField& fld = obr_.lookupObject
+                    <
+                        volScalarField
+                    >(fldName).internalField();
+
+                    writeGraphs
+                    (
+                        fldName,            // name of field
+                        alphaVol*fld,       // per cell field data
+
+                        regions,            // per cell the region(=droplet)
+                        sortedRegions,      // valid regions in sorted order
+                        1.0/sortedVols,     // per region normalisation
+
+                        indices,            // index of bin for each region
+                        binCount,           // per bin number of regions
+                        coords              // graph data for bins
+                    );
                 }
+            }
+            {
+                wordList vectorNames(obr_.names(volVectorField::typeName));
+                labelList selected = findStrings(fields_, vectorNames);
 
-                const coordSet coords
-                (
-                    "diameter",
-                    "x",
-                    xBin
-                );
-                */
-
-                pointField xBin(nBins_);
-                scalar dx = (maxDiam - minDiam)/nBins_;
-                scalar x = 0.5*dx;
-                forAll(bins.counts(), i)
+                forAll(selected, i)
                 {
-                    xBin[i] = point(x, 0, 0);
-                    x += dx;
+                    const word& fldName = vectorNames[selected[i]];
+                    Info<< "Vector field " << fldName << endl;
+
+                    const vectorField& fld = obr_.lookupObject
+                    <
+                        volVectorField
+                    >(fldName).internalField();
+
+
+                    // Components
+
+                    for (direction cmp = 0; cmp < vector::nComponents; cmp++)
+                    {
+                        writeGraphs
+                        (
+                            fldName + vector::componentNames[cmp],
+                            alphaVol*fld.component(cmp),// per cell field data
+
+                            regions,        // per cell the region(=droplet)
+                            sortedRegions,  // valid regions in sorted order
+                            1.0/sortedVols, // per region normalisation
+
+                            indices,        // index of bin for each region
+                            binCount,       // per bin number of regions
+                            coords          // graph data for bins
+                        );
+                    }
+
+                    // Magnitude
+                    writeGraphs
+                    (
+                        fldName + "mag",    // name of field
+                        alphaVol*mag(fld),  // per cell field data
+
+                        regions,            // per cell the region(=droplet)
+                        sortedRegions,      // valid regions in sorted order
+                        1.0/sortedVols,     // per region normalisation
+
+                        indices,            // index of bin for each region
+                        binCount,           // per bin number of regions
+                        coords              // graph data for bins
+                    );
                 }
-
-                scalarField normalisedCount(bins.counts().size());
-                forAll(bins.counts(), i)
-                {
-                    normalisedCount[i] = 1.0*bins.counts()[i];
-                }
-
-                const coordSet coords
-                (
-                    "diameter",
-                    "x",
-                    xBin,
-                    mag(xBin)
-                );
-                const wordList valNames(1, "count");
-
-
-                fileName outputPath;
-                if (Pstream::parRun())
-                {
-                    outputPath = mesh.time().path()/".."/name_;
-                }
-                else
-                {
-                    outputPath = mesh.time().path()/name_;
-                }
-
-                if (mesh.name() != fvMesh::defaultRegion)
-                {
-                    outputPath = outputPath/mesh.name();
-                }
-
-                mkDir(outputPath/mesh.time().timeName());
-                OFstream str
-                (
-                    outputPath
-                  / mesh.time().timeName()
-                  / formatterPtr_().getFileName(coords, valNames)
-                );
-                Info<< "Writing distribution to " << str.name() << endl;
-
-                List<const scalarField*> valPtrs(1);
-                valPtrs[0] = &normalisedCount;
-                formatterPtr_().write(coords, valNames, valPtrs, str);
             }
         }
     }
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
index 033b5f9aec..7802733249 100644
--- a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
@@ -25,11 +25,67 @@ Class
     Foam::regionSizeDistribution
 
 Description
-    Looks up a field, interpolates it to the faces and determines a connected
-    region from a patch where the field is above a certain value.
-    - Writes a field containing all regions starting at given patch
-      ('liquid core')
-    - All other regions are summed for volume and a histogram is calculated.
+    Droplet size distribution calculation.
+
+    Looks up a void-fraction (alpha) field and splits the mesh into regions
+    based on where the field is below the threshold value. These
+    regions ("droplets") can now be analysed.
+
+    Regions:
+    - (debug) write regions as a volScalarField
+    - (debug) print for all regions the sum of volume and alpha*volume
+    - print the regions connected to a user-defined set of patches.
+      (in spray calculation these form the liquid core)
+    - print the regions with too large volume. These are the 'background'
+      regions.
+
+    Fields:
+    - write volScalarField alpha_liquidCore : alpha with outside liquid core
+                                              set to 0.
+                           alpha_background : alpha with outside background
+                                              set to 0.
+
+    Histogram:
+    - determine histogram of diameter (given minDiameter, maxDiameter, nBins)
+    - write graph of number of droplets per bin
+    - write graph of sum, average and deviation of droplet volume per bin
+    - write graph of sum, average and deviation of user-defined fields. For
+      volVectorFields these are those of the 3 components and the magnitude.
+
+    Sample input:
+
+    functions
+    {
+        regionSizeDistribution
+        {
+            type            regionSizeDistribution;
+
+            outputControl   timeStep;
+            outputInterval  1;
+
+            // Field to determine regions from
+            field           alpha;
+            // Patches that provide the liquid core
+            patches         (inlet);
+            // Delimit alpha regions
+            threshold       0.4;
+
+            // Fields to sample (no need to include alpha)
+            fields          (p U);
+
+            // Number of bins for histogram
+            nBins           100;
+            // Max droplet diameter
+            maxDiameter     0.5e-4;
+            //// Min droplet diameter (default is 0)
+            //minDiameter     0;
+
+            // Writing format
+            setFormat       gnuplot;
+        }
+    }
+
+
 
 SourceFiles
     regionSizeDistribution.C
@@ -41,6 +97,9 @@ SourceFiles
 
 #include "pointFieldFwd.H"
 #include "writer.H"
+#include "Map.H"
+#include "volFieldsFwd.H"
+#include "wordReList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,6 +110,8 @@ namespace Foam
 class objectRegistry;
 class dictionary;
 class mapPolyMesh;
+class regionSplit;
+class polyMesh;
 
 /*---------------------------------------------------------------------------*\
                            Class regionSizeDistribution Declaration
@@ -72,23 +133,85 @@ class regionSizeDistribution
         word alphaName_;
 
         //- Patches to walk from
-        wordList patchNames_;
+        wordReList patchNames_;
 
         //- Clip value
         scalar threshold_;
 
-        //- Background region volFraction
-        scalar volFraction_;
+        //- Maximum droplet diameter
+        scalar maxDiam_;
+
+        //- Minimum droplet diameter
+        scalar minDiam_;
 
         //- Mumber of bins
         label nBins_;
 
+        //- Names of fields to sample on regions
+        wordReList fields_;
+
         //- Output formatter to write
         autoPtr<writer<scalar> > formatterPtr_;
 
 
     // Private Member Functions
 
+        template<class Type>
+        Map<Type> regionSum(const regionSplit&, const Field<Type>&) const;
+
+        //- Get data in order
+        template<class Type>
+        List<Type> extractData(const UList<label>& keys, const Map<Type>&)
+        const;
+
+        void writeGraph
+        (
+            const coordSet& coords,
+            const word& valueName,
+            const scalarField& values
+        ) const;
+
+        //- Write volfields with the parts of alpha which are not
+        //  droplets (liquidCore, backGround)
+        void writeAlphaFields
+        (
+            const regionSplit& regions,
+            const Map<label>& keepRegions,
+            const Map<scalar>& regionVolume,
+            const volScalarField& alpha
+        ) const;
+
+        //- Mark all regions starting at patches
+        Map<label> findPatchRegions(const polyMesh&, const regionSplit&) const;
+
+        //- Helper: divide if denom != 0
+        static tmp<scalarField> divide(const scalarField&, const scalarField&);
+
+        //- Given per-region data calculate per-bin average/deviation and graph
+        void writeGraphs
+        (
+            const word& fieldName,              // name of field
+            const labelList& indices,           // index of bin for each region
+            const scalarField& sortedField,     // per region field data
+            const scalarField& binCount,        // per bin number of regions
+            const coordSet& coords              // graph data for bins
+        ) const;
+
+        //- Given per-cell data calculate per-bin average/deviation and graph
+        void writeGraphs
+        (
+            const word& fieldName,              // name of field
+            const scalarField& cellField,       // per cell field data
+
+            const regionSplit& regions,         // per cell the region(=droplet)
+            const labelList& sortedRegions,     // valid regions in sorted order
+            const scalarField& sortedNormalisation,
+
+            const labelList& indices,           // index of bin for each region
+            const scalarField& binCount,        // per bin number of regions
+            const coordSet& coords              // graph data for bins
+        ) const;
+
         //- Disallow default bitwise copy construct
         regionSizeDistribution(const regionSizeDistribution&);
 
@@ -156,6 +279,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "regionSizeDistributionTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionTemplates.C b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionTemplates.C
new file mode 100644
index 0000000000..2f55434315
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionTemplates.C
@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "regionSizeDistribution.H"
+#include "regionSplit.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Map<Type> Foam::regionSizeDistribution::regionSum
+(
+    const regionSplit& regions,
+    const Field<Type>& fld
+) const
+{
+    // Per region the sum of fld
+    Map<Type> regionToSum(regions.nRegions()/Pstream::nProcs());
+
+    forAll(fld, cellI)
+    {
+        label regionI = regions[cellI];
+
+        typename Map<Type>::iterator fnd = regionToSum.find(regionI);
+        if (fnd == regionToSum.end())
+        {
+            regionToSum.insert(regionI, fld[cellI]);
+        }
+        else
+        {
+            fnd() += fld[cellI];
+        }
+    }
+    Pstream::mapCombineGather(regionToSum, plusEqOp<Type>());
+    Pstream::mapCombineScatter(regionToSum);
+
+    return regionToSum;
+}
+
+
+// Get data in sortedToc order
+template<class Type>
+Foam::List<Type> Foam::regionSizeDistribution::extractData
+(
+    const UList<label>& keys,
+    const Map<Type>& regionData
+) const
+{
+    List<Type> sortedData(keys.size());
+
+    forAll(keys, i)
+    {
+        sortedData[i] = regionData[keys[i]];
+    }
+    return sortedData;
+}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/forces/forces/forces.C b/src/postProcessing/functionObjects/forces/forces/forces.C
index 1dd26906b7..ff9ba3a4e3 100644
--- a/src/postProcessing/functionObjects/forces/forces/forces.C
+++ b/src/postProcessing/functionObjects/forces/forces/forces.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,14 +66,14 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
         const compressible::LESModel& les =
         obr_.lookupObject<compressible::LESModel>("LESProperties");
 
-        return les.devRhoBeff();
+        return les.devRhoReff();
     }
     else if (obr_.foundObject<incompressible::LESModel>("LESProperties"))
     {
         const incompressible::LESModel& les
             = obr_.lookupObject<incompressible::LESModel>("LESProperties");
 
-        return rho()*les.devBeff();
+        return rho()*les.devReff();
     }
     else if (obr_.foundObject<basicThermo>("thermophysicalProperties"))
     {
diff --git a/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C b/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C
index a512419c50..5e6c38045b 100644
--- a/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C
+++ b/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -170,6 +170,7 @@ namespace Foam
         {
             const symmTensor& v = values[elemI];
             os  << float(v[0]) << ' ' << float(v[1]) << ' ' << float(v[2])
+                << ' '
                 << float(v[3]) << ' ' << float(v[4]) << ' ' << float(v[5])
                 << nl;
 
@@ -190,7 +191,9 @@ namespace Foam
         {
             const tensor& v = values[elemI];
             os  << float(v[0]) << ' ' << float(v[1]) << ' ' << float(v[2])
+                << ' '
                 << float(v[3]) << ' ' << float(v[4]) << ' ' << float(v[5])
+                << ' '
                 << float(v[6]) << ' ' << float(v[7]) << ' ' << float(v[8])
                 << nl;
         }
diff --git a/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C b/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C
index 642d2aeee8..1b1d55c8a7 100644
--- a/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C
@@ -128,10 +128,14 @@ greyDiffusiveViewFactorFixedValueFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
 void Foam::radiation::greyDiffusiveViewFactorFixedValueFvPatchScalarField::
 updateCoeffs()
 {
+    if (this->updated())
+    {
+        return;
+    }
+
     // Do nothing
 
     if (debug)
@@ -149,6 +153,7 @@ updateCoeffs()
             << endl;
     }
 
+    fixedValueFvPatchScalarField::updateCoeffs();
 }
 
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/convectiveHeatTransfer/convectiveHeatTransferFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/convectiveHeatTransfer/convectiveHeatTransferFvPatchScalarField.C
index c747a97b84..116bfc9c25 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/convectiveHeatTransfer/convectiveHeatTransferFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/convectiveHeatTransfer/convectiveHeatTransferFvPatchScalarField.C
@@ -142,6 +142,8 @@ void convectiveHeatTransferFvPatchScalarField::updateCoeffs()
             htc[faceI] = 0.037*pow(Re, 0.8)*cbrt(Pr[faceI])*kappaw[faceI]/L_;
         }
     }
+
+    fixedValueFvPatchScalarField::updateCoeffs();
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C
index 0c905a91ce..69795c56d7 100644
--- a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C
+++ b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,17 +87,34 @@ tmp<volSymmTensorField> GenEddyVisc::B() const
 }
 
 
-tmp<volSymmTensorField> GenEddyVisc::devBeff() const
+tmp<volSymmTensorField> GenEddyVisc::devReff() const
 {
     return -nuEff()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> GenEddyVisc::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> GenEddyVisc::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nuEff(), U)
+      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> GenEddyVisc::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H
index 0746377eb7..62c216ee3e 100644
--- a/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H
+++ b/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -118,11 +118,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
         //- Return the deviatoric part of the effective sub-grid
         //  turbulence stress tensor including the laminar stress
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct Eddy-Viscosity and related properties
         virtual void correct(const tmp<volTensorField>& gradU);
diff --git a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.C b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.C
index 9d385c21eb..1618fb7f3b 100644
--- a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.C
+++ b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,7 +115,7 @@ GenSGSStress::GenSGSStress
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-tmp<volSymmTensorField> GenSGSStress::devBeff() const
+tmp<volSymmTensorField> GenSGSStress::devReff() const
 {
     return tmp<volSymmTensorField>
     (
@@ -135,7 +135,7 @@ tmp<volSymmTensorField> GenSGSStress::devBeff() const
 }
 
 
-tmp<fvVectorMatrix> GenSGSStress::divDevBeff
+tmp<fvVectorMatrix> GenSGSStress::divDevReff
 (
     volVectorField& U
 ) const
@@ -164,6 +164,38 @@ tmp<fvVectorMatrix> GenSGSStress::divDevBeff
 }
 
 
+tmp<fvVectorMatrix> GenSGSStress::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    if (couplingFactor_.value() > 0.0)
+    {
+        return
+        (
+            fvc::div(rho*B_ + couplingFactor_*rho*nuSgs_*fvc::grad(U))
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*rho*nuSgs_, U, "laplacian(muEff,U)"
+            )
+          - fvm::laplacian(muEff, U)
+        );
+    }
+    else
+    {
+        return
+        (
+            fvc::div(rho*B_)
+          + fvc::laplacian(rho*nuSgs_, U, "laplacian(muEff,U)")
+          - fvm::laplacian(muEff, U)
+        );
+    }
+}
+
+
 bool GenSGSStress::read()
 {
     if (LESModel::read())
diff --git a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H
index 36d84f9c77..0456d69a71 100644
--- a/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H
+++ b/src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -128,11 +128,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Returns div(B).
-        // This is the additional term due to the filtering of the NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Read LESProperties dictionary
         virtual bool read();
diff --git a/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H b/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H
index deca26e793..3688cb8dec 100644
--- a/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H
+++ b/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H
@@ -202,14 +202,6 @@ public:
         //- Return the sub-grid stress tensor.
         virtual tmp<volSymmTensorField> B() const = 0;
 
-        //- Return the deviatoric part of the effective sub-grid
-        //  turbulence stress tensor including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const = 0;
-
-        //- Returns div(dev(Beff)).
-        //  This is the additional term due to the filtering of the NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const = 0;
-
 
         // RAS compatibility functions for the turbulenceModel base class
 
@@ -225,18 +217,6 @@ public:
                 return B();
             }
 
-            //- Return the effective stress tensor including the laminar stress
-            virtual tmp<volSymmTensorField> devReff() const
-            {
-                return devBeff();
-            }
-
-            //- Return the source term for the momentum equation
-            virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const
-            {
-                return divDevBeff(U);
-            }
-
 
         //- Correct Eddy-Viscosity and related properties.
         //  This calls correct(const tmp<volTensorField>& gradU) by supplying
diff --git a/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.C b/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.C
index 97d32bce9d..84fc1fc7be 100644
--- a/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.C
+++ b/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,7 +81,7 @@ tmp<volSymmTensorField> Smagorinsky2::B() const
 }
 
 
-tmp<fvVectorMatrix> Smagorinsky2::divDevBeff
+tmp<fvVectorMatrix> Smagorinsky2::divDevReff
 (
     volVectorField& U
 ) const
@@ -101,6 +101,28 @@ tmp<fvVectorMatrix> Smagorinsky2::divDevBeff
 }
 
 
+tmp<fvVectorMatrix> Smagorinsky2::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volTensorField gradU(fvc::grad(U));
+
+    volSymmTensorField aniMuEff
+    (
+        "muEff",
+        I*(rho*nuEff()) + (cD2_*rho*delta())*symm(gradU)
+    );
+
+    return
+    (
+      - fvm::laplacian(aniMuEff, U)
+      - fvc::div(rho*nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool Smagorinsky2::read()
 {
     if (Smagorinsky::read())
diff --git a/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.H b/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.H
index 89c77912cd..c9133e4e87 100644
--- a/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.H
+++ b/src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -109,9 +109,17 @@ public:
         //- Return B.
         virtual tmp<volSymmTensorField> B() const;
 
-        //- Returns div(B).
-        // This is the additional term due to the filtering of the NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Read LESProperties dictionary
         virtual bool read();
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
index 6a34f90892..da8ab55e07 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -338,17 +338,34 @@ tmp<volSymmTensorField> SpalartAllmaras::B() const
 }
 
 
-tmp<volSymmTensorField> SpalartAllmaras::devBeff() const
+tmp<volSymmTensorField> SpalartAllmaras::devReff() const
 {
     return -nuEff()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> SpalartAllmaras::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> SpalartAllmaras::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nuEff(), U)
+      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
index 7ed507712a..68e5ef4c73 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -172,11 +172,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Return the deviatoric part of the divergence of Beff
-        //  i.e. the additional term in the filtered NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct nuTilda and related properties
         virtual void correct(const tmp<volTensorField>& gradU);
diff --git a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
index 5e6890d016..12a6f1b3d1 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -388,7 +388,7 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
           - fvm::Sp(fvc::div(phi()), omega_)
           - fvm::laplacian(DomegaEff(F1), omega_)
         ==
-            gamma(F1)*0.5*S2
+            gamma(F1)*S2
           - fvm::Sp(beta(F1)*omega_, omega_)
           - fvm::SuSp       // cross diffusion term
             (
@@ -425,17 +425,34 @@ tmp<volSymmTensorField> kOmegaSSTSAS::B() const
 }
 
 
-tmp<volSymmTensorField> kOmegaSSTSAS::devBeff() const
+tmp<volSymmTensorField> kOmegaSSTSAS::devReff() const
 {
     return -nuEff()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> kOmegaSSTSAS::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> kOmegaSSTSAS::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nuEff(), U)
+      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> kOmegaSSTSAS::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
index 2ca98380af..e61156f9c9 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -249,11 +249,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Return the deviatoric part of the divergence of Beff
-        //  i.e. the additional term in the filtered NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Solve the turbulence equations (k-w) and correct the turbulence
         //  viscosity
diff --git a/src/turbulenceModels/incompressible/LES/laminar/laminar.C b/src/turbulenceModels/incompressible/LES/laminar/laminar.C
index c88ff32beb..5a78c213ee 100644
--- a/src/turbulenceModels/incompressible/LES/laminar/laminar.C
+++ b/src/turbulenceModels/incompressible/LES/laminar/laminar.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -136,17 +136,34 @@ tmp<volSymmTensorField> laminar::B() const
 }
 
 
-tmp<volSymmTensorField> laminar::devBeff() const
+tmp<volSymmTensorField> laminar::devReff() const
 {
     return -nu()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> laminar::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> laminar::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nu(), U) - fvc::div(nu()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nu(), U)
+      - fvc::div(nu()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> laminar::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/laminar/laminar.H b/src/turbulenceModels/incompressible/LES/laminar/laminar.H
index 14de7a142a..777a9228a7 100644
--- a/src/turbulenceModels/incompressible/LES/laminar/laminar.H
+++ b/src/turbulenceModels/incompressible/LES/laminar/laminar.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,13 +104,21 @@ public:
         //- Return the sub-grid stress tensor B.
         virtual tmp<volSymmTensorField> B() const;
 
+        //- Return the effective sub-grid turbulence stress tensor
+        //  including the laminar stress
+        virtual tmp<volSymmTensorField> devReff() const;
+
         //- Return the deviatoric part of the effective sub-grid
         //  turbulence stress tensor including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
-        //- Return the deviatoric part of the divergence of Beff
-        //  i.e. the additional term in the filtered NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Read LESProperties dictionary
         bool read();
diff --git a/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.C b/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.C
index 4f2804577f..b17f4ea572 100644
--- a/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.C
+++ b/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,25 +91,39 @@ tmp<volSymmTensorField> mixedSmagorinsky::B() const
 }
 
 
-tmp<volSymmTensorField> mixedSmagorinsky::devBeff() const
+tmp<volSymmTensorField> mixedSmagorinsky::devReff() const
 {
     return
     (
-        scaleSimilarity::devBeff()
-      + Smagorinsky::devBeff()
+        scaleSimilarity::devReff()
+      + Smagorinsky::devReff()
     );
 }
 
 
-tmp<fvVectorMatrix> mixedSmagorinsky::divDevBeff
+tmp<fvVectorMatrix> mixedSmagorinsky::divDevReff
 (
     volVectorField& U
 ) const
 {
     return
     (
-        scaleSimilarity::divDevBeff(U)
-      + Smagorinsky::divDevBeff(U)
+        scaleSimilarity::divDevReff(U)
+      + Smagorinsky::divDevReff(U)
+    );
+}
+
+
+tmp<fvVectorMatrix> mixedSmagorinsky::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    return
+    (
+        scaleSimilarity::divDevRhoReff(rho, U)
+      + Smagorinsky::divDevRhoReff(rho, U)
     );
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.H b/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.H
index c8557db239..72dca5d354 100644
--- a/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.H
+++ b/src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -129,11 +129,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Implementation of div(B). This is necessary to override
-        // (and include) the div(B) terms from both the parent classes.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct Eddy-Viscosity and related properties
         virtual void correct(const tmp<volTensorField>& gradU);
diff --git a/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.C b/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.C
index b253573be1..61322de356 100644
--- a/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.C
+++ b/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,15 +79,25 @@ tmp<volSymmTensorField> scaleSimilarity::B() const
 }
 
 
-tmp<volSymmTensorField> scaleSimilarity::devBeff() const
+tmp<volSymmTensorField> scaleSimilarity::devReff() const
 {
     return dev(B());
 }
 
 
-tmp<fvVectorMatrix> scaleSimilarity::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> scaleSimilarity::divDevReff(volVectorField& U) const
 {
-    return fvm::Su(fvc::div(devBeff()), U);
+    return fvm::Su(fvc::div(devReff()), U);
+}
+
+
+tmp<fvVectorMatrix> scaleSimilarity::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    return fvm::Su(fvc::div(rho*devReff()), U);
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.H b/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.H
index 0dff31d87e..c757a8165e 100644
--- a/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.H
+++ b/src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -106,13 +106,21 @@ public:
         //- Return the sub-grid stress tensor.
         virtual tmp<volSymmTensorField> B() const;
 
+        //- Return the effective sub-grid turbulence stress tensor
+        //  including the laminar stress
+        virtual tmp<volSymmTensorField> devReff() const;
+
         //- Return the deviatoric part of the effective sub-grid
         //  turbulence stress tensor including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
-        //- Return the deviatoric part of the divergence of Beff
-        //  i.e. the additional term in the filtered NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct Eddy-Viscosity and related properties
         virtual void correct(const tmp<volTensorField>&);
diff --git a/src/turbulenceModels/incompressible/RAS/LRR/LRR.C b/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
index 30afa63ce9..c2c7e5efbe 100644
--- a/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
+++ b/src/turbulenceModels/incompressible/RAS/LRR/LRR.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -259,6 +259,44 @@ tmp<fvVectorMatrix> LRR::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LRR::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    if (couplingFactor_.value() > 0.0)
+    {
+        return
+        (
+            fvc::div
+            (
+                rho*R_ + couplingFactor_*(rho*nut_)*fvc::grad(U),
+                "div((rho*R))"
+            )
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*rho*nut_,
+                U,
+                "laplacian(muEff,U)"
+            )
+          - fvm::laplacian(muEff, U)
+        );
+    }
+    else
+    {
+        return
+        (
+            fvc::div(rho*R_)
+          + fvc::laplacian(rho*nut_, U, "laplacian(muEff,U)")
+          - fvm::laplacian(muEff, U)
+        );
+    }
+}
+
+
 bool LRR::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LRR/LRR.H b/src/turbulenceModels/incompressible/RAS/LRR/LRR.H
index 2f1a119b9c..e787356008 100644
--- a/src/turbulenceModels/incompressible/RAS/LRR/LRR.H
+++ b/src/turbulenceModels/incompressible/RAS/LRR/LRR.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -173,6 +173,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
index 1671c625a0..10de46f587 100644
--- a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -204,6 +204,22 @@ tmp<fvVectorMatrix> LamBremhorstKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LamBremhorstKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LamBremhorstKE::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.H b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.H
index 449055f4c5..3ae7cfda29 100644
--- a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.H
+++ b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -147,6 +147,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
index aafa62af19..49326b1ff5 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -276,7 +276,12 @@ tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevReff(volVectorField& U) const
         return
         (
             fvc::div(R_ + couplingFactor_*nut_*fvc::grad(U), "div(R)")
-          + fvc::laplacian((1.0-couplingFactor_)*nut_, U, "laplacian(nuEff,U)")
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*nut_,
+                U,
+                "laplacian(nuEff,U)"
+            )
           - fvm::laplacian(nuEff(), U)
         );
     }
@@ -292,6 +297,44 @@ tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    if (couplingFactor_.value() > 0.0)
+    {
+        return
+        (
+            fvc::div
+            (
+                rho*R_ + couplingFactor_*(rho*nut_)*fvc::grad(U),
+                "div((rho*R))"
+            )
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*rho*nut_,
+                U,
+                "laplacian(muEff,U)"
+            )
+          - fvm::laplacian(muEff, U)
+        );
+    }
+    else
+    {
+        return
+        (
+            fvc::div(rho*R_)
+          + fvc::laplacian(rho*nut_, U, "laplacian(muEff,U)")
+          - fvm::laplacian(muEff, U)
+        );
+    }
+}
+
+
 bool LaunderGibsonRSTM::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.H b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.H
index 0b6fc91b4b..528a7dd5d6 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.H
+++ b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -184,6 +184,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
index f0fa6e60c8..c514a378e2 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -210,6 +210,22 @@ tmp<fvVectorMatrix> LaunderSharmaKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LaunderSharmaKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LaunderSharmaKE::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.H b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.H
index 3f2338e18e..5e4ae36f5a 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.H
+++ b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -164,6 +164,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
index 03e1e5aa1a..20d8dee951 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
@@ -298,6 +298,23 @@ tmp<fvVectorMatrix> LienCubicKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LienCubicKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+        fvc::div(rho*nonlinearStress_)
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LienCubicKE::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.H b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.H
index ebe7c373d8..56696f4253 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.H
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.H
@@ -159,6 +159,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
index e0bd0f0df9..da81b66edd 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
@@ -365,6 +365,23 @@ tmp<fvVectorMatrix> LienCubicKELowRe::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LienCubicKELowRe::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+        fvc::div(rho*nonlinearStress_)
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LienCubicKELowRe::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.H b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.H
index 0efd63a454..8794af1fc4 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.H
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.H
@@ -186,6 +186,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
index 91e288f104..15296cdf5b 100644
--- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -243,12 +243,27 @@ tmp<fvVectorMatrix> LienLeschzinerLowRe::divDevReff(volVectorField& U) const
     return
     (
       - fvm::laplacian(nuEff(), U)
-    //- (fvc::grad(U) & fvc::grad(nuEff()))
       - fvc::div(nuEff()*T(fvc::grad(U)))
     );
 }
 
 
+tmp<fvVectorMatrix> LienLeschzinerLowRe::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LienLeschzinerLowRe::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.H b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.H
index 3fa1c790fe..5aebe9ea50 100644
--- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.H
+++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -156,6 +156,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
index 33cc03edac..4dbb7b7b68 100644
--- a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
+++ b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
@@ -287,6 +287,23 @@ tmp<fvVectorMatrix> NonlinearKEShih::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> NonlinearKEShih::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+        fvc::div(rho*nonlinearStress_)
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool NonlinearKEShih::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.H b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.H
index c7ad8a3713..f63f0c6565 100644
--- a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.H
+++ b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.H
@@ -162,6 +162,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
index 794d9f9b07..460a6609ac 100644
--- a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -220,6 +220,22 @@ tmp<fvVectorMatrix> RNGkEpsilon::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> RNGkEpsilon::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool RNGkEpsilon::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.H b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.H
index 62b60b9804..2d06ee9ec1 100644
--- a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.H
+++ b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -161,6 +161,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
index 1804ba5b01..4105cb1876 100644
--- a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
@@ -337,6 +337,22 @@ tmp<fvVectorMatrix> SpalartAllmaras::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool SpalartAllmaras::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H
index fbd1248181..4190fc7d80 100644
--- a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -184,6 +184,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
index 51f4aae20c..80aa64362a 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
@@ -116,7 +116,8 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
 
     z_ /= mag(z_);
 
-    evaluate();
+    const vectorField& c = patch().Cf();
+    scalarField::operator=(pow3(Ustar_)/(kappa_*((c & z_) - zGround_ + z0_)));
 }
 
 
@@ -169,14 +170,6 @@ void atmBoundaryLayerInletEpsilonFvPatchScalarField::rmap
 }
 
 
-void atmBoundaryLayerInletEpsilonFvPatchScalarField::updateCoeffs()
-{
-    const vectorField& c = patch().Cf();
-    tmp<scalarField> coord = (c & z_);
-    scalarField::operator=(pow3(Ustar_)/(kappa_*(coord - zGround_ + z0_)));
-}
-
-
 void atmBoundaryLayerInletEpsilonFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
index c72e35372a..bce9e01480 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
@@ -176,7 +176,7 @@ public:
         // Access
 
             //- Return max value
-            scalarField Ustar() const
+            const scalarField& Ustar() const
             {
                 return Ustar_;
             }
@@ -204,12 +204,6 @@ public:
             );
 
 
-        // Evaluation functions
-
-            //- Update coefficients
-            virtual void updateCoeffs();
-
-
         //- Write
         virtual void write(Ostream&) const;
 };
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
index 3a33677025..e41c3904fa 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
@@ -120,7 +120,25 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
         Ustar_[i] = kappa_*Uref_/(log((Href_  + z0_[i])/max(z0_[i] , 0.001)));
     }
 
-    evaluate();
+    const vectorField& c = patch().Cf();
+    const scalarField coord(c & z_);
+    scalarField Un(coord.size());
+
+    forAll(coord, i)
+    {
+        if ((coord[i] - zGround_[i]) < Href_)
+        {
+            Un[i] =
+                (Ustar_[i]/kappa_)
+              * log((coord[i] - zGround_[i] + z0_[i])/max(z0_[i], 0.001));
+        }
+        else
+        {
+            Un[i] = Uref_;
+        }
+    }
+
+    vectorField::operator=(n_*Un);
 }
 
 
@@ -174,32 +192,6 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::rmap
 }
 
 
-void atmBoundaryLayerInletVelocityFvPatchVectorField::updateCoeffs()
-{
-    const vectorField& c = patch().Cf();
-    const scalarField coord(c & z_);
-    scalarField Un(coord.size());
-
-    forAll(coord, i)
-    {
-        if ((coord[i] - zGround_[i]) < Href_)
-        {
-            Un[i] =
-                (Ustar_[i]/kappa_)
-              * log((coord[i] - zGround_[i] + z0_[i])/max(z0_[i], 0.001));
-        }
-        else
-        {
-            Un[i] = Uref_;
-        }
-    }
-
-    vectorField::operator=(n_*Un);
-
-    fixedValueFvPatchVectorField::updateCoeffs();
-}
-
-
 void atmBoundaryLayerInletVelocityFvPatchVectorField::write(Ostream& os) const
 {
     fvPatchVectorField::write(os);
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
index 68a6408f2b..0249b7111c 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
@@ -180,23 +180,26 @@ public:
 
     // Member functions
 
-        //- Return Ustar
-        scalarField& Ustar()
-        {
-            return Ustar_;
-        }
+        // Access
 
-        //- Return flow direction
-        vector& n()
-        {
-            return n_;
-        }
+            //- Return Ustar
+            const scalarField& Ustar() const
+            {
+                return Ustar_;
+            }
+
+            //- Return flow direction
+            const vector& n() const
+            {
+                return n_;
+            }
+
+            //- Return z direction
+            const vector& z() const
+            {
+                return z_;
+            }
 
-        //- Return z direction
-        vector& z()
-        {
-            return z_;
-        }
 
         // Mapping functions
 
@@ -214,12 +217,6 @@ public:
             );
 
 
-        // Evaluation functions
-
-            //- Update coefficients
-            virtual void updateCoeffs();
-
-
         //- Write
         virtual void write(Ostream&) const;
 };
diff --git a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
index 99b4188cc0..029b9f9db7 100644
--- a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -192,6 +192,22 @@ tmp<fvVectorMatrix> kEpsilon::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kEpsilon::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kEpsilon::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.H b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.H
index a3baf20513..89c94fe42f 100644
--- a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.H
+++ b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -155,6 +155,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
index 7106c0ea2c..51e1fbd491 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -201,6 +201,22 @@ tmp<fvVectorMatrix> kOmega::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kOmega::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kOmega::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
index 10729f334e..1e63fd2eca 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
+++ b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -188,6 +188,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
index c025076baa..e6c49db1c8 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -308,6 +308,22 @@ tmp<fvVectorMatrix> kOmegaSST::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kOmegaSST::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kOmegaSST::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
index bd4b3742b5..60f2775b1b 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -258,6 +258,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
index abf9b29cde..45f4f69d13 100644
--- a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
+++ b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C
@@ -570,6 +570,22 @@ tmp<fvVectorMatrix> kkLOmega::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kkLOmega::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kkLOmega::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H
index 83db5f82d5..dbd5ce0832 100644
--- a/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H
+++ b/src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -276,6 +276,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/laminar/laminar.C b/src/turbulenceModels/incompressible/RAS/laminar/laminar.C
index 68adca77b3..f0597b58da 100644
--- a/src/turbulenceModels/incompressible/RAS/laminar/laminar.C
+++ b/src/turbulenceModels/incompressible/RAS/laminar/laminar.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -177,6 +177,22 @@ tmp<fvVectorMatrix> laminar::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> laminar::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool laminar::read()
 {
     return RASModel::read();
diff --git a/src/turbulenceModels/incompressible/RAS/laminar/laminar.H b/src/turbulenceModels/incompressible/RAS/laminar/laminar.H
index 90afadff13..ccebe774ed 100644
--- a/src/turbulenceModels/incompressible/RAS/laminar/laminar.H
+++ b/src/turbulenceModels/incompressible/RAS/laminar/laminar.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,6 +105,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Correct the laminar viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
index 25940abb9e..1fd761b1d7 100644
--- a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
+++ b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -262,6 +262,22 @@ tmp<fvVectorMatrix> qZeta::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> qZeta::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool qZeta::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H
index ba8dc1b413..9347998a2e 100644
--- a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H
+++ b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -211,6 +211,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
index 71058caa7d..495b8930df 100644
--- a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
+++ b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -246,6 +246,22 @@ tmp<fvVectorMatrix> realizableKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> realizableKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool realizableKE::read()
 {
     if (RASModel::read())
diff --git a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.H b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.H
index 045f333390..ed2bc4841e 100644
--- a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.H
+++ b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -180,6 +180,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.C b/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.C
index c102b52cb1..a633c7a70e 100644
--- a/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.C
+++ b/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -203,6 +203,22 @@ tmp<fvVectorMatrix> laminar::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> laminar::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 void laminar::correct()
 {
     turbulenceModel::correct();
diff --git a/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.H b/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.H
index 72994e90fb..5b79ecd485 100644
--- a/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.H
+++ b/src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -111,6 +111,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Correct the laminar viscosity
         virtual void correct();
 
diff --git a/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H b/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H
index 405cfef0bf..c2eacdaf39 100644
--- a/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H
+++ b/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H
@@ -204,6 +204,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const = 0;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const = 0;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct() = 0;
 
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/alphaSgs b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/alphaSgs
index f50fb4a481..fbf52cc809 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/alphaSgs
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/alphaSgs
@@ -24,6 +24,7 @@ boundaryField
     ground
     {
         type        alphatWallFunction;
+        mut         muSgs;
         value       uniform 0;
     }
 
@@ -45,6 +46,7 @@ boundaryField
     "(region0_to.*)"
     {
         type        alphatWallFunction;
+        mut         muSgs;
         value       uniform 0;
     }
 }
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/combustionProperties b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/combustionProperties
index 08ddd33fd3..355d1e1b36 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/combustionProperties
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/combustionProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  infinitelyFastChemistry<psiCombustionModel,gasThermoPhysics>;
+combustionModel  infinitelyFastChemistry<psiThermoCombustion,gasThermoPhysics>;
 
 active true;
 
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/panelRegion/fvSolution b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/panelRegion/fvSolution
index 688f40489e..e1d5b8db0b 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/panelRegion/fvSolution
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/panelRegion/fvSolution
@@ -32,7 +32,7 @@ solvers
         relTol          0;
     }
 
-    "rho|rhot"
+    rhoThermo
     {
         solver          PCG;
         preconditioner  DIC;
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties
index ae740ee671..9973e9ad7b 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties
@@ -15,9 +15,9 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//combustionModel  noCombustion<psiCombustionModel>;
-combustionModel  infinitelyFastChemistry<psiCombustionModel,gasThermoPhysics>;
-//combustionModel  FSD<psiCombustionModel,gasThermoPhysics>;
+//combustionModel  noCombustion<psiThermoCombustion>;
+combustionModel  infinitelyFastChemistry<psiThermoCombustion,gasThermoPhysics>;
+//combustionModel  FSD<psiThermoCombustion,gasThermoPhysics>;
 
 active  true;
 
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/combustionProperties b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/combustionProperties
index 62192a7bd7..cfa2cd8f28 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/combustionProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/combustionProperties
@@ -16,7 +16,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  infinitelyFastChemistry<psiCombustionModel,gasThermoPhysics>;
+combustionModel  infinitelyFastChemistry<psiThermoCombustion,gasThermoPhysics>;
 
 active  on;
 
diff --git a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/polyMesh/blockMeshDict.m4 b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/polyMesh/blockMeshDict.m4
index 79da11e10a..2e865a7856 100644
--- a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/polyMesh/blockMeshDict.m4
+++ b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/polyMesh/blockMeshDict.m4
@@ -107,55 +107,79 @@ edges
 (
 );
 
-patches
+boundary
 (
     // is there no way of defining all my 'defaultFaces' to be 'wall'?
-    wall front
-    (
-    // inlet block
-    frontQuad(in1, join1, join2, in2)
-    // outlet block
-    frontQuad(poro1, out1, out2, poro2)
-    )
+    front
+    {
+        type wall;
+        faces
+        (
+            // inlet block
+            frontQuad(in1, join1, join2, in2)
+            // outlet block
+            frontQuad(poro1, out1, out2, poro2)
+        );
+    }
 
-    wall back
-    (
-    // inlet block
-    backQuad(in1, join1, join2, in2)
-    // outlet block
-    backQuad(poro1, out1, out2, poro2)
-    )
+    back
+    {
+        type wall;
+        faces
+        (
+            // inlet block
+            backQuad(in1, join1, join2, in2)
+            // outlet block
+            backQuad(poro1, out1, out2, poro2)
+        );
+    }
 
-    wall wall
-    (
-    // inlet block
-    quad2D(in1, join1)
-    quad2D(join2, in2)
-    // outlet block
-    quad2D(poro1, out1)
-    quad2D(out2, poro2)
-    )
+    wall
+    {
+        type wall;
+        faces
+        (
+            // inlet block
+            quad2D(in1, join1)
+            quad2D(join2, in2)
+            // outlet block
+            quad2D(poro1, out1)
+            quad2D(out2, poro2)
+        );
+    }
 
-    wall porosityWall
-    (
-    // porosity block
-    frontQuad(join1, poro1, poro2, join2)
-    // porosity block
-    backQuad(join1, poro1, poro2, join2)
-    // porosity block
-    quad2D(join1, poro1)
-    quad2D(poro2, join2)
-    )
+    porosityWall
+    {
+        type wall;
+        faces
+        (
+            // porosity block
+            frontQuad(join1, poro1, poro2, join2)
+            // porosity block
+            backQuad(join1, poro1, poro2, join2)
+            // porosity block
+            quad2D(join1, poro1)
+            quad2D(poro2, join2)
+        );
+    }
 
-    patch inlet
-    (
-    quad2D(in2, in1)
-    )
+    inlet
+    {
+        type patch;
+        faces
+        (
+            quad2D(in2, in1)
+        );
+    }
 
-    patch outlet
-    (
-    quad2D(out2, out1)
-    )
+    outlet
+    {
+        type patch;
+        faces
+        (
+            quad2D(out2, out1)
+        );
+    }
 );
 
 mergePatchPairs
diff --git a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/polyMesh/boundary b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/polyMesh/boundary
deleted file mode 100644
index 0abd1608ab..0000000000
--- a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/constant/polyMesh/boundary
+++ /dev/null
@@ -1,58 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       polyBoundaryMesh;
-    location    "constant/polyMesh";
-    object      boundary;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-6
-(
-    front
-    {
-        type            wall;
-        nFaces          700;
-        startFace       63400;
-    }
-    back
-    {
-        type            wall;
-        nFaces          700;
-        startFace       64100;
-    }
-    wall
-    {
-        type            wall;
-        nFaces          1400;
-        startFace       64800;
-    }
-    porosityWall
-    {
-        type            wall;
-        nFaces          1600;
-        startFace       66200;
-    }
-    inlet
-    {
-        type            patch;
-        nFaces          400;
-        startFace       67800;
-    }
-    outlet
-    {
-        type            patch;
-        nFaces          400;
-        startFace       68200;
-    }
-)
-
-// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
index a664893343..0abd1608ab 100644
--- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
+++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
@@ -8,7 +8,7 @@
 FoamFile
 {
     version     2.0;
-    format      binary;
+    format      ascii;
     class       polyBoundaryMesh;
     location    "constant/polyMesh";
     object      boundary;
diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/thermophysicalProperties b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/thermophysicalProperties
index dd688c1b3f..abc41fb6ce 100644
--- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/thermophysicalProperties
+++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType hePsiThermo<pureMixture<sutherlandTransport<specieThermo<eConstThermo<perfectGas>,sensibleInternalEnergy>>>>;
+thermoType heRhoThermo<pureMixture<sutherlandTransport<specieThermo<hConstThermo<perfectGas>,sensibleInternalEnergy>>>>;
 
 mixture
 {
@@ -26,7 +26,7 @@ mixture
     }
     thermodynamics
     {
-        Cv          719.3;
+        Cp          1005;
         Hf          0;
     }
     transport
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/constant/thermophysicalProperties b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/thermophysicalProperties
index dd688c1b3f..1523354896 100644
--- a/tutorials/compressible/rhoSimplecFoam/squareBend/constant/thermophysicalProperties
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType hePsiThermo<pureMixture<sutherlandTransport<specieThermo<eConstThermo<perfectGas>,sensibleInternalEnergy>>>>;
+thermoType hePsiThermo<pureMixture<sutherlandTransport<specieThermo<hConstThermo<perfectGas>,sensibleInternalEnergy>>>>;
 
 mixture
 {
@@ -26,7 +26,7 @@ mixture
     }
     thermodynamics
     {
-        Cv          719.3;
+        Cp          1007;
         Hf          0;
     }
     transport
diff --git a/tutorials/compressible/sonicFoam/laminar/shockTube/constant/thermophysicalProperties b/tutorials/compressible/sonicFoam/laminar/shockTube/constant/thermophysicalProperties
index 76fd63c9b7..f409900f93 100644
--- a/tutorials/compressible/sonicFoam/laminar/shockTube/constant/thermophysicalProperties
+++ b/tutorials/compressible/sonicFoam/laminar/shockTube/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hePsiThermo<pureMixture<constTransport<specieThermo<eConstThermo<perfectGas>,sensibleInternalEnergy>>>>;
+thermoType      hePsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>,sensibleInternalEnergy>>>>;
 
 mixture
 {
@@ -26,7 +26,7 @@ mixture
     }
     thermodynamics
     {
-        Cv              717.5;
+        Cp              1005;
         Hf              0;
     }
     transport
diff --git a/tutorials/compressible/sonicFoam/ras/nacaAirfoil/constant/thermophysicalProperties b/tutorials/compressible/sonicFoam/ras/nacaAirfoil/constant/thermophysicalProperties
index b302f7e040..e2a0824661 100644
--- a/tutorials/compressible/sonicFoam/ras/nacaAirfoil/constant/thermophysicalProperties
+++ b/tutorials/compressible/sonicFoam/ras/nacaAirfoil/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hePsiThermo<pureMixture<constTransport<specieThermo<eConstThermo<perfectGas>,sensibleInternalEnergy>>>>;
+thermoType      hePsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>,sensibleInternalEnergy>>>>;
 
 mixture
 {
@@ -26,7 +26,7 @@ mixture
     }
     thermodynamics
     {
-        Cv              717.5;
+        Cp              1005;
         Hf              0;
     }
     transport
diff --git a/tutorials/compressible/sonicFoam/ras/prism/constant/thermophysicalProperties b/tutorials/compressible/sonicFoam/ras/prism/constant/thermophysicalProperties
index b302f7e040..e2a0824661 100644
--- a/tutorials/compressible/sonicFoam/ras/prism/constant/thermophysicalProperties
+++ b/tutorials/compressible/sonicFoam/ras/prism/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hePsiThermo<pureMixture<constTransport<specieThermo<eConstThermo<perfectGas>,sensibleInternalEnergy>>>>;
+thermoType      hePsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>,sensibleInternalEnergy>>>>;
 
 mixture
 {
@@ -26,7 +26,7 @@ mixture
     }
     thermodynamics
     {
-        Cv              717.5;
+        Cp              1005;
         Hf              0;
     }
     transport
diff --git a/tutorials/incompressible/pimpleDyMFoam/mixerVesselAMI2D/constant/dynamicMeshDict b/tutorials/incompressible/pimpleDyMFoam/mixerVesselAMI2D/constant/dynamicMeshDict
index 6f91f2d72f..1df8467c2e 100644
--- a/tutorials/incompressible/pimpleDyMFoam/mixerVesselAMI2D/constant/dynamicMeshDict
+++ b/tutorials/incompressible/pimpleDyMFoam/mixerVesselAMI2D/constant/dynamicMeshDict
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                 |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/constant/dynamicMeshDict b/tutorials/incompressible/pimpleDyMFoam/propeller/constant/dynamicMeshDict
index 4ab93d683a..a657902851 100644
--- a/tutorials/incompressible/pimpleDyMFoam/propeller/constant/dynamicMeshDict
+++ b/tutorials/incompressible/pimpleDyMFoam/propeller/constant/dynamicMeshDict
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                 |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/system/changeDictionaryDict b/tutorials/incompressible/pimpleDyMFoam/propeller/system/changeDictionaryDict
deleted file mode 100644
index 3eb586f140..0000000000
--- a/tutorials/incompressible/pimpleDyMFoam/propeller/system/changeDictionaryDict
+++ /dev/null
@@ -1,42 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      changeDictionaryDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dictionaryReplacement
-{
-    boundary
-    {
-        AMI1
-        {
-            type            cyclicAMI;
-            neighbourPatch  AMI2;
-            transform       noOrdering;
-            surface         
-            {
-            }
-        }
-        AMI2
-        {
-            type            cyclicAMI;
-            neighbourPatch  AMI1;
-            transform       noOrdering;
-            surface         
-            {
-            }
-        }
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/Allrun b/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/Allrun
index 5aa430556a..3815c35895 100755
--- a/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/Allrun
+++ b/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/Allrun
@@ -11,13 +11,10 @@ cp $FOAM_TUTORIALS/resources/geometry/motorBike.obj.gz constant/triSurface/
 mkdir 0
 
 runApplication blockMesh
-
 cp system/decomposeParDict.hierarchical system/decomposeParDict
-
 runApplication decomposePar
 
 cp system/decomposeParDict.ptscotch system/decomposeParDict
-
 runParallel snappyHexMesh 8 -overwrite -parallel
 
 find . -type f -iname "*level*" -exec rm {} \;
diff --git a/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/controlDict b/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/controlDict
index d7759e4112..6569a100b7 100644
--- a/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/controlDict
+++ b/tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system/controlDict
@@ -17,7 +17,7 @@ FoamFile
 
 libs            ("libOpenFOAM.so" "libfieldFunctionObjects.so");
 
-application     pisoFoam;
+application     simpleFoam;
 
 startFrom       latestTime;
 
diff --git a/tutorials/incompressible/simpleFoam/motorBike/Allrun b/tutorials/incompressible/simpleFoam/motorBike/Allrun
index cb5a66db11..d8fdb9d007 100755
--- a/tutorials/incompressible/simpleFoam/motorBike/Allrun
+++ b/tutorials/incompressible/simpleFoam/motorBike/Allrun
@@ -1,4 +1,6 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
@@ -12,3 +14,5 @@ runApplication snappyHexMesh -overwrite
 
 runApplication potentialFoam -noFunctionObjects -writep
 runApplication `getApplication`
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0/U b/tutorials/incompressible/simpleFoam/turbineSiting/0.org/U
similarity index 100%
rename from tutorials/incompressible/simpleFoam/turbineSiting/0/U
rename to tutorials/incompressible/simpleFoam/turbineSiting/0.org/U
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0/epsilon b/tutorials/incompressible/simpleFoam/turbineSiting/0.org/epsilon
similarity index 100%
rename from tutorials/incompressible/simpleFoam/turbineSiting/0/epsilon
rename to tutorials/incompressible/simpleFoam/turbineSiting/0.org/epsilon
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0/include/ABLConditions b/tutorials/incompressible/simpleFoam/turbineSiting/0.org/include/ABLConditions
similarity index 100%
rename from tutorials/incompressible/simpleFoam/turbineSiting/0/include/ABLConditions
rename to tutorials/incompressible/simpleFoam/turbineSiting/0.org/include/ABLConditions
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0/include/fixedInlet b/tutorials/incompressible/simpleFoam/turbineSiting/0.org/include/fixedInlet
similarity index 100%
rename from tutorials/incompressible/simpleFoam/turbineSiting/0/include/fixedInlet
rename to tutorials/incompressible/simpleFoam/turbineSiting/0.org/include/fixedInlet
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0/include/initialConditions b/tutorials/incompressible/simpleFoam/turbineSiting/0.org/include/initialConditions
similarity index 100%
rename from tutorials/incompressible/simpleFoam/turbineSiting/0/include/initialConditions
rename to tutorials/incompressible/simpleFoam/turbineSiting/0.org/include/initialConditions
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0/include/sideAndTopPatches b/tutorials/incompressible/simpleFoam/turbineSiting/0.org/include/sideAndTopPatches
similarity index 93%
rename from tutorials/incompressible/simpleFoam/turbineSiting/0/include/sideAndTopPatches
rename to tutorials/incompressible/simpleFoam/turbineSiting/0.org/include/sideAndTopPatches
index 96611f16ff..bddc17121f 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/0/include/sideAndTopPatches
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/0.org/include/sideAndTopPatches
@@ -16,4 +16,9 @@ sides
     type slip;
 }
 
+"proc.*"
+{
+    type            processor;
+}
+
 // ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0/k b/tutorials/incompressible/simpleFoam/turbineSiting/0.org/k
similarity index 100%
rename from tutorials/incompressible/simpleFoam/turbineSiting/0/k
rename to tutorials/incompressible/simpleFoam/turbineSiting/0.org/k
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0/nut b/tutorials/incompressible/simpleFoam/turbineSiting/0.org/nut
similarity index 100%
rename from tutorials/incompressible/simpleFoam/turbineSiting/0/nut
rename to tutorials/incompressible/simpleFoam/turbineSiting/0.org/nut
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0/p b/tutorials/incompressible/simpleFoam/turbineSiting/0.org/p
similarity index 100%
rename from tutorials/incompressible/simpleFoam/turbineSiting/0/p
rename to tutorials/incompressible/simpleFoam/turbineSiting/0.org/p
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/Allclean b/tutorials/incompressible/simpleFoam/turbineSiting/Allclean
index 2e2bcf8132..dad1f81a5f 100755
--- a/tutorials/incompressible/simpleFoam/turbineSiting/Allclean
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/Allclean
@@ -4,22 +4,11 @@ cd ${0%/*} || exit 1    # run from this directory
 # Source tutorial clean functions
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
-cleanCase
-rm -rf VTK
-#rm -rf constant/cellToRegion constant/polyMesh/sets
-rm -rf constant/polyMesh/sets
-#rm -rf constant/cellLevel
-#rm -rf constant/cellZones
-#rm -rf constant/faceZones
-#rm -rf constant/faces
-#rm -rf constant/neighbour
-#rm -rf constant/owner
-#rm -rf constant/pointZones
-#rm -rf constant/points
-#rm -rf constant/refinementHistory
-#rm -rf constant/surfaceIndex
+rm -rf 0 > /dev/null 2>&1
 
-# Reset decomposeParDict
-cp system/decomposeParDict-nonPar system/decomposeParDict
+cleanCase
+
+# Remove decomposeParDict
+rm -f system/decomposeParDict
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/Allrun b/tutorials/incompressible/simpleFoam/turbineSiting/Allrun
index 81d296213c..9852cc1cf9 100755
--- a/tutorials/incompressible/simpleFoam/turbineSiting/Allrun
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/Allrun
@@ -4,36 +4,19 @@ cd ${0%/*} || exit 1    # run from this directory
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
+# Make dummy 0 directory
+mkdir 0
+
 runApplication blockMesh
-cp system/decomposeParDict-nonPar system/decomposeParDict
+cp system/decomposeParDict.hierarchical system/decomposeParDict
 runApplication decomposePar
 
-#runApplication snappyHexMesh -overwrite
-#runApplication setSet -batch makeZones
-#runApplication setsToZones -noFlipMap
-#runApplication `getApplication`
+cp system/decomposeParDict.ptscotch system/decomposeParDict
+runParallel snappyHexMesh 4 -overwrite
 
-cp system/decomposeParDict-par system/decomposeParDict
-runParallel snappyHexMesh 2 -overwrite
-# *ProcAddressing files written by decomposePar no longer valid
-rm -f processor*/constant/polyMesh/*ProcAddressing
+find . -type f -iname "*level*" -exec rm {} \;
 
-# Add wildcard entries for meshed patches since not preserved
-# by decomposePar. Notice -literalRE option to add wildcard itself
-# without evaluation.
-runParallel changeDictionary 2 -literalRE -enableFunctionEntries
-
-cp system/decomposeParDict-4proc system/decomposeParDict
-# Unset floating point trapping since creating processor directories
-unset FOAM_SIGFPE
-unset FOAM_SETNAN
-runParallel redistributePar 4 -overwrite
-runParallel renumberMesh 4 -overwrite
-
-# Add wildcard entries for meshes patches since not preserved
-# by decomposePar. Notice -literalRE option to add wildcard itself
-# without evaluation.
-#runParallel changeDictionary 4 -literalRE
+ls -d processor* | xargs -i cp -r 0.org/* ./{}/0/ $1
 
 runParallel topoSet 4
 runParallel `getApplication` 4
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/constant/polyMesh/boundary b/tutorials/incompressible/simpleFoam/turbineSiting/constant/polyMesh/boundary
index 7b72e85f59..c3da61bd8c 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/constant/polyMesh/boundary
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/constant/polyMesh/boundary
@@ -8,44 +8,50 @@
 FoamFile
 {
     version     2.0;
-    format      ascii;
+    format      binary;
     class       polyBoundaryMesh;
     location    "constant/polyMesh";
     object      boundary;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-5
+6
 (
     outlet
     {
         type            patch;
-        nFaces          600;
-        startFace       51900;
+        nFaces          922;
+        startFace       364825;
     }
     sides
     {
         type            patch;
-        nFaces          1200;
-        startFace       52500;
+        nFaces          1834;
+        startFace       365747;
     }
     inlet
     {
         type            patch;
-        nFaces          600;
-        startFace       53700;
+        nFaces          923;
+        startFace       367581;
     }
     ground
     {
         type            wall;
-        nFaces          900;
-        startFace       54300;
+        nFaces          0;
+        startFace       368504;
     }
     top
     {
         type            patch;
         nFaces          900;
-        startFace       55200;
+        startFace       368504;
+    }
+    terrain_patch0
+    {
+        type            wall;
+        nFaces          14400;
+        startFace       369404;
     }
 )
 
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/system/changeDictionaryDict b/tutorials/incompressible/simpleFoam/turbineSiting/system/changeDictionaryDict
deleted file mode 100644
index f0dc0c34f6..0000000000
--- a/tutorials/incompressible/simpleFoam/turbineSiting/system/changeDictionaryDict
+++ /dev/null
@@ -1,200 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      changeDictionaryDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include        "$FOAM_CASE/0/include/initialConditions"
-#include        "$FOAM_CASE/0/include/ABLConditions"
-
-dictionaryReplacement
-{
-
-    // Specify
-    // - all fvPatchFields with potential non-uniform values
-    // - all fvPatchFields originating from meshing
-    // - all fvPatchFields originating from mesh-redistribution
-
-    p
-    {
-        boundaryField
-        {
-            outlet
-            {
-                type            uniformFixedValue;
-                uniformValue    constant $pressure;
-            }
-            inlet
-            {
-                type            zeroGradient;
-            }
-            "terrain_.*"
-            {
-                type            zeroGradient;
-            }
-            ground
-            {
-                type            zeroGradient;
-            }
-            #include "$FOAM_CASE/0/include/sideAndTopPatches"
-            "procBoundary.*"
-            {
-                type            processor;
-            }
-        }
-    }
-
-    k
-    {
-        boundaryField
-        {
-            outlet
-            {
-                type            inletOutlet;
-                inletValue      uniform 0.0;
-                value           uniform $turbulentKE;
-            }
-            inlet
-            {
-                type            uniformFixedValue;
-                uniformValue    constant $turbulentKE;
-            }
-            "terrain_.*"
-            {
-                type            kqRWallFunction;
-                value           uniform 0.0;
-            }
-            ground
-            {
-                type            zeroGradient;
-            }
-            #include "$FOAM_CASE/0/include/sideAndTopPatches"
-            "procBoundary.*"
-            {
-                type            processor;
-            }
-        }
-    }
-
-    U
-    {
-        boundaryField
-        {
-            outlet
-            {
-                type            inletOutlet;
-                inletValue      uniform (0 0 0);
-                value           uniform $flowVelocity;
-            }
-            inlet
-            {
-                type            atmBoundaryLayerInletVelocity;
-                Uref            $Uref;
-                Href            $Href;
-                n               $windDirection;
-                z               $zDirection;
-                z0              $z0;
-                zGround         $zGround;
-                value           uniform $flowVelocity;
-            }
-            "terrain_.*"
-            {
-                type            uniformFixedValue;
-                uniformValue    constant $flowVelocity;
-            }
-            ground
-            {
-                type            uniformFixedValue;
-                uniformValue    constant $flowVelocity;
-            }
-            #include "$FOAM_CASE/0/include/sideAndTopPatches"
-            "procBoundary.*"
-            {
-                type            processor;
-            }
-        }
-    }
-
-    nut
-    {
-        boundaryField
-        {
-            outlet
-            {
-                type            calculated;
-                value           uniform 0;
-            }
-            inlet
-            {
-                type            calculated;
-                value           uniform 0;
-            }
-            "terrain_.*"
-            {
-                type            nutkAtmRoughWallFunction;
-                z0              $z0;
-                value           uniform 0.0;
-            }
-            ground
-            {
-                type            calculated;
-                value           uniform 0;
-            }
-            #include "$FOAM_CASE/0/include/sideAndTopPatches"
-            "procBoundary.*"
-            {
-                type            processor;
-            }
-        }
-    }
-
-    epsilon
-    {
-        boundaryField
-        {
-            outlet
-            {
-                type            zeroGradient;
-            }
-            inlet
-            {
-                type            atmBoundaryLayerInletEpsilon;
-                z               $zDirection;
-                z0              $z0;
-                zGround         $zGround;
-                Uref            $Uref;
-                Href            $Href;
-                value           uniform $turbulentEpsilon;
-            }
-            "terrain_.*"
-            {
-                type            epsilonWallFunction;
-                Cmu             0.09;
-                kappa           0.4;
-                E               9.8;
-                value           uniform $turbulentEpsilon;
-            }
-            ground
-            {
-                type            zeroGradient;
-            }
-            #include "$FOAM_CASE/0/include/sideAndTopPatches"
-            "procBoundary.*"
-            {
-                type            processor;
-            }
-        }
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict b/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict
deleted file mode 100644
index 42f1606ec9..0000000000
--- a/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict
+++ /dev/null
@@ -1,29 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      decomposeParDict;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-numberOfSubdomains 2;
-
-method          hierarchical;
-
-hierarchicalCoeffs
-{
-    n           (2 1 1);
-    delta       0.001;
-    order       xyz;
-}
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict-par b/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict-par
deleted file mode 100644
index ded6aceec0..0000000000
--- a/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict-par
+++ /dev/null
@@ -1,22 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      decomposeParDict;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-numberOfSubdomains 2;
-
-method          ptscotch;
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict-nonPar b/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict.hierarchical
similarity index 95%
rename from tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict-nonPar
rename to tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict.hierarchical
index 42f1606ec9..541de112df 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict-nonPar
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict.hierarchical
@@ -15,13 +15,13 @@ FoamFile
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-numberOfSubdomains 2;
+numberOfSubdomains 4;
 
 method          hierarchical;
 
 hierarchicalCoeffs
 {
-    n           (2 1 1);
+    n           (2 2 1);
     delta       0.001;
     order       xyz;
 }
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict-4proc b/tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict.ptscotch
similarity index 100%
rename from tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict-4proc
rename to tutorials/incompressible/simpleFoam/turbineSiting/system/decomposeParDict.ptscotch
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/system/fvSchemes b/tutorials/incompressible/simpleFoam/turbineSiting/system/fvSchemes
index e3d9959b1f..4c101301a1 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/system/fvSchemes
@@ -16,20 +16,18 @@ FoamFile
 
 ddtSchemes
 {
-    default            steadyState;
+    default             steadyState;
 }
 
 gradSchemes
 {
-    default            Gauss linear;
-    grad(p)             Gauss linear;
-    grad(U)             Gauss linear;
+    default             Gauss linear;
 }
 
 divSchemes
 {
     default             none;
-    div(phi,U)          Gauss upwind grad(U);
+    div(phi,U)          Gauss upwind;
     div((nuEff*dev(T(grad(U)))))    Gauss linear;
     div(phi,epsilon)    Gauss upwind;
     div(phi,k)          Gauss upwind;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
index 8e76e79670..cd9ebab8b2 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                 |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/patchifyObstacles b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/patchifyObstacles
index 9a580d2a66..bc9fe0a8a3 100755
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/patchifyObstacles
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/patchifyObstacles
@@ -1,11 +1,11 @@
 #!/bin/bash
-#--------------------------------*- C++ -*----------------------------------#
-# =========                 |                                               #
-# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox         #
-#  \\    /   O peration     | Version:  dev                                 #
-#   \\  /    A nd           | Web:      www.OpenFOAM.org                    #
-#    \\/     M anipulation  |                                               #
-#---------------------------------------------------------------------------#
+#--------------------------------*- C++ -*------------------------------------#
+# =========                 |                                                 #
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           #
+#  \\    /   O peration     | Version:  dev                                   #
+#   \\  /    A nd           | Web:      www.OpenFOAM.org                      #
+#    \\/     M anipulation  |                                                 #
+#-----------------------------------------------------------------------------#
 cd ${0%/*} || exit 1    # run from this directory
 
 x0=0.4
diff --git a/tutorials/multiphase/LTSInterFoam/wigleyHull/system/fvSchemes b/tutorials/multiphase/LTSInterFoam/wigleyHull/system/fvSchemes
index 22366f437b..7018c5eaa7 100644
--- a/tutorials/multiphase/LTSInterFoam/wigleyHull/system/fvSchemes
+++ b/tutorials/multiphase/LTSInterFoam/wigleyHull/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(phirb,alpha) Gauss interfaceCompression;
     div(phi,k)      Gauss upwind;
     div(phi,omega)  Gauss upwind;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/MRFInterFoam/mixerVessel2D/system/fvSchemes b/tutorials/multiphase/MRFInterFoam/mixerVessel2D/system/fvSchemes
index 41159c5be6..3827ac6548 100644
--- a/tutorials/multiphase/MRFInterFoam/mixerVessel2D/system/fvSchemes
+++ b/tutorials/multiphase/MRFInterFoam/mixerVessel2D/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss linear;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/fvSchemes b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/fvSchemes
index b4ad45237d..d1ed26cae3 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/fvSchemes
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss limitedLinearV 1;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
index 0c8bf54bb5..3d506aa1d7 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle/system/fvSchemes
@@ -31,6 +31,7 @@ divSchemes
     div(phiv,rho)   Gauss limitedLinear 0.2;
     div(phi,U)      Gauss filteredLinear2V 0.2 0;
     div(phiv,k)     Gauss filteredLinear2 0.2 0;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 gradSchemes
diff --git a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
index 0c8bf54bb5..3d506aa1d7 100644
--- a/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/les/throttle3D/system/fvSchemes
@@ -31,6 +31,7 @@ divSchemes
     div(phiv,rho)   Gauss limitedLinear 0.2;
     div(phi,U)      Gauss filteredLinear2V 0.2 0;
     div(phiv,k)     Gauss filteredLinear2 0.2 0;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 gradSchemes
diff --git a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
index 633e4a749e..32379fb2c2 100644
--- a/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
+++ b/tutorials/multiphase/cavitatingFoam/ras/throttle/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(phi,U)      Gauss limitedLinearV 1;
     div(phiv,omega) Gauss limitedLinear 1;
     div(phiv,k)     Gauss limitedLinear 1;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 gradSchemes
diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes
index c35af6cc4f..903d94d30c 100644
--- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes
+++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phid2,p_rgh) Gauss upwind;
     div(rho*phi,T)  Gauss upwind;
     div(phi,k)      Gauss vanLeer;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes
index c35af6cc4f..903d94d30c 100644
--- a/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes
+++ b/tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phid2,p_rgh) Gauss upwind;
     div(rho*phi,T)  Gauss upwind;
     div(phi,k)      Gauss vanLeer;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSchemes
index fa6ef27c43..ba444728e1 100644
--- a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss upwind;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
index 5ce6b8ee17..852f5e0c74 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(phirb,alpha) Gauss interfaceCompression;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSchemes
index a0a60315bd..b2c0af769d 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSchemes
index a0a60315bd..b2c0af769d 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSchemes
index a0a60315bd..b2c0af769d 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSchemes
index a0a60315bd..b2c0af769d 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D3DoF/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSchemes
index a0a60315bd..b2c0af769d 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank3D6DoF/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interDyMFoam/ras/testTubeMixer/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/testTubeMixer/system/fvSchemes
index a0a60315bd..b2c0af769d 100644
--- a/tutorials/multiphase/interDyMFoam/ras/testTubeMixer/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/testTubeMixer/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss vanLeer;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/laminar/capillaryRise/system/fvSchemes b/tutorials/multiphase/interFoam/laminar/capillaryRise/system/fvSchemes
index c2c563142c..2a580ed9eb 100644
--- a/tutorials/multiphase/interFoam/laminar/capillaryRise/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/laminar/capillaryRise/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss upwind;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSchemes b/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSchemes
index 7bb78c013c..986a20eebe 100644
--- a/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/laminar/damBreak/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss limitedLinearV 1;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSchemes b/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSchemes
index be6806ad73..26e06a656b 100644
--- a/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phi,B)      Gauss limitedLinear 1;
     div(B)          Gauss linear;
     div(phi,nuTilda) Gauss limitedLinear 1;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/ras/damBreak/system/fvSchemes b/tutorials/multiphase/interFoam/ras/damBreak/system/fvSchemes
index 28b3d7d20a..8ed8e03cc0 100644
--- a/tutorials/multiphase/interFoam/ras/damBreak/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/ras/damBreak/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss upwind;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSchemes b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSchemes
index 28b3d7d20a..8ed8e03cc0 100644
--- a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss upwind;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/ras/waterChannel/LTSInterFoam/system/fvSchemes b/tutorials/multiphase/interFoam/ras/waterChannel/LTSInterFoam/system/fvSchemes
index 22366f437b..7018c5eaa7 100644
--- a/tutorials/multiphase/interFoam/ras/waterChannel/LTSInterFoam/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/ras/waterChannel/LTSInterFoam/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(phirb,alpha) Gauss interfaceCompression;
     div(phi,k)      Gauss upwind;
     div(phi,omega)  Gauss upwind;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes
index 43847dd841..bb536c017c 100644
--- a/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/ras/waterChannel/system/fvSchemes
@@ -33,6 +33,7 @@ divSchemes
 
     div(phi,k)      Gauss upwind;
     div(phi,omega)  $div(phi,k);
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interFoam/ras/weirOverflow/system/fvSchemes b/tutorials/multiphase/interFoam/ras/weirOverflow/system/fvSchemes
index 28b3d7d20a..8ed8e03cc0 100644
--- a/tutorials/multiphase/interFoam/ras/weirOverflow/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/ras/weirOverflow/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss upwind;
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes
index ed937573ff..c80ccb314e 100644
--- a/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes
+++ b/tutorials/multiphase/interMixingFoam/laminar/damBreak/system/fvSchemes
@@ -30,6 +30,7 @@ divSchemes
     div(rho*phi,U)  Gauss limitedLinearV 1;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
index c1b7f71206..44313d502d 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
+++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
@@ -30,9 +30,9 @@ divSchemes
     div(rhoPhi,U)        Gauss linearUpwind grad(U);
     div(phi,omega)       Gauss linearUpwind grad(omega);
     div(phi,k)           Gauss linearUpwind grad(k);
-
     div(phi,alpha)       Gauss vanLeer;
     div(phirb,alpha)     Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 gradSchemes
diff --git a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes
index 6e3d6f1b84..e1436932da 100644
--- a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes
+++ b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(rho*phi,U)  Gauss upwind;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSchemes b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSchemes
index 6e3d6f1b84..e1436932da 100644
--- a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSchemes
+++ b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
     div(rho*phi,U)  Gauss upwind;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSchemes b/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSchemes
index 4bb0230588..35930b9203 100644
--- a/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSchemes
+++ b/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     div(rho*phi,U)  Gauss linear;
     div(phi,alpha1) Gauss vanLeer;
     div(phi,k)      Gauss limitedLinear 1;
-    div(((rho*nuEff)*dev(grad(U).T()))) Gauss linear;
+    div((muEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes