diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
index 926abe2a09..d5a5d20631 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H
@@ -178,7 +178,7 @@
     rho2 = rho20 + psi2*p;
 
     K1 = 0.5*magSqr(U1);
-    K2 = 0.5*magSqr(U1);
+    K2 = 0.5*magSqr(U2);
 
     dpdt = fvc::ddt(p);
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/Make/options b/applications/solvers/multiphase/multiphaseEulerFoam/Make/options
index 5799557cc9..7cd8f48ee4 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/Make/options
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/Make/options
@@ -20,4 +20,5 @@ EXE_LIBS = \
     -lincompressibleTransportModels \
     -lcompressibleMultiphaseEulerianInterfacialModels \
     -lincompressibleLESModels \
+    -lincompressibleRASModels \
     -lfiniteVolume
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index 443576601d..4eb19041aa 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -273,73 +273,165 @@ bool Foam::checkCoupledPoints
     const faceList& fcs = mesh.faces();
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    // Zero'th point on coupled faces
-    pointField nbrZeroPoint(fcs.size()-mesh.nInternalFaces(), vector::max);
-
-    // Exchange zero point
-    forAll(patches, patchI)
+    // Check size of faces
+    label maxSize = 0;
     {
-        if (patches[patchI].coupled())
+        labelList nbrSize(fcs.size()-mesh.nInternalFaces(), 0);
+
+        // Exchange size
+        forAll(patches, patchI)
         {
-            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
-            (
-                patches[patchI]
-            );
-
-            forAll(cpp, i)
+            if (patches[patchI].coupled())
             {
-                label bFaceI = cpp.start()+i-mesh.nInternalFaces();
-                const point& p0 = p[cpp[i][0]];
-                nbrZeroPoint[bFaceI] = p0;
-            }
-        }
-    }
-    syncTools::swapBoundaryFacePositions(mesh, nbrZeroPoint);
-
-    // Compare to local ones. Use same tolerance as for matching
-    label nErrorFaces = 0;
-    scalar avgMismatch = 0;
-    label nCoupledFaces = 0;
-
-    forAll(patches, patchI)
-    {
-        if (patches[patchI].coupled())
-        {
-            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
-            (
-                patches[patchI]
-            );
-
-            if (cpp.owner())
-            {
-                scalarField smallDist
+                const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
                 (
-                    cpp.calcFaceTol
-                    (
-                        //cpp.matchTolerance(),
-                        cpp,
-                        cpp.points(),
-                        cpp.faceCentres()
-                    )
+                    patches[patchI]
                 );
 
                 forAll(cpp, i)
                 {
                     label bFaceI = cpp.start()+i-mesh.nInternalFaces();
-                    const point& p0 = p[cpp[i][0]];
+                    nbrSize[bFaceI] = cpp[i].size();
+                    maxSize = max(maxSize, cpp[i].size());
+                }
+            }
+        }
+        syncTools::swapBoundaryFaceList(mesh, nbrSize);
 
-                    scalar d = mag(p0 - nbrZeroPoint[bFaceI]);
 
-                    if (d > smallDist[i])
+        // Check on owner
+        label nErrorFaces = 0;
+        forAll(patches, patchI)
+        {
+            if (patches[patchI].coupled())
+            {
+                const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
+                (
+                    patches[patchI]
+                );
+
+                if (cpp.owner())
+                {
+                    forAll(cpp, i)
                     {
-                        if (setPtr)
+                        label bFaceI = cpp.start()+i-mesh.nInternalFaces();
+
+                        if (cpp[i].size() != nbrSize[bFaceI])
                         {
-                            setPtr->insert(cpp.start()+i);
+                            if (setPtr)
+                            {
+                                setPtr->insert(cpp.start()+i);
+                            }
+                            nErrorFaces++;
+                        }
+                    }
+                }
+            }
+        }
+
+        reduce(nErrorFaces, sumOp<label>());
+        if (nErrorFaces > 0)
+        {
+            if (report)
+            {
+                Info<< "  **Error in coupled faces: "
+                    << nErrorFaces
+                    << " faces have different size "
+                    << " compared to their coupled equivalent." << endl;
+            }
+            return true;
+        }
+
+        reduce(maxSize, maxOp<label>());
+    }
+
+
+    label nErrorFaces = 0;
+    scalar avgMismatch = 0;
+    label nCoupledPoints = 0;
+
+    for (label index = 0; index < maxSize; index++)
+    {
+        // point at index on coupled faces
+        pointField nbrPoint(fcs.size()-mesh.nInternalFaces(), vector::max);
+
+        // Exchange point
+        forAll(patches, patchI)
+        {
+            if (patches[patchI].coupled())
+            {
+                const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
+                (
+                    patches[patchI]
+                );
+
+                forAll(cpp, i)
+                {
+                    const face& f = cpp[i];
+                    if (f.size() > index)
+                    {
+                        label bFaceI = cpp.start()+i-mesh.nInternalFaces();
+                        nbrPoint[bFaceI] = p[f[index]];
+                    }
+                }
+            }
+        }
+        syncTools::swapBoundaryFacePositions(mesh, nbrPoint);
+
+
+        // Compare to local ones. Use same tolerance as for matching
+
+        forAll(patches, patchI)
+        {
+            if (patches[patchI].coupled())
+            {
+                const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
+                (
+                    patches[patchI]
+                );
+
+                if (cpp.owner())
+                {
+                    scalarField smallDist
+                    (
+                        cpp.calcFaceTol
+                        (
+                            //cpp.matchTolerance(),
+                            cpp,
+                            cpp.points(),
+                            cpp.faceCentres()
+                        )
+                    );
+
+                    forAll(cpp, i)
+                    {
+                        const face& f = cpp[i];
+                        if (f.size() > index)
+                        {
+                            label bFaceI = cpp.start()+i-mesh.nInternalFaces();
+                            label reverseIndex = (f.size()-index)%f.size();
+                            scalar d = mag(p[f[reverseIndex]]-nbrPoint[bFaceI]);
+
+                            if (d > smallDist[i])
+                            {
+                                if (setPtr)
+                                {
+                                    // Avoid duplicate counting of faces
+                                    if (setPtr->insert(cpp.start()+i))
+                                    {
+                                        nErrorFaces++;
+                                    }
+                                }
+                                else
+                                {
+                                    // No checking on duplicates
+                                    nErrorFaces++;
+                                }
+                            }
+                            avgMismatch += d;
+                            nCoupledPoints++;
                         }
-                        nErrorFaces++;
                     }
-                    avgMismatch += d;
-                    nCoupledFaces++;
                 }
             }
         }
@@ -347,11 +439,11 @@ bool Foam::checkCoupledPoints
 
     reduce(nErrorFaces, sumOp<label>());
     reduce(avgMismatch, maxOp<scalar>());
-    reduce(nCoupledFaces, sumOp<label>());
+    reduce(nCoupledPoints, sumOp<label>());
 
-    if (nCoupledFaces > 0)
+    if (nCoupledPoints > 0)
     {
-        avgMismatch /= nCoupledFaces;
+        avgMismatch /= nCoupledPoints;
     }
 
     if (nErrorFaces > 0)
@@ -360,7 +452,7 @@ bool Foam::checkCoupledPoints
         {
             Info<< "  **Error in coupled point location: "
                 << nErrorFaces
-                << " faces have their 0th vertex not opposite"
+                << " faces have their 0th or consecutive vertex not opposite"
                 << " their coupled equivalent. Average mismatch "
                 << avgMismatch << "."
                 << endl;
@@ -581,7 +673,8 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
             if (nFaces > 0)
             {
                 Info<< "  <<Writing " << nFaces
-                    << " faces with incorrectly matched 0th vertex to set "
+                    << " faces with incorrectly matched 0th (or consecutive)"
+                    << " vertex to set "
                     << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
diff --git a/applications/utilities/mesh/manipulation/moveDynamicMesh/Make/options b/applications/utilities/mesh/manipulation/moveDynamicMesh/Make/options
index 36fba156e5..f7e149613f 100644
--- a/applications/utilities/mesh/manipulation/moveDynamicMesh/Make/options
+++ b/applications/utilities/mesh/manipulation/moveDynamicMesh/Make/options
@@ -1,9 +1,11 @@
 EXE_INC = \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -ldynamicFvMesh \
     -lmeshTools \
+    -lsampling \
     -ldynamicMesh
diff --git a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C
index 47f1dba097..00c97614de 100644
--- a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C
+++ b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.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
@@ -32,19 +32,98 @@ Description
 #include "argList.H"
 #include "Time.H"
 #include "dynamicFvMesh.H"
+#include "vtkSurfaceWriter.H"
+#include "cyclicAMIPolyPatch.H"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Dump patch + weights to vtk file
+void writeWeights
+(
+    const scalarField& wghtSum,
+    const primitivePatch& patch,
+    const fileName& folder,
+    const fileName& prefix,
+    const word& timeName
+)
+{
+    vtkSurfaceWriter writer;
+
+    writer.write
+    (
+        folder,
+        prefix + "_proc" + Foam::name(Pstream::myProcNo()) + "_" + timeName,
+        patch.localPoints(),
+        patch.localFaces(),
+        "weightsSum",
+        wghtSum,
+        false
+    );
+}
+
+
+void writeWeights(const polyMesh& mesh)
+{
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    const word tmName(mesh.time().timeName());
+
+    forAll(pbm, patchI)
+    {
+        if (isA<cyclicAMIPolyPatch>(pbm[patchI]))
+        {
+            const cyclicAMIPolyPatch& cpp =
+                refCast<const cyclicAMIPolyPatch>(pbm[patchI]);
+
+            if (cpp.owner())
+            {
+                const AMIPatchToPatchInterpolation& ami =
+                    cpp.AMI();
+                writeWeights
+                (
+                    ami.tgtWeightsSum(),
+                    cpp.neighbPatch(),
+                    "output",
+                    "tgt",
+                    tmName
+                );
+                writeWeights
+                (
+                    ami.srcWeightsSum(),
+                    cpp,
+                    "output",
+                    "src",
+                    tmName
+                );
+            }
+        }
+    }
+}
+
+
 // Main program:
 
 int main(int argc, char *argv[])
 {
+    argList::addBoolOption
+    (
+        "checkAMI",
+        "check AMI weights"
+    );
 
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createDynamicFvMesh.H"
 
+    const bool checkAMI  = args.optionFound("checkAMI");
+
+    if (checkAMI)
+    {
+        Info<< "Writing VTK files with weights of AMI patches." << nl << endl;
+    }
+
     while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << endl;
@@ -52,6 +131,11 @@ int main(int argc, char *argv[])
         mesh.update();
         mesh.checkMesh(true);
 
+        if (checkAMI)
+        {
+            writeWeights(mesh);
+        }
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/utilities/postProcessing/wall/wallHeatFlux/Make/options b/applications/utilities/postProcessing/wall/wallHeatFlux/Make/options
index 4d244f56ba..52ab08bdcc 100644
--- a/applications/utilities/postProcessing/wall/wallHeatFlux/Make/options
+++ b/applications/utilities/postProcessing/wall/wallHeatFlux/Make/options
@@ -1,15 +1,18 @@
 EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels \
-    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
-    -lcompressibleRASModels \
+    -lcompressibleTurbulenceModel \
     -lreactionThermophysicalModels \
     -lfiniteVolume \
     -lgenericPatchFields \
     -lspecie \
-    -lbasicThermophysicalModels
+    -lsolid \
+    -lbasicThermophysicalModels \
+    -lbasicSolidThermo
diff --git a/applications/utilities/postProcessing/wall/wallHeatFlux/createFields.H b/applications/utilities/postProcessing/wall/wallHeatFlux/createFields.H
index 4461a3aa55..f1096d657e 100644
--- a/applications/utilities/postProcessing/wall/wallHeatFlux/createFields.H
+++ b/applications/utilities/postProcessing/wall/wallHeatFlux/createFields.H
@@ -5,6 +5,7 @@ autoPtr<basicThermo> thermo
 
 const volScalarField& h = thermo->he();
 
+// Register copy of thermo density
 volScalarField rho
 (
     IOobject
@@ -16,28 +17,40 @@ volScalarField rho
     thermo->rho()
 );
 
-volVectorField U
-(
-    IOobject
+// Construct turbulence model (if fluid)
+autoPtr<volVectorField> UPtr;
+autoPtr<surfaceScalarField> phiPtr;
+autoPtr<compressible::turbulenceModel> turbulence;
+
+if (!isA<solidThermo>(thermo()))
+{
+    UPtr.reset
     (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
+        new volVectorField
+        (
+            IOobject
+            (
+                "U",
+                runTime.timeName(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh
+        )
+    );
+    const volVectorField& U = UPtr();
 
-#include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
+    // Copy phi to autoPtr. Rename to make sure copy is now registered as 'phi'.
+    phi.rename("phiFluid");
+    phiPtr.reset(new surfaceScalarField("phi", phi));
 
-autoPtr<compressible::RASModel> RASModel
-(
-    compressible::RASModel::New
+    turbulence = compressible::turbulenceModel::New
     (
         rho,
         U,
-        phi,
+        phiPtr(),
         thermo()
-    )
-);
+    );
+}
diff --git a/applications/utilities/postProcessing/wall/wallHeatFlux/wallHeatFlux.C b/applications/utilities/postProcessing/wall/wallHeatFlux/wallHeatFlux.C
index 18a0dbf6d4..137c4c0cfb 100644
--- a/applications/utilities/postProcessing/wall/wallHeatFlux/wallHeatFlux.C
+++ b/applications/utilities/postProcessing/wall/wallHeatFlux/wallHeatFlux.C
@@ -32,7 +32,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "RASModel.H"
+#include "turbulenceModel.H"
+#include "solidThermo.H"
 #include "wallFvPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -56,7 +57,14 @@ int main(int argc, char *argv[])
 
         surfaceScalarField heatFlux
         (
-            fvc::interpolate(RASModel->alphaEff())*fvc::snGrad(h)
+            fvc::interpolate
+            (
+                (
+                    turbulence.valid()
+                  ? turbulence->alphaEff()()
+                  : thermo->alpha()
+                )
+            )*fvc::snGrad(h)
         );
 
         const surfaceScalarField::GeometricBoundaryField& patchHeatFlux =
diff --git a/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C b/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C
index 7ed970aff4..d19ea04210 100644
--- a/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C
+++ b/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.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
@@ -93,6 +93,7 @@ int main(int argc, char *argv[])
             U[celli] *= ::pow(y[celli]/yblv, (1.0/7.0));
         }
     }
+    U.correctBoundaryConditions();
 
     Info<< "Writing U\n" << endl;
     U.write();
@@ -109,7 +110,7 @@ int main(int argc, char *argv[])
 
     if (args.optionFound("writenut"))
     {
-        Info<< "Writing nut" << endl;
+        Info<< "Writing " << nut.name() << nl << endl;
         nut.write();
     }
 
@@ -126,7 +127,7 @@ int main(int argc, char *argv[])
     k = sqr(nut/(ck0*min(y, ybl)));
     k.correctBoundaryConditions();
 
-    Info<< "Writing k\n" << endl;
+    Info<< "Writing " << k.name() << nl << endl;
     k.write();
 
 
@@ -137,7 +138,7 @@ int main(int argc, char *argv[])
     epsilon = ce0*k*sqrt(k)/min(y, ybl);
     epsilon.correctBoundaryConditions();
 
-    Info<< "Writing epsilon\n" << endl;
+    Info<< "Writing " << epsilon.name() << nl << endl;
     epsilon.write();
 
 
diff --git a/applications/utilities/preProcessing/applyBoundaryLayer/createFields.H b/applications/utilities/preProcessing/applyBoundaryLayer/createFields.H
index d401808244..d74cadc9bc 100644
--- a/applications/utilities/preProcessing/applyBoundaryLayer/createFields.H
+++ b/applications/utilities/preProcessing/applyBoundaryLayer/createFields.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
@@ -41,9 +41,9 @@ License
 
     singlePhaseTransportModel laminarTransport(U, phi);
 
-    autoPtr<incompressible::RASModel> turbulence
+    autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::RASModel::New(U, phi, laminarTransport)
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
     );
 
     Info<< "Calculating wall distance field" << endl;
diff --git a/bin/tools/CleanFunctions b/bin/tools/CleanFunctions
index be1c3723b4..758627a29a 100644
--- a/bin/tools/CleanFunctions
+++ b/bin/tools/CleanFunctions
@@ -88,7 +88,7 @@ cleanCase()
             rm -rf \
             allOwner* cell* face* meshModifiers* \
             owner* neighbour* point* edge* \
-            cellLevel* pointLevel* refinementHistory* level0Edge surfaceIndex* sets \
+            cellLevel* pointLevel* refinementHistory* level0Edge* surfaceIndex* sets \
             > /dev/null 2>&1 \
         )
     fi
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
index 523cb99b0b..fd9b0a0828 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
@@ -172,7 +172,7 @@ Foam::scalarField Foam::coupledPolyPatch::calcFaceTol
             maxLenSqr = max(maxLenSqr, magSqr(pt - cc));
             maxCmpt = max(maxCmpt, cmptMax(cmptMag(pt)));
         }
-        tols[faceI] = max(SMALL*maxCmpt, Foam::sqrt(maxLenSqr));
+        tols[faceI] = max(SMALL, max(SMALL*maxCmpt, Foam::sqrt(maxLenSqr)));
     }
     return tols;
 }
@@ -200,12 +200,40 @@ Foam::label Foam::coupledPolyPatch::getRotation
         }
     }
 
-    if (anchorFp == -1 || mag(minDistSqr) > tol)
+    if (anchorFp == -1 || Foam::sqrt(minDistSqr) > tol)
     {
         return -1;
     }
     else
     {
+        // Check that anchor is unique.
+        forAll(f, fp)
+        {
+            scalar distSqr = magSqr(anchor - points[f[fp]]);
+
+            if (distSqr == minDistSqr && fp != anchorFp)
+            {
+                WarningIn
+                (
+                    "label coupledPolyPatch::getRotation\n"
+                    "(\n"
+                    "    const pointField&,\n"
+                    "    const face&,\n"
+                    "    const point&,\n"
+                    "    const scalar\n"
+                    ")"
+                )   << "Cannot determine unique anchor point on face "
+                    << UIndirectList<point>(points, f)
+                    << endl
+                    << "Both at index " << anchorFp << " and " << fp
+                    << " the vertices have the same distance "
+                    << Foam::sqrt(minDistSqr)
+                    << " to the anchor " << anchor
+                    << ". Continuing but results might be wrong."
+                    << endl;
+            }
+        }
+
         // Positive rotation
         return (f.size() - anchorFp) % f.size();
     }
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
index 3d35250608..8de9076603 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.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
@@ -118,8 +118,19 @@ void Foam::primitiveMesh::makeFaceCentresAndAreas
                 sumAc += a*c;
             }
 
-            fCtrs[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
-            fAreas[facei] = 0.5*sumN;
+
+            if (sumA < ROOTVSMALL)
+            {
+                // Sum of area too small. No chance of reliably calculating
+                // centroid so fallback to average.
+                fCtrs[facei] = fCentre;
+                fAreas[facei] = 0.5*sumN;
+            }
+            else
+            {
+                fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
+                fAreas[facei] = 0.5*sumN;
+            }
         }
     }
 }
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcVolumeIntegrate.C b/src/finiteVolume/finiteVolume/fvc/fvcVolumeIntegrate.C
index f01c72ce6c..783a20482e 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcVolumeIntegrate.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcVolumeIntegrate.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
@@ -49,6 +49,7 @@ volumeIntegrate
     return vf.mesh().V()*vf.internalField();
 }
 
+
 template<class Type>
 tmp<Field<Type> >
 volumeIntegrate
@@ -62,6 +63,23 @@ volumeIntegrate
 }
 
 
+template<class Type>
+tmp<Field<Type> > volumeIntegrate(const DimensionedField<Type, volMesh>& df)
+{
+    return df.mesh().V()*df.field();
+}
+
+
+template<class Type>
+tmp<Field<Type> >
+volumeIntegrate(const tmp<DimensionedField<Type, volMesh> >& tdf)
+{
+    tmp<Field<Type> > tdidf = tdf().mesh().V()*tdf().field();
+    tdf.clear();
+    return tdidf;
+}
+
+
 template<class Type>
 dimensioned<Type>
 domainIntegrate
@@ -77,9 +95,9 @@ domainIntegrate
     );
 }
 
+
 template<class Type>
-dimensioned<Type>
-domainIntegrate
+dimensioned<Type> domainIntegrate
 (
     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tvf
 )
@@ -90,6 +108,33 @@ domainIntegrate
 }
 
 
+template<class Type>
+dimensioned<Type> domainIntegrate
+(
+    const DimensionedField<Type, volMesh>& df
+)
+{
+    return dimensioned<Type>
+    (
+        "domainIntegrate(" + df.name() + ')',
+        dimVol*df.dimensions(),
+        gSum(fvc::volumeIntegrate(df))
+    );
+}
+
+
+template<class Type>
+dimensioned<Type> domainIntegrate
+(
+    const tmp<DimensionedField<Type, volMesh> >& tdf
+)
+{
+    dimensioned<Type> integral = domainIntegrate(tdf());
+    tdf.clear();
+    return integral;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fvc
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcVolumeIntegrate.H b/src/finiteVolume/finiteVolume/fvc/fvcVolumeIntegrate.H
index 741d1451b4..784087ee30 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcVolumeIntegrate.H
+++ b/src/finiteVolume/finiteVolume/fvc/fvcVolumeIntegrate.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
@@ -67,6 +67,19 @@ namespace fvc
     );
 
 
+    template<class Type>
+    tmp<Field<Type> > volumeIntegrate
+    (
+        const DimensionedField<Type, volMesh>&
+    );
+
+    template<class Type>
+    tmp<Field<Type> > volumeIntegrate
+    (
+        const tmp<DimensionedField<Type, volMesh> >&
+    );
+
+
     template<class Type>
     dimensioned<Type> domainIntegrate
     (
@@ -78,6 +91,19 @@ namespace fvc
     (
         const tmp<GeometricField<Type, fvPatchField, volMesh> >&
     );
+
+
+    template<class Type>
+    dimensioned<Type> domainIntegrate
+    (
+        const DimensionedField<Type, volMesh>&
+    );
+
+    template<class Type>
+    dimensioned<Type> domainIntegrate
+    (
+        const tmp<DimensionedField<Type, volMesh> >&
+    );
 }
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index 1bc817a4c1..eb3a3fc11c 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -336,7 +336,6 @@ void Foam::ReactingCloud<CloudType>::evolve()
 }
 
 
-
 template<class CloudType>
 void Foam::ReactingCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
 {
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 8fa42ac0ec..58217adbf5 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -169,14 +169,15 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
     const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
     const forceSuSp Feff = Fcp + Fncp;
+    const scalar massEff = forces.massEff(p, mass);
 
 
     // New particle velocity
     //~~~~~~~~~~~~~~~~~~~~~~
 
     // Update velocity - treat as 3-D
-    const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/mass;
-    const scalar bp = Feff.Sp()/mass;
+    const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/massEff;
+    const scalar bp = Feff.Sp()/massEff;
 
     Spu = dt*Feff.Sp();
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index ba48774d69..4e4e0bf7be 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -67,19 +67,31 @@ inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
 )
 :
     dict_(parentDict.subOrEmptyDict("constantProperties")),
-    parcelTypeId_(-1),
-    rhoMin_(0.0),
+    parcelTypeId_(1),
+    rhoMin_(1e-15),
     rho0_(0.0),
-    minParticleMass_(0.0),
+    minParticleMass_(1e-15),
     youngsModulus_(0.0),
     poissonsRatio_(0.0)
 {
     if (readFields)
     {
-        dict_.lookup("parcelTypeId") >> parcelTypeId_;
-        dict_.lookup("rhoMin") >> rhoMin_;
+        if (dict_.readIfPresent("parcelTypeId", parcelTypeId_))
+        {
+            Info<< "    employing parcel parcelTypeId of " << parcelTypeId_
+                << endl;
+        }
+        if (dict_.readIfPresent("rhoMin", rhoMin_))
+        {
+            Info<< "    employing parcel rhoMin of " << rhoMin_ << endl;
+        }
+        if (dict_.readIfPresent("minParticleMass", minParticleMass_))
+        {
+            Info<< "    employing parcel minParticleMass of "
+                << minParticleMass_ << endl;
+        }
+
         dict_.lookup("rho0") >> rho0_;
-        dict_.lookup("minParticleMass") >> minParticleMass_;
         dict_.lookup("youngsModulus") >> youngsModulus_;
         dict_.lookup("poissonsRatio") >> poissonsRatio_;
     }
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index 230a4c1cd3..0bb8f3704e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -115,6 +115,7 @@ public:
                 const scalar poissonsRatio,
                 const scalar T0,
                 const scalar TMin,
+                const scalar TMax,
                 const scalar Cp0,
                 const scalar epsilon0,
                 const scalar f0,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
index 04b85b99dc..e6c140973f 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
@@ -59,14 +59,18 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
 )
 :
     ParcelType::constantProperties(parentDict, readFields),
-    pMin_(0.0),
+    pMin_(1000.0),
     constantVolume_(false),
     Tvap_(0.0),
     Tbp_(0.0)
 {
     if (readFields)
     {
-        this->dict().lookup("pMin") >> pMin_;
+        if (this->dict().readIfPresent("pMin", pMin_))
+        {
+            Info<< "    employing parcel pMin of " << pMin_ << endl;
+        }
+
         this->dict().lookup("constantVolume") >> constantVolume_;
         this->dict().lookup("Tvap") >> Tvap_;
         this->dict().lookup("Tbp") >> Tbp_;
@@ -85,6 +89,7 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
     const scalar poissonsRatio,
     const scalar T0,
     const scalar TMin,
+    const scalar TMax,
     const scalar Cp0,
     const scalar epsilon0,
     const scalar f0,
@@ -105,6 +110,7 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
         poissonsRatio,
         T0,
         TMin,
+        TMax,
         Cp0,
         epsilon0,
         f0,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index cb9b714204..7baa54c1eb 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.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
@@ -324,7 +324,16 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     IntegrationScheme<scalar>::integrationResult Tres =
         td.cloud().TIntegrator().integrate(T_, dt, ap*bp, bp);
 
-    scalar Tnew = max(Tres.value(), td.cloud().constProps().TMin());
+    scalar Tnew =
+        min
+        (
+            max
+            (
+                Tres.value(),
+                td.cloud().constProps().TMin()
+            ),
+            td.cloud().constProps().TMax()
+        );
 
     Sph = dt*htc*As;
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index e7b5a14a49..232b45d471 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -82,6 +82,9 @@ public:
             //- Minimum temperature [K]
             scalar TMin_;
 
+            //- Maximum temperature [K]
+            scalar TMax_;
+
             //- Particle specific heat capacity [J/(kg.K)]
             scalar Cp0_;
 
@@ -123,6 +126,7 @@ public:
                 const scalar poissonsRatio,
                 const scalar T0,
                 const scalar TMin,
+                const scalar TMax,
                 const scalar Cp0,
                 const scalar epsilon0,
                 const scalar f0,
@@ -140,6 +144,9 @@ public:
                 //- Return const access to minimum temperature [K]
                 inline scalar TMin() const;
 
+                //- Return const access to maximum temperature [K]
+                inline scalar TMax() const;
+
                 //- Return const access to the particle specific heat capacity
                 //  [J/(kg.K)]
                 inline scalar Cp0() const;
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
index e5c7e349c5..85f7a0c3d5 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
@@ -31,6 +31,7 @@ inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties()
     ParcelType::constantProperties(),
     T0_(0.0),
     TMin_(0.0),
+    TMax_(VGREAT),
     Cp0_(0.0),
     epsilon0_(0.0),
     f0_(0.0),
@@ -47,6 +48,7 @@ inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
     ParcelType::constantProperties(cp),
     T0_(cp.T0_),
     TMin_(cp.TMin_),
+    TMax_(cp.TMax_),
     Cp0_(cp.Cp0_),
     epsilon0_(cp.epsilon0_),
     f0_(cp.f0_),
@@ -63,7 +65,8 @@ inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
 :
     ParcelType::constantProperties(parentDict, readFields),
     T0_(0.0),
-    TMin_(0.0),
+    TMin_(200),
+    TMax_(5000),
     Cp0_(0.0),
     epsilon0_(0.0),
     f0_(0.0),
@@ -71,8 +74,16 @@ inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
 {
     if (readFields)
     {
+        if (this->dict().readIfPresent("TMin", TMin_))
+        {
+            Info<< "    employing parcel TMin of " << TMin_ << endl;
+        }
+        if (this->dict().readIfPresent("TMax", TMax_))
+        {
+            Info<< "    employing parcel TMax of " << TMax_ << endl;
+        }
+
         this->dict().lookup("T0") >> T0_;
-        this->dict().lookup("TMin") >> TMin_;
         this->dict().lookup("Cp0") >> Cp0_;
         this->dict().lookup("epsilon0") >> epsilon0_;
         this->dict().lookup("f0") >> f0_;
@@ -92,6 +103,7 @@ inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
     const scalar poissonsRatio,
     const scalar T0,
     const scalar TMin,
+    const scalar TMax,
     const scalar Cp0,
     const scalar epsilon0,
     const scalar f0,
@@ -109,6 +121,7 @@ inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
     ),
     T0_(T0),
     TMin_(TMin),
+    TMax_(TMax),
     Cp0_(Cp0),
     epsilon0_(epsilon0),
     f0_(f0),
@@ -195,6 +208,14 @@ Foam::ThermoParcel<ParcelType>::constantProperties::TMin() const
 }
 
 
+template<class ParcelType>
+inline Foam::scalar
+Foam::ThermoParcel<ParcelType>::constantProperties::TMax() const
+{
+    return TMax_;
+}
+
+
 template<class ParcelType>
 inline Foam::scalar
 Foam::ThermoParcel<ParcelType>::constantProperties::Cp0() const
diff --git a/src/lagrangian/intermediate/parcels/include/makeParcelForces.H b/src/lagrangian/intermediate/parcels/include/makeParcelForces.H
index 81e1658000..7658900598 100644
--- a/src/lagrangian/intermediate/parcels/include/makeParcelForces.H
+++ b/src/lagrangian/intermediate/parcels/include/makeParcelForces.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
@@ -36,6 +36,7 @@ License
 #include "ParamagneticForce.H"
 #include "PressureGradientForce.H"
 #include "SRFForce.H"
+#include "VirtualMassForce.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,7 +49,8 @@ License
     makeParticleForceModelType(NonInertialFrameForce, CloudType);             \
     makeParticleForceModelType(ParamagneticForce, CloudType);                 \
     makeParticleForceModelType(PressureGradientForce, CloudType);             \
-    makeParticleForceModelType(SRFForce, CloudType);
+    makeParticleForceModelType(SRFForce, CloudType);                          \
+    makeParticleForceModelType(VirtualMassForce, CloudType);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/intermediate/parcels/include/makeThermoParcelForces.H b/src/lagrangian/intermediate/parcels/include/makeThermoParcelForces.H
index f71ad43d5e..dff1f464a1 100644
--- a/src/lagrangian/intermediate/parcels/include/makeThermoParcelForces.H
+++ b/src/lagrangian/intermediate/parcels/include/makeThermoParcelForces.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
@@ -37,6 +37,7 @@ License
 #include "ParamagneticForce.H"
 #include "PressureGradientForce.H"
 #include "SRFForce.H"
+#include "VirtualMassForce.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,7 +51,8 @@ License
     makeParticleForceModelType(NonInertialFrameForce, CloudType);             \
     makeParticleForceModelType(ParamagneticForce, CloudType);                 \
     makeParticleForceModelType(PressureGradientForce, CloudType);             \
-    makeParticleForceModelType(SRFForce, CloudType);
+    makeParticleForceModelType(SRFForce, CloudType);                          \
+    makeParticleForceModelType(VirtualMassForce, CloudType);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C
index c4ad6aa5d7..04016823b4 100644
--- a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C
+++ b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.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
@@ -180,4 +180,21 @@ Foam::forceSuSp Foam::ParticleForceList<CloudType>::calcNonCoupled
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::ParticleForceList<CloudType>::massEff
+(
+    const typename CloudType::parcelType& p,
+    const scalar mass
+) const
+{
+    scalar massEff = mass;
+    forAll(*this, i)
+    {
+        massEff += this->operator[](i).massAdd(p, mass);
+    }
+
+    return massEff;
+}
+
+
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H
index 5eda11a0f5..fe5f561268 100644
--- a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H
+++ b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H
@@ -130,6 +130,13 @@ public:
                 const scalar Re,
                 const scalar muc
             ) const;
+
+            //- Return the effective mass
+            virtual scalar massEff
+            (
+                const typename CloudType::parcelType& p,
+                const scalar mass
+            ) const;
 };
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.C
index fc4ba1314b..7adbb47ebd 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.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
@@ -27,7 +27,6 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-
 template<class CloudType>
 Foam::ParticleForce<CloudType>::ParticleForce
 (
@@ -120,6 +119,17 @@ Foam::forceSuSp Foam::ParticleForce<CloudType>::calcNonCoupled
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::ParticleForce<CloudType>::massAdd
+(
+    const typename CloudType::parcelType& p,
+    const scalar mass
+) const
+{
+    return 0.0;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "ParticleForceNew.C"
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.H
index c5a6cf01d7..a72abc3217 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.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
@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                      Class ParticleForce Declaration
+                        Class ParticleForce Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class CloudType>
@@ -169,6 +169,13 @@ public:
                 const scalar Re,
                 const scalar muc
             ) const;
+
+            //- Return the added mass
+            virtual scalar massAdd
+            (
+                const typename CloudType::parcelType& p,
+                const scalar mass
+            ) const;
 };
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C
new file mode 100644
index 0000000000..3dca602dd4
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C
@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "VirtualMassForce.H"
+#include "fvcDdt.H"
+#include "fvcGrad.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::VirtualMassForce<CloudType>::VirtualMassForce
+(
+    CloudType& owner,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
+    UName_(this->coeffs().template lookupOrDefault<word>("U", "U")),
+    Cvm_(readScalar(this->coeffs().lookup("Cvm"))),
+    DUcDtPtr_(NULL),
+    DUcDtInterpPtr_(NULL)
+{}
+
+
+template<class CloudType>
+Foam::VirtualMassForce<CloudType>::VirtualMassForce
+(
+    const VirtualMassForce& vmf
+)
+:
+    ParticleForce<CloudType>(vmf),
+    UName_(vmf.UName_),
+    Cvm_(vmf.Cvm_),
+    DUcDtPtr_(NULL),
+    DUcDtInterpPtr_(NULL)
+{}
+
+
+// * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::VirtualMassForce<CloudType>::~VirtualMassForce()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::VirtualMassForce<CloudType>::cacheFields(const bool store)
+{
+    if (store && !DUcDtPtr_)
+    {
+        const volVectorField& Uc = this->mesh().template
+            lookupObject<volVectorField>(UName_);
+
+        DUcDtPtr_ = new volVectorField
+        (
+            "DUcDt",
+            fvc::ddt(Uc) + (Uc & fvc::grad(Uc))
+        );
+
+        DUcDtInterpPtr_.reset
+        (
+            interpolation<vector>::New
+            (
+                this->owner().solution().interpolationSchemes(),
+                *DUcDtPtr_
+            ).ptr()
+        );
+    }
+    else
+    {
+        DUcDtInterpPtr_.clear();
+
+        if (DUcDtPtr_)
+        {
+            delete DUcDtPtr_;
+            DUcDtPtr_ = NULL;
+        }
+    }
+}
+
+
+template<class CloudType>
+Foam::forceSuSp Foam::VirtualMassForce<CloudType>::calcCoupled
+(
+    const typename CloudType::parcelType& p,
+    const scalar dt,
+    const scalar mass,
+    const scalar Re,
+    const scalar muc
+) const
+{
+    forceSuSp value(vector::zero, 0.0);
+
+    vector DUcDt =
+        DUcDtInterp().interpolate(p.position(), p.currentTetIndices());
+
+    value.Su() = mass*p.rhoc()/p.rho()*Cvm_*DUcDt;
+
+    return value;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::VirtualMassForce<CloudType>::massAdd
+(
+    const typename CloudType::parcelType& p,
+    const scalar mass
+) const
+{
+    return mass*p.rhoc()/p.rho()*Cvm_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H
new file mode 100644
index 0000000000..c5277ec696
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H
@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::VirtualMassForce
+
+Description
+    Calculates particle virtual mass force
+
+SourceFiles
+    VirtualMassForceI.H
+    VirtualMassForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef VirtualMassForce_H
+#define VirtualMassForce_H
+
+#include "ParticleForce.H"
+#include "volFields.H"
+#include "interpolation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class VirtualMassForce Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class VirtualMassForce
+:
+    public ParticleForce<CloudType>
+{
+    // Private data
+
+        //- Name of velocity field
+        const word UName_;
+
+        //- Virtual mass coefficient - typically 0.5
+        scalar Cvm_;
+
+        //- Rate of change of carrier phase velocity
+        volVectorField* DUcDtPtr_;
+
+        //- Rate of change of carrier phase velocity interpolator
+        autoPtr<interpolation<vector> > DUcDtInterpPtr_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("virtualMass");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        VirtualMassForce
+        (
+            CloudType& owner,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct copy
+        VirtualMassForce(const VirtualMassForce& pgf);
+
+        //- Construct and return a clone
+        virtual autoPtr<ParticleForce<CloudType> > clone() const
+        {
+            return autoPtr<ParticleForce<CloudType> >
+            (
+                new VirtualMassForce<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~VirtualMassForce();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the rate of change of carrier phase velocity
+            inline const volVectorField& DUcDt() const;
+
+            //- Return the rate of change of carrier phase velocity interpolator
+            inline const interpolation<vector>& DUcDtInterp() const;
+
+
+        // Evaluation
+
+            //- Cache fields
+            virtual void cacheFields(const bool store);
+
+            //- Calculate the non-coupled force
+            virtual forceSuSp calcCoupled
+            (
+                const typename CloudType::parcelType& p,
+                const scalar dt,
+                const scalar mass,
+                const scalar Re,
+                const scalar muc
+            ) const;
+
+            //- Return the added mass
+            virtual scalar massAdd
+            (
+                const typename CloudType::parcelType& p,
+                const scalar mass
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "VirtualMassForceI.H"
+
+#ifdef NoRepository
+    #include "VirtualMassForce.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForceI.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForceI.H
new file mode 100644
index 0000000000..e156c58d14
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForceI.H
@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+const Foam::volVectorField& Foam::VirtualMassForce<CloudType>::DUcDt() const
+{
+    if (DUcDtPtr_)
+    {
+        return *DUcDtPtr_;
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "const volVectorField& VirtualMassForce<CloudType>::DUcDt()"
+            "const"
+        )   << "DUcDt field not allocated" << abort(FatalError);
+
+        return *reinterpret_cast<const volVectorField*>(0);
+    }
+}
+
+
+template<class CloudType>
+inline const Foam::interpolation<Foam::vector>&
+Foam::VirtualMassForce<CloudType>::DUcDtInterp() const
+{
+    if (!DUcDtInterpPtr_.valid())
+    {
+        FatalErrorIn
+        (
+            "inline const Foam::interpolation<Foam::vector>&"
+            "Foam::VirtualMassForce<CloudType>::DUcDtInterp() const"
+        )   << "Carrier pahase DUcDt interpolation object not set"
+            << abort(FatalError);
+    }
+
+    return DUcDtInterpPtr_();
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
index d1f66fb6e9..95d9bf0dd4 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.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
@@ -28,209 +28,14 @@ License
 #include "unitConversion.H"
 #include "refinementSurfaces.H"
 #include "searchableSurfaces.H"
-#include "regExp.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 const Foam::scalar Foam::layerParameters::defaultConcaveAngle = 90;
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-// Read the number of layers from dictionary. Per patch 0 or the number
-// of layers.
-Foam::labelList Foam::layerParameters::readNumLayers
-(
-    const PtrList<dictionary>& surfaceDicts,
-    const refinementSurfaces& refineSurfaces,
-    const labelList& globalToPatch,
-    const polyBoundaryMesh& boundaryMesh
-)
-{
-    // Per surface the number of layers
-    labelList globalSurfLayers(surfaceDicts.size());
-    // Per surface, per region the number of layers
-    List<Map<label> > regionSurfLayers(surfaceDicts.size());
-
-    const labelList& surfaceIndices = refineSurfaces.surfaces();
-
-    forAll(surfaceDicts, surfI)
-    {
-        const dictionary& dict = surfaceDicts[surfI];
-
-        globalSurfLayers[surfI] = readLabel(dict.lookup("surfaceLayers"));
-
-        if (dict.found("regions"))
-        {
-            // Per-region layer information
-
-            PtrList<dictionary> regionDicts(dict.lookup("regions"));
-
-            const wordList& regionNames =
-                refineSurfaces.geometry()[surfaceIndices[surfI]].regions();
-
-            forAll(regionDicts, dictI)
-            {
-                const dictionary& regionDict = regionDicts[dictI];
-
-                const word regionName(regionDict.lookup("name"));
-
-                label regionI = findIndex(regionNames, regionName);
-
-                label nLayers = readLabel(regionDict.lookup("surfaceLayers"));
-
-                Info<< "    region " << regionName << ':'<< nl
-                    << "        surface layers:" << nLayers << nl;
-
-                regionSurfLayers[surfI].insert(regionI, nLayers);
-            }
-        }
-    }
-
-
-    // Transfer per surface/region information into patchwise region info
-
-    labelList nLayers(boundaryMesh.size(), 0);
-
-    forAll(surfaceIndices, surfI)
-    {
-        const wordList& regionNames =
-            refineSurfaces.geometry()[surfaceIndices[surfI]].regions();
-
-        forAll(regionNames, regionI)
-        {
-            const word& regionName = regionNames[regionI];
-
-            label global = refineSurfaces.globalRegion(surfI, regionI);
-
-            label patchI = globalToPatch[global];
-
-            // Initialise to surface-wise layers
-            nLayers[patchI] = globalSurfLayers[surfI];
-
-            // Override with region specific data if available
-            Map<label>::const_iterator iter =
-                regionSurfLayers[surfI].find(regionI);
-
-            if (iter != regionSurfLayers[surfI].end())
-            {
-                nLayers[patchI] = iter();
-            }
-
-            // Check
-            if (nLayers[patchI] < 0)
-            {
-                FatalErrorIn
-                (
-                    "layerParameters::readNumLayers(..)"
-                )   << "Illegal number of layers " << nLayers[patchI]
-                    << " for surface "
-                    << refineSurfaces.names()[surfI]
-                    << " region " << regionName << endl
-                    << exit(FatalError);
-            }
-        }
-    }
-    return nLayers;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-//// Construct from dictionary
-//Foam::layerParameters::layerParameters
-//(
-//    const PtrList<dictionary>& surfaceDicts,
-//    const refinementSurfaces& refineSurfaces,
-//    const labelList& globalToPatch,
-//    const dictionary& dict,
-//    const polyBoundaryMesh& boundaryMesh
-//)
-//:
-//    numLayers_
-//    (
-//        readNumLayers
-//        (
-//            surfaceDicts,
-//            refineSurfaces,
-//            globalToPatch,
-//            boundaryMesh
-//        )
-//    ),
-//    expansionRatio_
-//    (
-//        numLayers_.size(),
-//        readScalar(dict.lookup("expansionRatio"))
-//    ),
-//    relativeSizes_(false),
-//    finalLayerThickness_
-//    (
-//        numLayers_.size(),
-//        readScalar(dict.lookup("finalLayerRatio"))
-//    ),
-//    minThickness_
-//    (
-//        numLayers_.size(),
-//        readScalar(dict.lookup("minThickness"))
-//    ),
-//    featureAngle_(readScalar(dict.lookup("featureAngle"))),
-//    concaveAngle_
-//    (
-//        dict.lookupOrDefault("concaveAngle", defaultConcaveAngle)
-//    ),
-//    nGrow_(readLabel(dict.lookup("nGrow"))),
-//    nSmoothSurfaceNormals_
-//    (
-//        readLabel(dict.lookup("nSmoothSurfaceNormals"))
-//    ),
-//    nSmoothNormals_(readLabel(dict.lookup("nSmoothNormals"))),
-//    nSmoothThickness_(readLabel(dict.lookup("nSmoothThickness"))),
-//    maxFaceThicknessRatio_
-//    (
-//        readScalar(dict.lookup("maxFaceThicknessRatio"))
-//    ),
-//    layerTerminationCos_
-//    (
-//        Foam::cos(degToRad(0.5*featureAngle_))
-//    ),
-//    maxThicknessToMedialRatio_
-//    (
-//        readScalar(dict.lookup("maxThicknessToMedialRatio"))
-//    ),
-//    minMedianAxisAngleCos_
-//    (
-//        Foam::cos(degToRad(readScalar(dict.lookup("minMedianAxisAngle"))))
-//    ),
-//    nBufferCellsNoExtrude_
-//    (
-//        readLabel(dict.lookup("nBufferCellsNoExtrude"))
-//    ),
-//    nSnap_(readLabel(dict.lookup("nSnap"))),
-//    nLayerIter_(readLabel(dict.lookup("nLayerIter"))),
-//    nRelaxedIter_(labelMax)
-//{
-//    if (nGrow_ > 0)
-//    {
-//        WarningIn("layerParameters::layerParameters(..)")
-//            << "The nGrow parameter effect has changed with respect to 1.6.x."
-//            << endl
-//            << "Please set nGrow=0 for 1.6.x behaviour."
-//            << endl;
-//    }
-//
-//    dict.readIfPresent("nRelaxedIter", nRelaxedIter_);
-//
-//    if (nLayerIter_ < 0 || nRelaxedIter_ < 0)
-//    {
-//        FatalErrorIn("layerParameters::layerParameters(..)")
-//            << "Layer iterations should be >= 0." << endl
-//            << "nLayerIter:" << nLayerIter_
-//            << " nRelaxedIter:" << nRelaxedIter_
-//            << exit(FatalError);
-//    }
-//}
-
-
 // Construct from dictionary
 Foam::layerParameters::layerParameters
 (
@@ -315,71 +120,51 @@ Foam::layerParameters::layerParameters
 
     const dictionary& layersDict = dict.subDict("layers");
 
-    forAll(boundaryMesh, patchI)
+    forAllConstIter(dictionary, layersDict, iter)
     {
-        const word& patchName = boundaryMesh[patchI].name();
-
-        if (layersDict.found(patchName))
+        if (iter().isDict())
         {
-            const dictionary& layerDict = layersDict.subDict(patchName);
-
-            numLayers_[patchI] =
-                readLabel(layerDict.lookup("nSurfaceLayers"));
-
-            layerDict.readIfPresent
+            const word& key = iter().keyword();
+            const labelHashSet patchIDs
             (
-                "expansionRatio",
-                expansionRatio_[patchI]
+                boundaryMesh.patchSet(List<wordRe>(1, key))
             );
-            layerDict.readIfPresent
-            (
-                "finalLayerThickness",
-                finalLayerThickness_[patchI]
-            );
-            layerDict.readIfPresent
-            (
-                "minThickness",
-                minThickness_[patchI]
-            );
-        }
-    }
 
-
-    // Check whether layer specification matches any patches
-    const List<keyType> wildCards = layersDict.keys(true);
-
-    forAll(wildCards, i)
-    {
-        regExp re(wildCards[i]);
-
-        bool hasMatch = false;
-        forAll(boundaryMesh, patchI)
-        {
-            if (re.match(boundaryMesh[patchI].name()))
+            if (patchIDs.size() == 0)
             {
-                hasMatch = true;
-                break;
+                IOWarningIn("layerParameters::layerParameters(..)", layersDict)
+                    << "Layer specification for " << key
+                    << " does not match any patch." << endl
+                    << "Valid patches are " << boundaryMesh.names() << endl;
             }
-        }
-        if (!hasMatch)
-        {
-            IOWarningIn("layerParameters::layerParameters(..)", layersDict)
-                << "Wildcard layer specification for " << wildCards[i]
-                << " does not match any patch." << endl
-                << "Valid patches are " << boundaryMesh.names() << endl;
-        }
-    }
+            else
+            {
+                const dictionary& layerDict = iter().dict();
 
-    const List<keyType> nonWildCards = layersDict.keys(false);
+                forAllConstIter(labelHashSet, patchIDs, patchIter)
+                {
+                    label patchI = patchIter.key();
 
-    forAll(nonWildCards, i)
-    {
-        if (boundaryMesh.findPatchID(nonWildCards[i]) == -1)
-        {
-            IOWarningIn("layerParameters::layerParameters(..)", layersDict)
-                << "Layer specification for " << nonWildCards[i]
-                << " does not match any patch." << endl
-                << "Valid patches are " << boundaryMesh.names() << endl;
+                    numLayers_[patchI] =
+                        readLabel(layerDict.lookup("nSurfaceLayers"));
+
+                    layerDict.readIfPresent
+                    (
+                        "expansionRatio",
+                        expansionRatio_[patchI]
+                    );
+                    layerDict.readIfPresent
+                    (
+                        "finalLayerThickness",
+                        finalLayerThickness_[patchI]
+                    );
+                    layerDict.readIfPresent
+                    (
+                        "minThickness",
+                        minThickness_[patchI]
+                    );
+                }
+            }
         }
     }
 }
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
index e6324131e9..1a7913a8ef 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.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,15 +109,6 @@ class layerParameters
 
     // Private Member Functions
 
-        //- Extract patch-wise number of layers
-        static labelList readNumLayers
-        (
-            const PtrList<dictionary>& surfaceDicts,
-            const refinementSurfaces& refineSurfaces,
-            const labelList& globalToPatch,
-            const polyBoundaryMesh& boundaryMesh
-        );
-
         //- Disallow default bitwise copy construct
         layerParameters(const layerParameters&);
 
@@ -129,17 +120,7 @@ public:
 
     // Constructors
 
-        ////- Construct from dictionary - old syntax
-        //layerParameters
-        //(
-        //    const PtrList<dictionary>& surfaceDicts,
-        //    const refinementSurfaces& refineSurfaces,
-        //    const labelList& globalToPatch,
-        //    const dictionary& dict,
-        //    const polyBoundaryMesh& boundaryMesh
-        //);
-
-        //- Construct from dictionary - new syntax
+        //- Construct from dictionary
         layerParameters(const dictionary& dict, const polyBoundaryMesh&);
 
 
diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
index 56a669b119..8736987ca9 100644
--- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
+++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
@@ -129,17 +129,9 @@ extern "C"
 
 // Hack: scotch generates floating point errors so need to switch of error
 //       trapping!
-#if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
-#    define LINUX
-#endif
-
-#if defined(LINUX) && defined(__GNUC__)
-#    define LINUX_GNUC
-#endif
-
-#ifdef LINUX_GNUC
-#   ifndef __USE_GNU
-#       define __USE_GNU
+#ifdef __GLIBC__
+#   ifndef _GNU_SOURCE
+#       define _GNU_SOURCE
 #   endif
 #   include <fenv.h>
 #endif
@@ -651,7 +643,7 @@ Foam::label Foam::ptscotchDecomp::decompose
 
 
     // Hack:switch off fpu error trapping
-#   ifdef LINUX_GNUC
+#   ifdef  FE_NOMASK_ENV
     int oldExcepts = fedisableexcept
     (
         FE_DIVBYZERO
@@ -681,7 +673,7 @@ Foam::label Foam::ptscotchDecomp::decompose
         "SCOTCH_graphMap"
     );
 
-#   ifdef LINUX_GNUC
+#   ifdef  FE_NOMASK_ENV
     feenableexcept(oldExcepts);
 #   endif
 
diff --git a/src/parallel/decompose/scotchDecomp/scotchDecomp.C b/src/parallel/decompose/scotchDecomp/scotchDecomp.C
index d81af37b00..e8bfe1aa10 100644
--- a/src/parallel/decompose/scotchDecomp/scotchDecomp.C
+++ b/src/parallel/decompose/scotchDecomp/scotchDecomp.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
@@ -135,17 +135,9 @@ extern "C"
 
 // Hack: scotch generates floating point errors so need to switch of error
 //       trapping!
-#if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
-#    define LINUX
-#endif
-
-#if defined(LINUX) && defined(__GNUC__)
-#    define LINUX_GNUC
-#endif
-
-#ifdef LINUX_GNUC
-#   ifndef __USE_GNU
-#       define __USE_GNU
+#ifdef __GLIBC__
+#   ifndef _GNU_SOURCE
+#       define _GNU_SOURCE
 #   endif
 #   include <fenv.h>
 #endif
@@ -503,7 +495,7 @@ Foam::label Foam::scotchDecomp::decomposeOneProc
 
 
     // Hack:switch off fpu error trapping
-#   ifdef LINUX_GNUC
+#   ifdef FE_NOMASK_ENV
     int oldExcepts = fedisableexcept
     (
         FE_DIVBYZERO
@@ -526,7 +518,7 @@ Foam::label Foam::scotchDecomp::decomposeOneProc
         "SCOTCH_graphMap"
     );
 
-#   ifdef LINUX_GNUC
+#   ifdef FE_NOMASK_ENV
     feenableexcept(oldExcepts);
 #   endif
 
diff --git a/src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModel.H b/src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModel.H
index 7e4d0f4cc6..57d8caa3c4 100644
--- a/src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModel.H
+++ b/src/regionModels/pyrolysisModels/pyrolysisModel/pyrolysisModel.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,7 +161,7 @@ public:
 
     // Selectors
 
-        //- Return a reference to the selected pyrolysis film model
+        //- Return a reference to the selected pyrolysis model
         static autoPtr<pyrolysisModel> New(const fvMesh& mesh);
 
         //- Return a reference to a named selected pyrolysis model
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
index d20552f368..05c6c7b9a2 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
@@ -153,10 +153,10 @@ void reactingOneDim::updatePhiGas()
         tmp<volScalarField> tHsiGas =
             solidChemistry_->gasHs(solidThermo_.p(), solidThermo_.T(), gasI);
 
-        tmp<volScalarField> tRRiGas = solidChemistry_->RRg(gasI);
-
         const volScalarField& HsiGas = tHsiGas();
-        const volScalarField& RRiGas = tRRiGas();
+
+        const DimensionedField<scalar, volMesh>& RRiGas =
+            solidChemistry_->RRg(gasI);
 
         label totalFaceId = 0;
         forAll(intCoupledPatchIDs_, i)
diff --git a/src/regionModels/regionModel/regionModel/regionModel.H b/src/regionModels/regionModel/regionModel/regionModel.H
index b65b3a23c7..6fbad3d372 100644
--- a/src/regionModels/regionModel/regionModel/regionModel.H
+++ b/src/regionModels/regionModel/regionModel/regionModel.H
@@ -262,7 +262,7 @@ public:
             //- Evolve the region
             virtual void evolveRegion();
 
-            //- Evolve the film
+            //- Evolve the region
             virtual void evolve();
 
 
diff --git a/src/thermophysicalModels/basicSolidThermo/solidThermo/makeSolidThermo.H b/src/thermophysicalModels/basicSolidThermo/solidThermo/makeSolidThermo.H
index c6fed868b5..4a911e8c0f 100644
--- a/src/thermophysicalModels/basicSolidThermo/solidThermo/makeSolidThermo.H
+++ b/src/thermophysicalModels/basicSolidThermo/solidThermo/makeSolidThermo.H
@@ -33,6 +33,7 @@ Description
 #define makesolidThermo_H
 
 #include "addToRunTimeSelectionTable.H"
+#include "basicThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
@@ -127,6 +128,12 @@ addToRunTimeSelectionTable                                                    \
     BaseThermo,                                                               \
     Cthermo##Mixture##Transport##Radiation##Type##Thermo##Rho##BaseThermo,    \
     mesh                                                                      \
+);                                                                            \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    basicThermo,                                                              \
+    Cthermo##Mixture##Transport##Radiation##Type##Thermo##Rho##BaseThermo,    \
+    fvMesh                                                                    \
 );                                                                            \
                                                                               \
 addToRunTimeSelectionTable                                                    \
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C
index c18cccb712..522afe5161 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C
@@ -64,7 +64,19 @@ Foam::ODEChemistryModel<CompType, ThermoType>::ODEChemistryModel
         RR_.set
         (
             fieldI,
-            new scalarField(mesh.nCells(), 0.0)
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    "RR::" + Y_[fieldI].name(),
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+            )
         );
     }
 
@@ -687,6 +699,11 @@ Foam::label Foam::ODEChemistryModel<CompType, ThermoType>::nEqns() const
 template<class CompType, class ThermoType>
 void Foam::ODEChemistryModel<CompType, ThermoType>::calculate()
 {
+    if (!this->chemistry_)
+    {
+        return;
+    }
+
     const volScalarField rho
     (
         IOobject
@@ -701,36 +718,24 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::calculate()
         this->thermo().rho()
     );
 
-    if (this->mesh().changing())
+    forAll(rho, celli)
     {
+        const scalar rhoi = rho[celli];
+        const scalar Ti = this->thermo().T()[celli];
+        const scalar pi = this->thermo().p()[celli];
+
+        scalarField c(nSpecie_, 0.0);
         for (label i=0; i<nSpecie_; i++)
         {
-            RR_[i].setSize(rho.size());
-            RR_[i] = 0.0;
+            const scalar Yi = Y_[i][celli];
+            c[i] = rhoi*Yi/specieThermo_[i].W();
         }
-    }
 
-    if (this->chemistry_)
-    {
-        forAll(rho, celli)
+        const scalarField dcdt(omega(c, Ti, pi));
+
+        for (label i=0; i<nSpecie_; i++)
         {
-            const scalar rhoi = rho[celli];
-            const scalar Ti = this->thermo().T()[celli];
-            const scalar pi = this->thermo().p()[celli];
-
-            scalarField c(nSpecie_, 0.0);
-            for (label i=0; i<nSpecie_; i++)
-            {
-                const scalar Yi = Y_[i][celli];
-                c[i] = rhoi*Yi/specieThermo_[i].W();
-            }
-
-            const scalarField dcdt(omega(c, Ti, pi));
-
-            for (label i=0; i<nSpecie_; i++)
-            {
-                RR_[i][celli] = dcdt[i]*specieThermo_[i].W();
-            }
+            RR_[i][celli] = dcdt[i]*specieThermo_[i].W();
         }
     }
 }
@@ -747,6 +752,11 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
 
     scalar deltaTMin = GREAT;
 
+    if (!this->chemistry_)
+    {
+        return deltaTMin;
+    }
+
     const volScalarField rho
     (
         IOobject
@@ -761,21 +771,6 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
         this->thermo().rho()
     );
 
-    if (this->mesh().changing())
-    {
-        for (label i = 0; i < nSpecie_; i++)
-        {
-            RR_[i].setSize(this->mesh().nCells());
-            RR_[i] = 0.0;
-        }
-    }
-
-    if (!this->chemistry_)
-    {
-        return deltaTMin;
-    }
-
-
     tmp<volScalarField> thc = this->thermo().hc();
     const scalarField& hc = thc();
 
@@ -856,7 +851,7 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
             "const scalar, "
             "const scalar, "
             "const scalar"
-        ")"
+        ") const"
     );
 
     return (0);
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.H
index 9128e41557..da560dae87 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.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
@@ -42,6 +42,7 @@ SourceFiles
 #include "ODE.H"
 #include "volFieldsFwd.H"
 #include "simpleMatrix.H"
+#include "DimensionedField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -90,14 +91,14 @@ protected:
         label nReaction_;
 
         //- List of reaction rate per specie [kg/m3/s]
-        PtrList<scalarField> RR_;
+        PtrList<DimensionedField<scalar, volMesh> > RR_;
 
 
     // Protected Member Functions
 
         //- Write access to chemical source terms
         //  (e.g. for multi-chemistry model)
-        inline PtrList<scalarField>& RR();
+        inline PtrList<DimensionedField<scalar, volMesh> >& RR();
 
 
 public:
@@ -205,8 +206,11 @@ public:
         // Chemistry model functions (overriding abstract functions in
         // basicChemistryModel.H)
 
-            //- Return const access to the chemical source terms
-            inline tmp<volScalarField> RR(const label i) const;
+            //- Return const access to the chemical source terms for specie, i
+            inline const DimensionedField<scalar, volMesh>& RR
+            (
+                const label i
+            ) const;
 
             //- Solve the reaction system for the given start time and time
             //  step and return the characteristic time
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModelI.H
index 8031b68ba4..9f06a064a7 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModelI.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModelI.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
@@ -29,7 +29,7 @@ License
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CompType, class ThermoType>
-inline Foam::PtrList<Foam::scalarField>&
+inline Foam::PtrList<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >&
 Foam::ODEChemistryModel<CompType, ThermoType>::RR()
 {
     return RR_;
@@ -69,36 +69,13 @@ Foam::ODEChemistryModel<CompType, ThermoType>::nReaction() const
 
 
 template<class CompType, class ThermoType>
-inline Foam::tmp<Foam::volScalarField>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
 Foam::ODEChemistryModel<CompType, ThermoType>::RR
 (
     const label i
 ) const
 {
-    tmp<volScalarField> tRR
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "RR(" + this->Y_[i].name() + ')',
-                this->time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            this->mesh(),
-            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-
-    if (this->chemistry_)
-    {
-        tRR().internalField() = RR_[i];
-        tRR().correctBoundaryConditions();
-    }
-    return tRR;
+    return RR_[i];
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C
index 38efa9e858..1182c22b92 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C
@@ -38,11 +38,7 @@ namespace Foam
 
 void Foam::basicChemistryModel::correct()
 {
-    if (mesh_.changing())
-    {
-        deltaTChem_.setSize(mesh_.nCells());
-        deltaTChem_ = deltaTChemIni_;
-    }
+    // do nothing
 }
 
 
@@ -64,7 +60,19 @@ Foam::basicChemistryModel::basicChemistryModel(const fvMesh& mesh)
     mesh_(mesh),
     chemistry_(lookup("chemistry")),
     deltaTChemIni_(readScalar(lookup("initialChemicalTimeStep"))),
-    deltaTChem_(mesh.nCells(), deltaTChemIni_)
+    deltaTChem_
+    (
+        IOobject
+        (
+            "deltaTChem",
+            mesh.time().constant(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("deltaTChem0", dimTime, deltaTChemIni_)
+    )
 {}
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
index 2f27d59f70..da54d4c916 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.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
@@ -40,6 +40,8 @@ SourceFiles
 #include "Switch.H"
 #include "scalarField.H"
 #include "volFieldsFwd.H"
+#include "volMesh.H"
+#include "DimensionedField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,14 +82,14 @@ protected:
         const scalar deltaTChemIni_;
 
         //- Latest estimation of integration step
-        scalarField deltaTChem_;
+        DimensionedField<scalar, volMesh> deltaTChem_;
 
 
     // Protected Member Functions
 
         //- Return non-const access to the latest estimation of integration
         //  step, e.g. for multi-chemistry model
-        scalarField& deltaTChem();
+        inline DimensionedField<scalar, volMesh>& deltaTChem();
 
         //- Correct function - updates due to mesh changes
         void correct();
@@ -118,7 +120,7 @@ public:
         inline Switch chemistry() const;
 
         //- Return the latest estimation of integration step
-        inline const scalarField& deltaTChem() const;
+        inline const DimensionedField<scalar, volMesh>& deltaTChem() const;
 
 
         // Functions to be derived in derived classes
@@ -126,7 +128,10 @@ public:
             // Fields
 
                 //- Return const access to chemical source terms [kg/m3/s]
-                virtual tmp<volScalarField> RR(const label i) const = 0;
+                virtual const DimensionedField<scalar, volMesh>& RR
+                (
+                    const label i
+                ) const = 0;
 
 
             // Chemistry solution
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModelI.H
index eeec0af749..6455309e15 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModelI.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModelI.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
@@ -37,13 +37,15 @@ inline Foam::Switch Foam::basicChemistryModel::chemistry() const
 }
 
 
-inline const Foam::scalarField& Foam::basicChemistryModel::deltaTChem() const
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::basicChemistryModel::deltaTChem() const
 {
     return deltaTChem_;
 }
 
 
-inline Foam::scalarField& Foam::basicChemistryModel::deltaTChem()
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::basicChemistryModel::deltaTChem()
 {
     return deltaTChem_;
 }
diff --git a/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C
index 093c6d8d84..23ad0e62a0 100644
--- a/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C
@@ -71,7 +71,19 @@ ODESolidChemistryModel
         RRs_.set
         (
             fieldI,
-            new scalarField(mesh.nCells(), 0.0)
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    "RRs::" + Ys_[fieldI].name(),
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+            )
         );
 
 
@@ -134,7 +146,8 @@ ODESolidChemistryModel
                     Y0Default
                 )
             );
-        // Calculate inital values of Ysi0 = rho*delta*Yi
+
+            // Calculate inital values of Ysi0 = rho*delta*Yi
             Ys0_[fieldI].internalField() =
                 this->solid().rho()
                *max(Ys_[fieldI], scalar(0.001))*mesh.V();
@@ -143,7 +156,23 @@ ODESolidChemistryModel
 
     forAll(RRg_, fieldI)
     {
-        RRg_.set(fieldI, new scalarField(mesh.nCells(), 0.0));
+        RRg_.set
+        (
+            fieldI,
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    "RRg::" + pyrolisisGases_[fieldI],
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+            )
+        );
     }
 
     forAll(gasThermo_, gasI)
@@ -519,9 +548,12 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::nEqns() const
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
-calculate()
+void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::calculate()
 {
+    if (!this->chemistry_)
+    {
+        return;
+    }
 
     const volScalarField rho
     (
@@ -537,58 +569,43 @@ calculate()
         this->solid().rho()
     );
 
-    if (this->mesh().changing())
-    {
-        forAll(RRs_, i)
-        {
-            RRs_[i].setSize(rho.size());
-        }
-        forAll(RRg_, i)
-        {
-            RRg_[i].setSize(rho.size());
-        }
-    }
-
     forAll(RRs_, i)
     {
-        RRs_[i] = 0.0;
+        RRs_[i].field() = 0.0;
     }
     forAll(RRg_, i)
     {
-        RRg_[i] = 0.0;
+        RRg_[i].field() = 0.0;
     }
 
-    if (this->chemistry_)
+    forAll(rho, celli)
     {
-        forAll(rho, celli)
+        cellCounter_ = celli;
+
+        const scalar delta = this->mesh().V()[celli];
+
+        if (reactingCells_[celli])
         {
-            cellCounter_ = celli;
+            scalar rhoi = rho[celli];
+            scalar Ti = this->solid().T()[celli];
+            scalar pi = this->solid().p()[celli];
 
-            const scalar delta = this->mesh().V()[celli];
-
-            if (reactingCells_[celli])
+            scalarField c(nSpecie_, 0.0);
+            for (label i=0; i<nSolids_; i++)
             {
-                scalar rhoi = rho[celli];
-                scalar Ti = this->solid().T()[celli];
-                scalar pi = this->solid().p()[celli];
+                c[i] = rhoi*Ys_[i][celli]*delta;
+            }
 
-                scalarField c(nSpecie_, 0.0);
-                for (label i=0; i<nSolids_; i++)
-                {
-                    c[i] = rhoi*Ys_[i][celli]*delta;
-                }
+            const scalarField dcdt = omega(c, Ti, pi, true);
 
-                const scalarField dcdt = omega(c, Ti, pi, true);
+            forAll(RRs_, i)
+            {
+                RRs_[i][celli] = dcdt[i]/delta;
+            }
 
-                forAll(RRs_, i)
-                {
-                    RRs_[i][celli] = dcdt[i]/delta;
-                }
-
-                forAll(RRg_, i)
-                {
-                    RRg_[i][celli] = dcdt[nSolids_ + i]/delta;
-                }
+            forAll(RRg_, i)
+            {
+                RRg_[i][celli] = dcdt[nSolids_ + i]/delta;
             }
         }
     }
@@ -603,6 +620,13 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
     const scalar deltaT
 )
 {
+    scalar deltaTMin = GREAT;
+
+    if (!this->chemistry_)
+    {
+        return deltaTMin;
+    }
+
     const volScalarField rho
     (
         IOobject
@@ -617,33 +641,15 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
         this->solid().rho()
     );
 
-    if (this->mesh().changing())
-    {
-        forAll(RRs_, i)
-        {
-            RRs_[i].setSize(rho.size());
-        }
-        forAll(RRg_, i)
-        {
-            RRg_[i].setSize(rho.size());
-        }
-    }
-
     forAll(RRs_, i)
     {
-        RRs_[i] = 0.0;
+        RRs_[i].field() = 0.0;
     }
     forAll(RRg_, i)
     {
-        RRg_[i] = 0.0;
+        RRg_[i].field() = 0.0;
     }
 
-    if (!this->chemistry_)
-    {
-        return GREAT;
-    }
-
-    scalar deltaTMin = GREAT;
 
     forAll(rho, celli)
     {
@@ -798,7 +804,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
             "const scalar, "
             "const scalar, "
             "const scalar"
-        ")"
+        ") const"
     );
     return (0);
 }
diff --git a/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.H
index 5dc785c9c7..df48750381 100644
--- a/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.H
+++ b/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.H
@@ -41,6 +41,7 @@ SourceFiles
 #include "solidReaction.H"
 #include "ODE.H"
 #include "volFieldsFwd.H"
+#include "DimensionedField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -96,19 +97,19 @@ protected:
         label nReaction_;
 
         //- List of reaction rate per solid [kg/m3/s]
-        PtrList<scalarField> RRs_;
+        PtrList<DimensionedField<scalar, volMesh> > RRs_;
 
         //- List of reaction rate per gas [kg/m3/s]
-        PtrList<scalarField> RRg_;
+        PtrList<DimensionedField<scalar, volMesh> > RRg_;
 
 
     // Protected Member Functions
 
         //- Write access to source terms for solids
-        inline PtrList<scalarField>& RRs();
+        inline PtrList<DimensionedField<scalar, volMesh> >& RRs();
 
         //- Write access to source terms for gases
-        inline PtrList<scalarField>& RRg();
+        inline PtrList<DimensionedField<scalar, volMesh> >& RRg();
 
 
 private:
@@ -203,19 +204,28 @@ public:
         // Chemistry model functions
 
             //- Return const access to the chemical source terms for solids
-            inline tmp<volScalarField> RRs(const label i) const;
+            inline const DimensionedField<scalar, volMesh>& RRs
+            (
+                const label i
+            ) const;
 
             //- Return const access to the chemical source terms for gases
-            inline tmp<volScalarField> RRg(const label i) const;
+            inline const DimensionedField<scalar, volMesh>& RRg
+            (
+                const label i
+            ) const;
 
             //- Return total gas source term
-            inline tmp<volScalarField> RRg() const;
+            inline tmp<DimensionedField<scalar, volMesh> > RRg() const;
 
             //- Return total solid source term
-            inline tmp<volScalarField> RRs() const;
+            inline tmp<DimensionedField<scalar, volMesh> > RRs() const;
 
             //- Return const access to the total source terms
-            inline tmp<volScalarField> RR(const label i) const;
+            inline const DimensionedField<scalar, volMesh>& RR
+            (
+                const label i
+            ) const;
 
             //- Return sensible enthalpy for gas i [J/Kg]
             virtual tmp<volScalarField> gasHs
diff --git a/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModelI.H b/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModelI.H
index 802345f3a9..bcb6a89f99 100644
--- a/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModelI.H
+++ b/src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModelI.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
@@ -29,7 +29,7 @@ License
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CompType, class SolidThermo, class GasThermo>
-inline Foam::PtrList<Foam::scalarField>&
+inline Foam::PtrList<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >&
 Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRs()
 {
     return RRs_;
@@ -37,7 +37,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRs()
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-inline Foam::PtrList<Foam::scalarField>&
+inline Foam::PtrList<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >&
 Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg()
 {
     return RRg_;
@@ -87,80 +87,34 @@ nReaction() const
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-inline Foam::tmp<Foam::volScalarField>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
 Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRs
 (
     const label i
 ) const
 {
-    tmp<volScalarField> tRRs
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "RRs(" + Ys_[i].name() + ')',
-                this->time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            this->mesh(),
-            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-
-    if (this->chemistry_)
-    {
-        tRRs().internalField() = RRs_[i];
-        tRRs().correctBoundaryConditions();
-    }
-    return tRRs;
+    return RRs_[i];
 }
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-inline Foam::tmp<Foam::volScalarField>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
 Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg
 (
     const label i
 ) const
 {
-    tmp<volScalarField> tRRg
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "RRg(" + this->pyrolisisGases_[i] + ')',
-                this->time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            this->mesh(),
-            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-
-    if (this->chemistry_)
-    {
-        tRRg().internalField() = RRg_[i];
-        tRRg().correctBoundaryConditions();
-    }
-    return tRRg;
+    return RRg_[i];
 }
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-inline Foam::tmp<Foam::volScalarField>
+inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
 Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg() const
 {
-    tmp<volScalarField> tRRg
+    tmp<DimensionedField<scalar, volMesh> > tRRg
     (
-        new volScalarField
+        new DimensionedField<scalar, volMesh>
         (
             IOobject
             (
@@ -171,30 +125,29 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg() const
                 IOobject::NO_WRITE
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0),
-            zeroGradientFvPatchScalarField::typeName
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
         )
     );
 
     if (this->chemistry_)
     {
+        DimensionedField<scalar, volMesh>& RRg = tRRg();
         for (label i=0; i < nGases_; i++)
         {
-            tRRg().internalField() += RRg_[i];
+            RRg += RRg_[i];
         }
-        tRRg().correctBoundaryConditions();
     }
     return tRRg;
 }
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-inline Foam::tmp<Foam::volScalarField>
+inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
 Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRs() const
 {
-    tmp<volScalarField> tRRs
+    tmp<DimensionedField<scalar, volMesh> > tRRs
     (
-        new volScalarField
+        new DimensionedField<scalar, volMesh>
         (
             IOobject
             (
@@ -205,32 +158,31 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRs() const
                 IOobject::NO_WRITE
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0),
-            zeroGradientFvPatchScalarField::typeName
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
         )
     );
 
     if (this->chemistry_)
     {
+        DimensionedField<scalar, volMesh>& RRs = tRRs();
         for (label i=0; i < nSolids_; i++)
         {
-            tRRs().internalField() += RRs_[i];
+            RRs += RRs_[i];
         }
-        tRRs().correctBoundaryConditions();
     }
     return tRRs;
 }
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-inline Foam::tmp<Foam::volScalarField>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
 Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RR
 (
     const label i
 ) const
 {
     notImplemented("ODESolidChemistryModel::RR(const label)");
-    return (Foam::volScalarField::null());
+    return (DimensionedField<scalar, volMesh>::null());
 }
 
 
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
index 9114f046ff..d9580162ef 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
@@ -124,16 +124,22 @@ public:
         inline const solidReactionThermo& solid() const;
 
         //- Return total gases mass source term [kg/m3/s]
-        virtual tmp<volScalarField> RRg() const = 0;
+        virtual tmp<DimensionedField<scalar, volMesh> > RRg() const = 0;
 
         //- Return total solids mass source term [kg/m3/s]
-        virtual tmp<volScalarField> RRs() const = 0;
+        virtual tmp<DimensionedField<scalar, volMesh> > RRs() const = 0;
 
         //- Return chemical source terms for solids [kg/m3/s]
-        virtual tmp<volScalarField> RRs(const label i) const = 0;
+        virtual const DimensionedField<scalar, volMesh>& RRs
+        (
+            const label i
+        ) const = 0;
 
         //- Return chemical source terms for gases [kg/m3/s]
-        virtual tmp<volScalarField> RRg(const label i) const = 0;
+        virtual const DimensionedField<scalar, volMesh>& RRg
+        (
+            const label i
+        ) const = 0;
 
         //- Return sensible enthalpy for gas i [J/Kg]
         virtual tmp<volScalarField> gasHs
diff --git a/tutorials/combustion/chemFoam/nc7h16/validation/createGraph b/tutorials/combustion/chemFoam/nc7h16/validation/createGraph
old mode 100755
new mode 100644
index 634beebec7..2302a2fe42
--- a/tutorials/combustion/chemFoam/nc7h16/validation/createGraph
+++ b/tutorials/combustion/chemFoam/nc7h16/validation/createGraph
@@ -11,10 +11,10 @@ gnuplot<<EOF
     set xlabel "Time / [s]" font "Helvetica,24"
     set ylabel "Temperature / [K]" font "Helvetica,24"
     set grid
-    set key left top 
+    set key left top
     set xrange [0:0.001]
     set yrange [750:2750]
-    set ytics 250 
+    set ytics 250
     plot \
         "../chemFoam.out" u 1:2 t "OpenFOAM" with points lt 1 pt 6 ps 1.5,\
         "chemkinII" with lines title "Chemkin II" lt -1
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties
index f1fec22f39..91b20eee70 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties
@@ -53,13 +53,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              300;
     Cp0             4187;
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties
index 9973e9ad7b..1576592b91 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/combustionProperties
@@ -23,7 +23,8 @@ active  true;
 
 infinitelyFastChemistryCoeffs
 {
-    C       5.0;
+    semiImplicit no;
+    C           5.0;
 }
 
 FSDCoeffs
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties
index f1fec22f39..91b20eee70 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties
@@ -53,13 +53,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              300;
     Cp0             4187;
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/combustionProperties b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/combustionProperties
index cfa2cd8f28..689dc54c62 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/combustionProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/combustionProperties
@@ -22,6 +22,7 @@ active  on;
 
 infinitelyFastChemistryCoeffs
 {
+    semiImplicit no;
     C       5.0;
 }
 
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties
index f1fec22f39..91b20eee70 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties
@@ -53,13 +53,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              300;
     Cp0             4187;
diff --git a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties b/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
index c155594c01..5c93f36555 100644
--- a/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
@@ -60,13 +60,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              350;
     Cp0             4100;
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
index 8f07ac40e7..fc5042f2c1 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
@@ -54,13 +54,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              300;
     Cp0             4187;
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
index f149873427..29cfd8d563 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
@@ -53,10 +53,6 @@ constantProperties
 {
     parcelTypeId    2;
 
-    rhoMin          1e-15;
-    TMin            200;
-    minParticleMass 1e-15;
-
     rho0            2500;
     T0              300;
     Cp0             900;
diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
index 3ea1510347..3601470655 100644
--- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
+++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
@@ -45,11 +45,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    minParticleMass 1e-15;
-
     rho0            964;
     youngsModulus   6e8;
     poissonsRatio   0.35;
diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
index fac715eceb..fbd8c9f174 100644
--- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
+++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
@@ -37,11 +37,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    minParticleMass 1e-15;
-
     rho0            964;
     youngsModulus   6e8;
     poissonsRatio   0.35;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties
index 36af7910eb..613b407de1 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties
@@ -54,13 +54,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              300;
     Cp0             4100;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactingCloud1Properties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
index c2fb9a550b..10bb7ccc33 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
@@ -54,13 +54,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              350;
     Cp0             4100;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
index c05dc6a728..aac578722d 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
@@ -54,13 +54,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              350;
     Cp0             4100;
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties
index 4b520415a0..875e31a59c 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties
@@ -53,13 +53,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              300;
     Cp0             4187;
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties
index adef8f34d4..78f4e7bb5d 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties
@@ -53,13 +53,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              300;
     Cp0             4187;
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
index 6c2f59c4b4..6704aad914 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
@@ -53,13 +53,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              300;
     Cp0             4187;
diff --git a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties
index 778983882a..79ede769ec 100644
--- a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties
@@ -54,13 +54,6 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     rho0            1000;
     T0              300;
     Cp0             4187;
diff --git a/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties b/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties
index 987dd33eaa..153a5fc0b2 100644
--- a/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties
+++ b/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties
@@ -54,16 +54,9 @@ solution
 
 constantProperties
 {
-    parcelTypeId    1;
-
-    rhoMin          1e-15;
-    TMin            200;
-    pMin            1000;
-    minParticleMass 1e-15;
-
     T0              320;
 
-    // place holders for rho0 and T0 - reset from liquid props using T0
+    // place holders for rho0 and Cp0 - reset from liquid props using T0
     rho0            1000;
     Cp0             4187;
 
diff --git a/tutorials/mesh/cvMesh/blob/system/cvMeshDict b/tutorials/mesh/cvMesh/blob/system/cvMeshDict
index 27fff06883..af0db03508 100644
--- a/tutorials/mesh/cvMesh/blob/system/cvMeshDict
+++ b/tutorials/mesh/cvMesh/blob/system/cvMeshDict
@@ -151,7 +151,7 @@ motionControl
         {
             priority            1;
             mode                bothSides;
-            
+
             surfaceCellSizeFunction uniformValue;
             uniformValueCoeffs
             {
diff --git a/tutorials/mesh/cvMesh/simpleShapes/system/cvMeshDict b/tutorials/mesh/cvMesh/simpleShapes/system/cvMeshDict
index 2d7c9faed2..26952c55ce 100644
--- a/tutorials/mesh/cvMesh/simpleShapes/system/cvMeshDict
+++ b/tutorials/mesh/cvMesh/simpleShapes/system/cvMeshDict
@@ -139,7 +139,7 @@ surfaceConformation
             // conformation, initial protrusion tests.
             surfacePtReplaceDistCoeff   0.5;
 
-            
+
             surfacePtExclusionDistanceCoeff 0.5;
         }
 
diff --git a/tutorials/mesh/snappyHexMesh/Allrun b/tutorials/mesh/snappyHexMesh/Allrun
index 057d8bf8dd..6d1366f448 100755
--- a/tutorials/mesh/snappyHexMesh/Allrun
+++ b/tutorials/mesh/snappyHexMesh/Allrun
@@ -1,12 +1,12 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # run from this directory
 
-                                                                                                                                                                                                                                                                                                                          
-(                                                                                                                                                                                                                                                                                                                         
-    cd flange || exit                                                                                                                                                                                                                                                                                                   
-    ./Allrun                                                                                                                                                                                                                                                                                                              
-)                                                                                                                                                                                                                                                                                                                         
-                                                                                                                                                                                                                                                                                                                          
+
+(
+    cd flange || exit
+    ./Allrun
+)
+
 exit 0
 
 # These cases are links to solver test cases and are run when the Allrun