diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
index 8a7fb290b8..f6794e92b2 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@@ -38,7 +38,6 @@ Description
 #include "CoalCloud.H"
 #include "psiChemistryModel.H"
 #include "chemistrySolver.H"
-#include "thermoPhysicsTypes.H"
 #include "timeActivatedExplicitCellSource.H"
 #include "radiationModel.H"
 
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H b/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H
index 0a5fd91b75..dbd7cf9659 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/createClouds.H
@@ -1,5 +1,5 @@
 Info<< "\nConstructing coal cloud" << endl;
-CoalCloud<gasThermoPhysics> coalParcels
+thermoCoalCloud coalParcels
 (
     "coalCloud1",
     rho,
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H
index 728e605e57..2accb8e1c9 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H
@@ -1,5 +1,5 @@
 Info<< "\nConstructing reacting cloud" << endl;
-BasicReactingCloud<icoPoly8ThermoPhysics> parcels
+icoPoly8ThermoReactingCloud parcels
 (
     "reactingCloud1",
     rho,
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
index 4242452440..bf0a16af28 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
@@ -46,7 +46,6 @@ Description
 #include "BasicReactingCloud.H"
 #include "rhoChemistryModel.H"
 #include "chemistrySolver.H"
-#include "thermoPhysicsTypes.H"
 #include "radiationModel.H"
 #include "porousZones.H"
 #include "timeActivatedExplicitMulticomponentPointSource.H"
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
index 4ae0633ab7..f3e043143c 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
@@ -1,5 +1,5 @@
 Info<< "\nConstructing reacting cloud" << endl;
-BasicReactingCloud<gasThermoPhysics> parcels
+thermoReactingCloud parcels
 (
     "reactingCloud1",
     rho,
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index 62f8f4834a..4dfd367d09 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -37,7 +37,6 @@ Description
 #include "BasicReactingCloud.H"
 #include "psiChemistryModel.H"
 #include "chemistrySolver.H"
-#include "thermoPhysicsTypes.H"
 #include "radiationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index cd7c13d783..92c92d72e5 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -280,11 +280,14 @@ addLayersControls
     nBufferCellsNoExtrude 0;
 
 
-    // Overall max number of layer addition iterations
+    // Overall max number of layer addition iterations. The mesher will exit
+    // if it reaches this number of iterations; possibly with an illegal
+    // mesh.
     nLayerIter 50;
 
     // Max number of iterations after which relaxed meshQuality controls
-    // get used.
+    // get used. Up to nRelaxIter it uses the settings in meshQualityControls,
+    // after nRelaxIter it uses the values in meshQualityControls::relaxed.
     nRelaxedIter 20;
 }
 
diff --git a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C
index 3d23cce5c3..c6fa201218 100644
--- a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C
+++ b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C
@@ -125,6 +125,7 @@ label findPatchID(const polyMesh& mesh, const word& name)
 
 int main(int argc, char *argv[])
 {
+#   include "addRegionOption.H"
     argList::validArgs.append("faceZone");
     argList::validArgs.append("patch");
     argList::validOptions.insert("additionalPatches", "(patch2 .. patchN)");
@@ -134,7 +135,7 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
     runTime.functionObjects().off();
-#   include "createMesh.H"
+#   include "createNamedMesh.H"
     const word oldInstance = mesh.pointsInstance();
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatch.C b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
index 40ba0efade..6374722e46 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatch.C
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
@@ -516,11 +516,13 @@ void syncPoints
 
 int main(int argc, char *argv[])
 {
+#   include "addRegionOption.H"
     argList::validOptions.insert("overwrite", "");
 
 #   include "setRootCase.H"
 #   include "createTime.H"
     runTime.functionObjects().off();
+#   include "createNamedPolyMesh.H"
 
     const bool overwrite = args.optionFound("overwrite");
 
@@ -532,6 +534,11 @@ int main(int argc, char *argv[])
         (
             "createPatchDict",
             runTime.system(),
+            (
+                regionName != polyMesh::defaultRegion
+              ? regionName
+              : word::null
+            ),
             runTime,
             IOobject::MUST_READ,
             IOobject::NO_WRITE,
@@ -551,7 +558,6 @@ int main(int argc, char *argv[])
     coupledPolyPatch::matchTol = tol;
 
 
-#   include "createPolyMesh.H"
     const word oldInstance = mesh.pointsInstance();
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
diff --git a/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C b/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C
index e12ab576fd..b1c064a6ee 100644
--- a/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C
+++ b/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C
@@ -222,13 +222,14 @@ labelList findBaffles(const polyMesh& mesh, const labelList& boundaryFaces)
 
 int main(int argc, char *argv[])
 {
+#   include "addRegionOption.H"
     argList::validOptions.insert("split", "");
     argList::validOptions.insert("overwrite", "");
     argList::validOptions.insert("detectOnly", "");
 #   include "setRootCase.H"
 #   include "createTime.H"
     runTime.functionObjects().off();
-#   include "createMesh.H"
+#   include "createNamedMesh.H"
     const word oldInstance = mesh.pointsInstance();
 
     bool split      = args.optionFound("split");
diff --git a/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C b/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C
index a369a195b8..df419e118b 100644
--- a/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C
+++ b/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C
@@ -60,6 +60,7 @@ int main(int argc, char *argv[])
 {
     argList::validOptions.insert("noFlipMap", "");
 
+#   include "addRegionOption.H"
 #   include "addTimeOptions.H"
 #   include "setRootCase.H"
 #   include "createTime.H"
@@ -77,7 +78,7 @@ int main(int argc, char *argv[])
 
     runTime.setTime(Times[startTime], startTime);
 
-#   include "createPolyMesh.H"
+#   include "createNamedPolyMesh.H"
 
     // Search for list of objects for the time of the mesh
     IOobjectList objects
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
index aad15ee459..e1620d7201 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict
@@ -22,6 +22,12 @@ numberOfSubdomains  4;
 //- Keep owner and neighbour on same processor for faces in zones:
 // preserveFaceZones (heater solid1 solid3);
 
+
+//- Keep owner and neighbour on same processor for faces in patches:
+//  (makes sense only for cyclic patches)
+//preservePatches (cyclic_left_right);
+
+
 method          scotch;
 // method          hierarchical;
 // method          simple;
diff --git a/applications/utilities/parallelProcessing/decomposePar/distributeCells.C b/applications/utilities/parallelProcessing/decomposePar/distributeCells.C
index af47af244f..ffd54ac685 100644
--- a/applications/utilities/parallelProcessing/decomposePar/distributeCells.C
+++ b/applications/utilities/parallelProcessing/decomposePar/distributeCells.C
@@ -45,6 +45,35 @@ void domainDecomposition::distributeCells()
 
     labelHashSet sameProcFaces;
 
+    if (decompositionDict_.found("preservePatches"))
+    {
+        wordList pNames(decompositionDict_.lookup("preservePatches"));
+
+        Info<< "Keeping owner of faces in patches " << pNames
+            << " on same processor. This only makes sense for cyclics." << endl;
+
+        const polyBoundaryMesh& patches = boundaryMesh();
+
+        forAll(pNames, i)
+        {
+            label patchI = patches.findPatchID(pNames[i]);
+
+            if (patchI == -1)
+            {
+                FatalErrorIn("domainDecomposition::distributeCells()")
+                    << "Unknown preservePatch " << pNames[i]
+                    << endl << "Valid patches are " << patches.names()
+                    << exit(FatalError);
+            }
+
+            const polyPatch& pp = patches[patchI];
+
+            forAll(pp, i)
+            {
+                sameProcFaces.insert(pp.start() + i);
+            }
+        }
+    }
     if (decompositionDict_.found("preserveFaceZones"))
     {
         wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
diff --git a/src/OpenFOAM/db/Time/Time.C b/src/OpenFOAM/db/Time/Time.C
index f0c67306be..7405cfb5b6 100644
--- a/src/OpenFOAM/db/Time/Time.C
+++ b/src/OpenFOAM/db/Time/Time.C
@@ -115,7 +115,14 @@ void Foam::Time::setControls()
         {
             if (timeDirs.size())
             {
-                startTime_ = timeDirs[0].value();
+                if (timeDirs[0].name() == constant() && timeDirs.size() >= 2)
+                {
+                    startTime_ = timeDirs[1].value();
+                }
+                else
+                {
+                    startTime_ = timeDirs[0].value();
+                }
             }
         }
         else if (startFrom == "latestTime")
diff --git a/src/OpenFOAM/dimensionSet/dimensionSet.C b/src/OpenFOAM/dimensionSet/dimensionSet.C
index 9fba472d98..e60a3998d4 100644
--- a/src/OpenFOAM/dimensionSet/dimensionSet.C
+++ b/src/OpenFOAM/dimensionSet/dimensionSet.C
@@ -110,6 +110,7 @@ Foam::scalar Foam::dimensionSet::operator[](const dimensionType type) const
     return exponents_[type];
 }
 
+
 Foam::scalar& Foam::dimensionSet::operator[](const dimensionType type)
 {
     return exponents_[type];
@@ -130,6 +131,7 @@ bool Foam::dimensionSet::operator==(const dimensionSet& ds) const
     return equall;
 }
 
+
 bool Foam::dimensionSet::operator!=(const dimensionSet& ds) const
 {
     return !(operator==(ds));
@@ -163,6 +165,7 @@ bool Foam::dimensionSet::operator+=(const dimensionSet& ds) const
     return true;
 }
 
+
 bool Foam::dimensionSet::operator-=(const dimensionSet& ds) const
 {
     if (dimensionSet::debug && *this != ds)
@@ -176,6 +179,7 @@ bool Foam::dimensionSet::operator-=(const dimensionSet& ds) const
     return true;
 }
 
+
 bool Foam::dimensionSet::operator*=(const dimensionSet& ds)
 {
     reset((*this)*ds);
@@ -183,6 +187,7 @@ bool Foam::dimensionSet::operator*=(const dimensionSet& ds)
     return true;
 }
 
+
 bool Foam::dimensionSet::operator/=(const dimensionSet& ds)
 {
     reset((*this)/ds);
@@ -206,6 +211,7 @@ Foam::dimensionSet Foam::max(const dimensionSet& ds1, const dimensionSet& ds2)
     return ds1;
 }
 
+
 Foam::dimensionSet Foam::min(const dimensionSet& ds1, const dimensionSet& ds2)
 {
     if (dimensionSet::debug && ds1 != ds2)
@@ -256,6 +262,7 @@ Foam::dimensionSet Foam::pow(const dimensionSet& ds, const scalar p)
     return dimPow;
 }
 
+
 Foam::dimensionSet Foam::pow
 (
     const dimensionSet& ds,
@@ -283,6 +290,7 @@ Foam::dimensionSet Foam::pow
     return dimPow;
 }
 
+
 Foam::dimensionSet Foam::pow
 (
     const dimensionedScalar& dS,
@@ -309,61 +317,79 @@ Foam::dimensionSet Foam::sqr(const dimensionSet& ds)
     return pow(ds, 2);
 }
 
+
 Foam::dimensionSet Foam::pow3(const dimensionSet& ds)
 {
     return pow(ds, 3);
 }
 
+
 Foam::dimensionSet Foam::pow4(const dimensionSet& ds)
 {
     return pow(ds, 4);
 }
 
+
 Foam::dimensionSet Foam::pow5(const dimensionSet& ds)
 {
     return pow(ds, 5);
 }
 
+
 Foam::dimensionSet Foam::pow6(const dimensionSet& ds)
 {
     return pow(ds, 6);
 }
 
+
+Foam::dimensionSet Foam::pow025(const dimensionSet& ds)
+{
+    return sqrt(sqrt(ds));
+}
+
+
 Foam::dimensionSet Foam::sqrt(const dimensionSet& ds)
 {
     return pow(ds, 0.5);
 }
 
+
 Foam::dimensionSet Foam::magSqr(const dimensionSet& ds)
 {
     return pow(ds, 2);
 }
 
+
 Foam::dimensionSet Foam::mag(const dimensionSet& ds)
 {
     return ds;
 }
 
+
 Foam::dimensionSet Foam::sign(const dimensionSet&)
 {
     return dimless;
 }
 
+
 Foam::dimensionSet Foam::pos(const dimensionSet&)
 {
     return dimless;
 }
 
+
 Foam::dimensionSet Foam::neg(const dimensionSet&)
 {
     return dimless;
 }
 
+
 Foam::dimensionSet Foam::inv(const dimensionSet& ds)
 {
     return dimless/ds;
 }
 
+
 Foam::dimensionSet Foam::trans(const dimensionSet& ds)
 {
     if (dimensionSet::debug && !ds.dimensionless())
@@ -376,6 +402,7 @@ Foam::dimensionSet Foam::trans(const dimensionSet& ds)
     return ds;
 }
 
+
 Foam::dimensionSet Foam::transform(const dimensionSet& ds)
 {
     return ds;
@@ -389,6 +416,7 @@ Foam::dimensionSet Foam::operator-(const dimensionSet& ds)
     return ds;
 }
 
+
 Foam::dimensionSet Foam::operator+
 (
     const dimensionSet& ds1,
@@ -409,6 +437,7 @@ Foam::dimensionSet Foam::operator+
     return dimSum;
 }
 
+
 Foam::dimensionSet Foam::operator-
 (
     const dimensionSet& ds1,
@@ -429,6 +458,7 @@ Foam::dimensionSet Foam::operator-
     return dimDifference;
 }
 
+
 Foam::dimensionSet Foam::operator*
 (
     const dimensionSet& ds1,
@@ -445,6 +475,7 @@ Foam::dimensionSet Foam::operator*
     return dimProduct;
 }
 
+
 Foam::dimensionSet Foam::operator/
 (
     const dimensionSet& ds1,
@@ -471,6 +502,7 @@ Foam::dimensionSet Foam::operator&
     return ds1*ds2;
 }
 
+
 Foam::dimensionSet Foam::operator^
 (
     const dimensionSet& ds1,
@@ -480,6 +512,7 @@ Foam::dimensionSet Foam::operator^
     return ds1*ds2;
 }
 
+
 Foam::dimensionSet Foam::operator&&
 (
     const dimensionSet& ds1,
diff --git a/src/OpenFOAM/dimensionSet/dimensionSet.H b/src/OpenFOAM/dimensionSet/dimensionSet.H
index e61d3ec4df..13b0a899a9 100644
--- a/src/OpenFOAM/dimensionSet/dimensionSet.H
+++ b/src/OpenFOAM/dimensionSet/dimensionSet.H
@@ -70,6 +70,7 @@ dimensionSet pow3(const dimensionSet&);
 dimensionSet pow4(const dimensionSet&);
 dimensionSet pow5(const dimensionSet&);
 dimensionSet pow6(const dimensionSet&);
+dimensionSet pow025(const dimensionSet&);
 
 dimensionSet sqrt(const dimensionSet&);
 dimensionSet magSqr(const dimensionSet&);
@@ -226,6 +227,7 @@ public:
         friend dimensionSet pow4(const dimensionSet&);
         friend dimensionSet pow5(const dimensionSet&);
         friend dimensionSet pow6(const dimensionSet&);
+        friend dimensionSet pow025(const dimensionSet&);
 
         friend dimensionSet sqrt(const dimensionSet&);
         friend dimensionSet magSqr(const dimensionSet&);
diff --git a/src/OpenFOAM/dimensionedTypes/dimensionedScalar/dimensionedScalar.C b/src/OpenFOAM/dimensionedTypes/dimensionedScalar/dimensionedScalar.C
index 8761396a0f..054d1eb43b 100644
--- a/src/OpenFOAM/dimensionedTypes/dimensionedScalar/dimensionedScalar.C
+++ b/src/OpenFOAM/dimensionedTypes/dimensionedScalar/dimensionedScalar.C
@@ -38,32 +38,38 @@ dimensionedScalar operator+(const dimensionedScalar& ds1, const scalar s2)
     return ds1 + dimensionedScalar(s2);
 }
 
+
 dimensionedScalar operator+(const scalar s1, const dimensionedScalar& ds2)
 {
     return dimensionedScalar(s1) + ds2;
 }
 
+
 dimensionedScalar operator-(const dimensionedScalar& ds1, const scalar s2)
 {
     return ds1 - dimensionedScalar(s2);
 }
 
+
 dimensionedScalar operator-(const scalar s1, const dimensionedScalar& ds2)
 {
     return dimensionedScalar(s1) - ds2;
 }
 
+
 dimensionedScalar operator*(const dimensionedScalar& ds1, const scalar s2)
 {
     return ds1 * dimensionedScalar(s2);
 }
 
+
 dimensionedScalar operator/(const scalar s1, const dimensionedScalar& ds2)
 {
     return dimensionedScalar(s1)/ds2;
 }
 
 
+
 dimensionedScalar pow
 (
     const dimensionedScalar& ds,
@@ -78,6 +84,7 @@ dimensionedScalar pow
     );
 }
 
+
 dimensionedScalar pow3(const dimensionedScalar& ds)
 {
     return dimensionedScalar
@@ -88,6 +95,7 @@ dimensionedScalar pow3(const dimensionedScalar& ds)
     );
 }
 
+
 dimensionedScalar pow4(const dimensionedScalar& ds)
 {
     return dimensionedScalar
@@ -98,6 +106,7 @@ dimensionedScalar pow4(const dimensionedScalar& ds)
     );
 }
 
+
 dimensionedScalar pow5(const dimensionedScalar& ds)
 {
     return dimensionedScalar
@@ -108,6 +117,7 @@ dimensionedScalar pow5(const dimensionedScalar& ds)
     );
 }
 
+
 dimensionedScalar pow6(const dimensionedScalar& ds)
 {
     return dimensionedScalar
@@ -118,6 +128,18 @@ dimensionedScalar pow6(const dimensionedScalar& ds)
     );
 }
 
+
+dimensionedScalar pow025(const dimensionedScalar& ds)
+{
+    return dimensionedScalar
+    (
+        "pow025(" + ds.name() + ')',
+        pow025(ds.dimensions()),
+        pow025(ds.value())
+    );
+}
+
+
 dimensionedScalar sqrt(const dimensionedScalar& ds)
 {
     return dimensionedScalar
@@ -128,6 +150,7 @@ dimensionedScalar sqrt(const dimensionedScalar& ds)
     );
 }
 
+
 dimensionedScalar cbrt(const dimensionedScalar& ds)
 {
     return dimensionedScalar
@@ -138,6 +161,7 @@ dimensionedScalar cbrt(const dimensionedScalar& ds)
     );
 }
 
+
 dimensionedScalar hypot
 (
     const dimensionedScalar& x,
@@ -152,6 +176,7 @@ dimensionedScalar hypot
     );
 }
 
+
 dimensionedScalar sign(const dimensionedScalar& ds)
 {
     return dimensionedScalar
@@ -162,6 +187,7 @@ dimensionedScalar sign(const dimensionedScalar& ds)
     );
 }
 
+
 dimensionedScalar pos(const dimensionedScalar& ds)
 {
     return dimensionedScalar
@@ -172,6 +198,7 @@ dimensionedScalar pos(const dimensionedScalar& ds)
     );
 }
 
+
 dimensionedScalar neg(const dimensionedScalar& ds)
 {
     return dimensionedScalar
diff --git a/src/OpenFOAM/dimensionedTypes/dimensionedScalar/dimensionedScalar.H b/src/OpenFOAM/dimensionedTypes/dimensionedScalar/dimensionedScalar.H
index 492d5064ef..c6e46168c9 100644
--- a/src/OpenFOAM/dimensionedTypes/dimensionedScalar/dimensionedScalar.H
+++ b/src/OpenFOAM/dimensionedTypes/dimensionedScalar/dimensionedScalar.H
@@ -61,6 +61,7 @@ dimensionedScalar pow3(const dimensionedScalar&);
 dimensionedScalar pow4(const dimensionedScalar&);
 dimensionedScalar pow5(const dimensionedScalar&);
 dimensionedScalar pow6(const dimensionedScalar&);
+dimensionedScalar pow025(const dimensionedScalar&);
 
 dimensionedScalar sqrt(const dimensionedScalar&);
 dimensionedScalar cbrt(const dimensionedScalar&);
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.C
index ade367e8ff..4aad204f79 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.C
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.C
@@ -376,6 +376,7 @@ UNARY_FUNCTION(scalar, scalar, pow3, pow3)
 UNARY_FUNCTION(scalar, scalar, pow4, pow4)
 UNARY_FUNCTION(scalar, scalar, pow5, pow5)
 UNARY_FUNCTION(scalar, scalar, pow6, pow6)
+UNARY_FUNCTION(scalar, scalar, pow025, pow025)
 UNARY_FUNCTION(scalar, scalar, sqrt, sqrt)
 UNARY_FUNCTION(scalar, scalar, sign, sign)
 UNARY_FUNCTION(scalar, scalar, pos, pos)
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.H b/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.H
index 1b5e5d22bf..e2d565283e 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.H
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.H
@@ -84,6 +84,7 @@ UNARY_FUNCTION(scalar, scalar, pow3, pow3)
 UNARY_FUNCTION(scalar, scalar, pow4, pow4)
 UNARY_FUNCTION(scalar, scalar, pow5, pow5)
 UNARY_FUNCTION(scalar, scalar, pow6, pow6)
+UNARY_FUNCTION(scalar, scalar, pow025, pow025)
 UNARY_FUNCTION(scalar, scalar, sqrt, sqrt)
 UNARY_FUNCTION(scalar, scalar, sign, sign)
 UNARY_FUNCTION(scalar, scalar, pos, pos)
diff --git a/src/OpenFOAM/fields/FieldFields/scalarFieldField/scalarFieldField.C b/src/OpenFOAM/fields/FieldFields/scalarFieldField/scalarFieldField.C
index fe9067a8f5..8a308adbb5 100644
--- a/src/OpenFOAM/fields/FieldFields/scalarFieldField/scalarFieldField.C
+++ b/src/OpenFOAM/fields/FieldFields/scalarFieldField/scalarFieldField.C
@@ -104,6 +104,7 @@ UNARY_FUNCTION(scalar, scalar, pow3)
 UNARY_FUNCTION(scalar, scalar, pow4)
 UNARY_FUNCTION(scalar, scalar, pow5)
 UNARY_FUNCTION(scalar, scalar, pow6)
+UNARY_FUNCTION(scalar, scalar, pow025)
 UNARY_FUNCTION(scalar, scalar, sqrt)
 UNARY_FUNCTION(scalar, scalar, sign)
 UNARY_FUNCTION(scalar, scalar, pos)
diff --git a/src/OpenFOAM/fields/FieldFields/scalarFieldField/scalarFieldField.H b/src/OpenFOAM/fields/FieldFields/scalarFieldField/scalarFieldField.H
index 45f1e86ea0..019b731979 100644
--- a/src/OpenFOAM/fields/FieldFields/scalarFieldField/scalarFieldField.H
+++ b/src/OpenFOAM/fields/FieldFields/scalarFieldField/scalarFieldField.H
@@ -57,6 +57,7 @@ void stabilise
     const scalar s
 );
 
+
 template<template<class> class Field>
 tmp<FieldField<Field, scalar> > stabilise
 (
@@ -64,6 +65,7 @@ tmp<FieldField<Field, scalar> > stabilise
     const scalar s
 );
 
+
 template<template<class> class Field>
 tmp<FieldField<Field, scalar> > stabilise
 (
@@ -95,6 +97,7 @@ UNARY_FUNCTION(scalar, scalar, pow3)
 UNARY_FUNCTION(scalar, scalar, pow4)
 UNARY_FUNCTION(scalar, scalar, pow5)
 UNARY_FUNCTION(scalar, scalar, pow6)
+UNARY_FUNCTION(scalar, scalar, pow025)
 UNARY_FUNCTION(scalar, scalar, sqrt)
 UNARY_FUNCTION(scalar, scalar, sign)
 UNARY_FUNCTION(scalar, scalar, pos)
diff --git a/src/OpenFOAM/fields/Fields/scalarField/scalarField.C b/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
index 4f2eccfa9c..ed8e2da5e3 100644
--- a/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
+++ b/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
@@ -109,6 +109,7 @@ UNARY_FUNCTION(scalar, scalar, pow3)
 UNARY_FUNCTION(scalar, scalar, pow4)
 UNARY_FUNCTION(scalar, scalar, pow5)
 UNARY_FUNCTION(scalar, scalar, pow6)
+UNARY_FUNCTION(scalar, scalar, pow025)
 UNARY_FUNCTION(scalar, scalar, sqrt)
 UNARY_FUNCTION(scalar, scalar, sign)
 UNARY_FUNCTION(scalar, scalar, pos)
diff --git a/src/OpenFOAM/fields/Fields/scalarField/scalarField.H b/src/OpenFOAM/fields/Fields/scalarField/scalarField.H
index a1d0865d20..bd25d64a42 100644
--- a/src/OpenFOAM/fields/Fields/scalarField/scalarField.H
+++ b/src/OpenFOAM/fields/Fields/scalarField/scalarField.H
@@ -96,6 +96,7 @@ UNARY_FUNCTION(scalar, scalar, pow3)
 UNARY_FUNCTION(scalar, scalar, pow4)
 UNARY_FUNCTION(scalar, scalar, pow5)
 UNARY_FUNCTION(scalar, scalar, pow6)
+UNARY_FUNCTION(scalar, scalar, pow025)
 UNARY_FUNCTION(scalar, scalar, sqrt)
 UNARY_FUNCTION(scalar, scalar, sign)
 UNARY_FUNCTION(scalar, scalar, pos)
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
index fa6a2a58d0..3aaf5ddf3b 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
@@ -447,6 +447,7 @@ UNARY_FUNCTION(scalar, scalar, pow3, pow3)
 UNARY_FUNCTION(scalar, scalar, pow4, pow4)
 UNARY_FUNCTION(scalar, scalar, pow5, pow5)
 UNARY_FUNCTION(scalar, scalar, pow6, pow6)
+UNARY_FUNCTION(scalar, scalar, pow025, pow025)
 UNARY_FUNCTION(scalar, scalar, sqrt, sqrt)
 UNARY_FUNCTION(scalar, scalar, sign, sign)
 UNARY_FUNCTION(scalar, scalar, pos, pos)
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.H b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.H
index 4a2d9166b4..3b50684200 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.H
@@ -92,6 +92,7 @@ UNARY_FUNCTION(scalar, scalar, pow3, pow3)
 UNARY_FUNCTION(scalar, scalar, pow4, pow4)
 UNARY_FUNCTION(scalar, scalar, pow5, pow5)
 UNARY_FUNCTION(scalar, scalar, pow6, pow6)
+UNARY_FUNCTION(scalar, scalar, pow025, pow025)
 UNARY_FUNCTION(scalar, scalar, sqrt, sqrt)
 UNARY_FUNCTION(scalar, scalar, sign, sign)
 UNARY_FUNCTION(scalar, scalar, pos, pos)
diff --git a/src/OpenFOAM/primitives/Scalar/Scalar.H b/src/OpenFOAM/primitives/Scalar/Scalar.H
index cf5d9cb31d..69961ab2cc 100644
--- a/src/OpenFOAM/primitives/Scalar/Scalar.H
+++ b/src/OpenFOAM/primitives/Scalar/Scalar.H
@@ -90,56 +90,67 @@ inline Scalar& setComponent(Scalar& s, const direction)
     return s;
 }
 
+
 inline Scalar component(const Scalar s, const direction)
 {
     return s;
 }
 
+
 inline Scalar sign(const Scalar s)
 {
     return (s >= 0)? 1: -1;
 }
 
+
 inline Scalar pos(const Scalar s)
 {
     return (s >= 0)? 1: 0;
 }
 
+
 inline Scalar neg(const Scalar s)
 {
     return (s < 0)? 1: 0;
 }
 
+
 inline bool equal(const Scalar& s1, const Scalar& s2)
 {
     return mag(s1 - s2) <= ScalarVSMALL;
 }
 
+
 inline bool notEqual(const Scalar s1, const Scalar s2)
 {
     return mag(s1 - s2) > ScalarVSMALL;
 }
 
+
 inline Scalar limit(const Scalar s1, const Scalar s2)
 {
     return (mag(s1) < mag(s2)) ? s1: 0.0;
 }
 
+
 inline Scalar minMod(const Scalar s1, const Scalar s2)
 {
     return (mag(s1) < mag(s2)) ? s1: s2;
 }
 
+
 inline Scalar magSqr(const Scalar s)
 {
     return s*s;
 }
 
+
 inline Scalar sqr(const Scalar s)
 {
     return s*s;
 }
 
+
 inline Scalar sqrtSumSqr(const Scalar a, const Scalar b)
 {
     Scalar maga = mag(a);
@@ -155,61 +166,79 @@ inline Scalar sqrtSumSqr(const Scalar a, const Scalar b)
     }
 }
 
+
 inline Scalar pow3(const Scalar s)
 {
     return s*sqr(s);
 }
 
+
 inline Scalar pow4(const Scalar s)
 {
     return sqr(sqr(s));
 }
 
+
 inline Scalar pow5(const Scalar s)
 {
     return s*pow4(s);
 }
 
+
 inline Scalar pow6(const Scalar s)
 {
     return pow3(sqr(s));
 }
 
+
+inline Scalar pow025(const Scalar s)
+{
+    return sqrt(sqrt(s));
+}
+
+
 inline Scalar inv(const Scalar s)
 {
     return 1.0/s;
 }
 
+
 inline Scalar dot(const Scalar s1, const Scalar s2)
 {
     return s1*s2;
 }
 
+
 inline Scalar cmptMultiply(const Scalar s1, const Scalar s2)
 {
     return s1*s2;
 }
 
+
 inline Scalar cmptDivide(const Scalar s1, const Scalar s2)
 {
     return s1/s2;
 }
 
+
 inline Scalar cmptMax(const Scalar s)
 {
     return s;
 }
 
+
 inline Scalar cmptMin(const Scalar s)
 {
     return s;
 }
 
+
 inline Scalar cmptAv(const Scalar s)
 {
     return s;
 }
 
+
 inline Scalar cmptMag(const Scalar s)
 {
     return mag(s);
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 27ca388fe3..3d0274bc33 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -2661,10 +2661,22 @@ void Foam::autoLayerDriver::addLayers
     {
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-        Info<< nl
-            << "patch               faces    layers avg thickness[m]" << nl
-            << "                                    near-wall overall" << nl
-            << "-----               -----    ------ --------- -------" << endl;
+        // Find maximum length of a patch name, for a nicer output
+        label maxPatchNameLen = 0;
+        forAll(meshMover.adaptPatchIDs(), i)
+        {
+            label patchI = meshMover.adaptPatchIDs()[i];
+            word patchName = patches[patchI].name();
+            maxPatchNameLen = max(maxPatchNameLen,label(patchName.size()));
+        }
+
+        Info<< nl 
+            << setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
+            << setw(0) << " faces    layers avg thickness[m]" << nl
+            << setf(ios_base::left) << setw(maxPatchNameLen) << " "
+            << setw(0) << "                 near-wall overall" << nl
+            << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
+            << setw(0) << " -----    ------ --------- -------" << endl;
 
         forAll(meshMover.adaptPatchIDs(), i)
         {
@@ -2704,18 +2716,24 @@ void Foam::autoLayerDriver::addLayers
 
             label totNPoints = returnReduce(meshPoints.size(), sumOp<label>());
 
-            //reduce(maxThickness, maxOp<scalar>());
-            //reduce(minThickness, minOp<scalar>());
-            scalar avgThickness =
-                returnReduce(sumThickness, sumOp<scalar>())
-              / totNPoints;
-            scalar avgNearWallThickness =
-                returnReduce(sumNearWallThickness, sumOp<scalar>())
-              / totNPoints;
+            // For empty patches, totNPoints is 0.
+            scalar avgThickness = 0;
+            scalar avgNearWallThickness = 0;
 
-            Info<< setf(ios_base::left) << setw(19) << patches[patchI].name();
-            //Sout.unsetf(ios_base::left);
-            Info<< setprecision(3)
+            if (totNPoints > 0)
+            {
+                //reduce(maxThickness, maxOp<scalar>());
+                //reduce(minThickness, minOp<scalar>());
+                avgThickness =
+                    returnReduce(sumThickness, sumOp<scalar>())
+                  / totNPoints;
+                avgNearWallThickness =
+                    returnReduce(sumNearWallThickness, sumOp<scalar>())
+                  / totNPoints;
+            }
+
+            Info<< setf(ios_base::left) << setw(maxPatchNameLen)
+                << patches[patchI].name() << setprecision(3)
                 << " " << setw(8)
                 << returnReduce(patches[patchI].size(), sumOp<scalar>())
                 << " " << setw(6) << layerParams.numLayers()[patchI]
diff --git a/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C
index c38e0494ae..d5c98d52b5 100644
--- a/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -141,9 +141,18 @@ Foam::labelList Foam::decompositionMethod::decompose
     const pointField& coarsePoints
 )
 {
-    scalarField coarseWeights(0);
+    // Decompose based on agglomerated points
+    labelList coarseDistribution(decompose(coarsePoints));
 
-    return decompose(fineToCoarse, coarsePoints, coarseWeights);
+    // Rework back into decomposition for original mesh_
+    labelList fineDistribution(fineToCoarse.size());
+
+    forAll(fineDistribution, i)
+    {
+        fineDistribution[i] = coarseDistribution[fineToCoarse[i]];
+    }
+
+    return fineDistribution;
 }
 
 
diff --git a/src/lagrangian/basic/Cloud/CloudIO.C b/src/lagrangian/basic/Cloud/CloudIO.C
index c1e4b30f80..ac23db92e4 100644
--- a/src/lagrangian/basic/Cloud/CloudIO.C
+++ b/src/lagrangian/basic/Cloud/CloudIO.C
@@ -139,7 +139,13 @@ void Foam::Cloud<ParticleType>::readFields()
 
 template<class ParticleType>
 void Foam::Cloud<ParticleType>::writeFields() const
-{}
+{
+    if (this->size())
+    {
+        const ParticleType& p = *this->first();
+        ParticleType::writeFields(p.cloud());
+    }
+}
 
 
 template<class ParticleType>
diff --git a/src/lagrangian/basic/Particle/ParticleI.H b/src/lagrangian/basic/Particle/ParticleI.H
index 65ac64e525..2e6d8f868f 100644
--- a/src/lagrangian/basic/Particle/ParticleI.H
+++ b/src/lagrangian/basic/Particle/ParticleI.H
@@ -25,6 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "polyMesh.H"
+#include "wallPolyPatch.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -46,21 +47,25 @@ inline Foam::scalar Foam::Particle<ParticleType>::lambda
         Sf /= mag(Sf);
         vector Cf = mesh.faceCentres()[facei];
 
-        // move reference point for wall
+        // patch interaction
         if (!cloud_.internalFace(facei))
         {
-            const vector& C = mesh.cellCentres()[celli_];
-            scalar CCf = mag((C - Cf) & Sf);
-            // check if distance between cell centre and face centre
-            // is larger than wallImpactDistance
-            if
-            (
-                CCf
-              > static_cast<const ParticleType&>(*this).wallImpactDistance(Sf)
-            )
+            label patchi = patch(facei);
+            const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+            // move reference point for wall
+            if (isA<wallPolyPatch>(patch))
             {
-                Cf -= static_cast<const ParticleType&>(*this)
-                    .wallImpactDistance(Sf)*Sf;
+                const vector& C = mesh.cellCentres()[celli_];
+                scalar CCf = mag((C - Cf) & Sf);
+                // check if distance between cell centre and face centre
+                // is larger than wallImpactDistance
+                const ParticleType& p =
+                    static_cast<const ParticleType&>(*this);
+                if (CCf > p.wallImpactDistance(Sf))
+                {
+                    Cf -=p.wallImpactDistance(Sf)*Sf;
+                }
             }
         }
 
@@ -190,21 +195,24 @@ inline Foam::scalar Foam::Particle<ParticleType>::lambda
     Sf /= mag(Sf);
     vector Cf = mesh.faceCentres()[facei];
 
-    // move reference point for wall
+    // patch interaction
     if (!cloud_.internalFace(facei))
     {
-        const vector& C = mesh.cellCentres()[celli_];
-        scalar CCf = mag((C - Cf) & Sf);
-        // check if distance between cell centre and face centre
-        // is larger than wallImpactDistance
-        if
-        (
-            CCf
-          > static_cast<const ParticleType&>(*this).wallImpactDistance(Sf)
-        )
+        label patchi = patch(facei);
+        const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+        // move reference point for wall
+        if (isA<wallPolyPatch>(patch))
         {
-            Cf -= static_cast<const ParticleType&>(*this)
-                .wallImpactDistance(Sf)*Sf;
+            const vector& C = mesh.cellCentres()[celli_];
+            scalar CCf = mag((C - Cf) & Sf);
+            // check if distance between cell centre and face centre
+            // is larger than wallImpactDistance
+            const ParticleType& p = static_cast<const ParticleType&>(*this);
+            if (CCf > p.wallImpactDistance(Sf))
+            {
+                Cf -=p.wallImpactDistance(Sf)*Sf;
+            }
         }
     }
 
diff --git a/src/lagrangian/basic/indexedParticle/indexedParticle.H b/src/lagrangian/basic/indexedParticle/indexedParticle.H
index 7023f0ba6d..be2ecd65d9 100644
--- a/src/lagrangian/basic/indexedParticle/indexedParticle.H
+++ b/src/lagrangian/basic/indexedParticle/indexedParticle.H
@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class indexedParticle Declaration
+                      Class indexedParticle Declaration
 \*---------------------------------------------------------------------------*/
 
 class indexedParticle
diff --git a/src/lagrangian/basic/indexedParticle/indexedParticleCloud.C b/src/lagrangian/basic/indexedParticle/indexedParticleCloud.C
index a9ebc140ec..81ce7e3322 100644
--- a/src/lagrangian/basic/indexedParticle/indexedParticleCloud.C
+++ b/src/lagrangian/basic/indexedParticle/indexedParticleCloud.C
@@ -45,20 +45,17 @@ defineTemplateTypeNameAndDebug(Cloud<indexedParticle>, 0);
 Foam::indexedParticleCloud::indexedParticleCloud
 (
     const polyMesh& mesh,
-    const word& cloudName
+    const word& cloudName,
+    bool readFields
 )
 :
     Cloud<indexedParticle>(mesh, cloudName, false)
 {
-    indexedParticle::readFields(*this);
+    if (readFields)
+    {
+        indexedParticle::readFields(*this);
+    }
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::indexedParticleCloud::writeFields() const
-{
-    indexedParticle::writeFields(*this);
-}
-
 // ************************************************************************* //
diff --git a/src/lagrangian/basic/indexedParticle/indexedParticleCloud.H b/src/lagrangian/basic/indexedParticle/indexedParticleCloud.H
index 65033e0a9d..0f8f39edbb 100644
--- a/src/lagrangian/basic/indexedParticle/indexedParticleCloud.H
+++ b/src/lagrangian/basic/indexedParticle/indexedParticleCloud.H
@@ -45,7 +45,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class indexedParticleCloud Declaration
+                    Class indexedParticleCloud Declaration
 \*---------------------------------------------------------------------------*/
 
 class indexedParticleCloud
@@ -69,14 +69,9 @@ public:
         indexedParticleCloud
         (
             const polyMesh&,
-            const word& cloudName = "defaultCloud"
+            const word& cloudName = "defaultCloud",
+            bool readFields = true
         );
-
-
-    // Member Functions
-
-        //- Write fields
-        virtual void writeFields() const;
 };
 
 
diff --git a/src/lagrangian/basic/passiveParticle/passiveParticleCloud.C b/src/lagrangian/basic/passiveParticle/passiveParticleCloud.C
index ebb67d26de..79f937aaa0 100644
--- a/src/lagrangian/basic/passiveParticle/passiveParticleCloud.C
+++ b/src/lagrangian/basic/passiveParticle/passiveParticleCloud.C
@@ -43,12 +43,16 @@ defineTemplateTypeNameAndDebug(Cloud<passiveParticle>, 0);
 Foam::passiveParticleCloud::passiveParticleCloud
 (
     const polyMesh& mesh,
-    const word& cloudName
+    const word& cloudName,
+    bool readFields
 )
 :
     Cloud<passiveParticle>(mesh, cloudName, false)
 {
-    readFields();
+    if (readFields)
+    {
+        passiveParticle::readFields(*this);
+    }
 }
 
 
@@ -63,18 +67,4 @@ Foam::passiveParticleCloud::passiveParticleCloud
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::passiveParticleCloud::readFields()
-{
-    passiveParticle::readFields(*this);
-}
-
-
-void Foam::passiveParticleCloud::writeFields() const
-{
-    passiveParticle::writeFields(*this);
-}
-
-
 // ************************************************************************* //
diff --git a/src/lagrangian/basic/passiveParticle/passiveParticleCloud.H b/src/lagrangian/basic/passiveParticle/passiveParticleCloud.H
index e3c0266fa3..1ef383867b 100644
--- a/src/lagrangian/basic/passiveParticle/passiveParticleCloud.H
+++ b/src/lagrangian/basic/passiveParticle/passiveParticleCloud.H
@@ -45,7 +45,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class passiveParticleCloud Declaration
+                   Class passiveParticleCloud Declaration
 \*---------------------------------------------------------------------------*/
 
 class passiveParticleCloud
@@ -72,7 +72,8 @@ public:
         passiveParticleCloud
         (
             const polyMesh&,
-            const word& cloudName = "defaultCloud"
+            const word& cloudName = "defaultCloud",
+            bool readFields = true
         );
 
         //- Construct from mesh, cloud name, and a list of particles
@@ -82,15 +83,6 @@ public:
             const word& cloudName,
             const IDLList<passiveParticle>& particles
         );
-
-
-    // Member Functions
-
-        //- Read fields
-        virtual void readFields();
-
-        //- Write fields
-        virtual void writeFields() const;
 };
 
 
diff --git a/src/lagrangian/coalCombustion/CoalCloud/CoalCloud.C b/src/lagrangian/coalCombustion/CoalCloud/CoalCloud.C
deleted file mode 100644
index 0f274d70bf..0000000000
--- a/src/lagrangian/coalCombustion/CoalCloud/CoalCloud.C
+++ /dev/null
@@ -1,70 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "CoalCloud.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-Foam::CoalCloud<ThermoType>::CoalCloud
-(
-    const word& cloudName,
-    const volScalarField& rho,
-    const volVectorField& U,
-    const dimensionedVector& g,
-    basicThermo& thermo
-)
-:
-    ReactingMultiphaseCloud<CoalParcel<ThermoType> >
-    (
-        cloudName,
-        rho,
-        U,
-        g,
-        thermo
-    )
-{
-    CoalParcel<ThermoType>::readFields(*this);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-Foam::CoalCloud<ThermoType>::~CoalCloud()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class ThermoType>
-void Foam::CoalCloud<ThermoType>::writeFields() const
-{
-    CoalParcel<ThermoType>::writeFields(*this);
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/coalCombustion/CoalCloud/CoalCloud.H b/src/lagrangian/coalCombustion/CoalCloud/CoalCloud.H
index 8516508d8b..0ca7aa497c 100644
--- a/src/lagrangian/coalCombustion/CoalCloud/CoalCloud.H
+++ b/src/lagrangian/coalCombustion/CoalCloud/CoalCloud.H
@@ -26,9 +26,7 @@ Class
     CoalCloud
 
 Description
-
-SourceFiles
-    CoalCloud.C
+    Coal cloud templated on the type of carrier phase thermodynamics
 
 \*---------------------------------------------------------------------------*/
 
@@ -37,69 +35,21 @@ SourceFiles
 
 #include "ReactingMultiphaseCloud.H"
 #include "CoalParcel.H"
+#include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+    typedef ReactingMultiphaseCloud<CoalParcel<constGasThermoPhysics> >
+        constThermoCoalCloud;
 
-/*---------------------------------------------------------------------------*\
-                          Class CoalCloud Declaration
-\*---------------------------------------------------------------------------*/
+    typedef ReactingMultiphaseCloud<CoalParcel<gasThermoPhysics> >
+        thermoCoalCloud;
 
-template<class ThermoType>
-class CoalCloud
-:
-    public ReactingMultiphaseCloud<CoalParcel<ThermoType> >
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        CoalCloud(const CoalCloud&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const CoalCloud&);
-
-
-public:
-
-    //-Runtime type information
-    TypeName("CoalCloud");
-
-
-    // Constructors
-
-        //- Construct given carrier gas fields
-        CoalCloud
-        (
-            const word& cloudName,
-            const volScalarField& rho,
-            const volVectorField& U,
-            const dimensionedVector& g,
-            basicThermo& thermo
-        );
-
-
-    //- Destructor
-    ~CoalCloud();
-
-
-    // Member Functions
-
-        //- Write fields
-        virtual void writeFields() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "CoalCloud.C"
-#endif
+    typedef ReactingMultiphaseCloud<CoalParcel<icoPoly8ThermoPhysics> >
+        icoPoly8ThermoCoalCloud;
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/coalCombustion/CoalCloud/defineCoalCloud.C b/src/lagrangian/coalCombustion/CoalCloud/defineCoalCloud.C
deleted file mode 100644
index cc0fa34f49..0000000000
--- a/src/lagrangian/coalCombustion/CoalCloud/defineCoalCloud.C
+++ /dev/null
@@ -1,38 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "createReactingCloudTypes.H"
-#include "CoalCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    createReactingCloudType(CoalCloud);
-};
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/coalCombustion/Make/files b/src/lagrangian/coalCombustion/Make/files
index 2c7aac61ab..e00ab17131 100644
--- a/src/lagrangian/coalCombustion/Make/files
+++ b/src/lagrangian/coalCombustion/Make/files
@@ -1,6 +1,3 @@
-/* Coal cloud */
-CoalCloud/defineCoalCloud.C
-
 /* Coal parcel and sub-models */
 CoalParcel/defineCoalParcel.C
 CoalParcel/makeCoalParcelSubmodels.C
diff --git a/src/lagrangian/coalCombustion/include/createCoalParcelTypes.H b/src/lagrangian/coalCombustion/include/createCoalParcelTypes.H
index 7d7f95da36..f81595ad59 100644
--- a/src/lagrangian/coalCombustion/include/createCoalParcelTypes.H
+++ b/src/lagrangian/coalCombustion/include/createCoalParcelTypes.H
@@ -92,8 +92,7 @@ License
     (                                                                         \
         ReactingMultiphaseCloud<ParcelType##ThermoType>,                      \
         0                                                                     \
-    );                                                                        \
-    defineParcelTypeNameAndDebug(CoalCloud<ParcelType##ThermoType>, 0);
+    );
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/dieselSpray/parcel/parcelIO.C b/src/lagrangian/dieselSpray/parcel/parcelIO.C
index 4ae7bbc90a..f9c0293383 100644
--- a/src/lagrangian/dieselSpray/parcel/parcelIO.C
+++ b/src/lagrangian/dieselSpray/parcel/parcelIO.C
@@ -96,10 +96,7 @@ Foam::parcel::parcel
 }
 
 
-void Foam::parcel::readFields
-(
-    Cloud<parcel>& c
-)
+void Foam::parcel::readFields(Cloud<parcel>& c)
 {
     if (!c.size())
     {
@@ -195,10 +192,7 @@ void Foam::parcel::readFields
 }
 
 
-void Foam::parcel::writeFields
-(
-    const Cloud<parcel>& c
-)
+void Foam::parcel::writeFields(const Cloud<parcel>& c)
 {
     Particle<parcel>::writeFields(c);
 
diff --git a/src/lagrangian/dieselSpray/spray/spray.C b/src/lagrangian/dieselSpray/spray/spray.C
index 1f1945ef5d..8dc5ad5fe0 100644
--- a/src/lagrangian/dieselSpray/spray/spray.C
+++ b/src/lagrangian/dieselSpray/spray/spray.C
@@ -59,7 +59,8 @@ Foam::spray::spray
     const basicMultiComponentMixture& composition,
     const PtrList<gasThermoPhysics>& gasProperties,
     const dictionary&,
-    const dimensionedVector& g
+    const dimensionedVector& g,
+    bool readFields
 )
 :
     Cloud<parcel>(U.mesh(), false), // suppress className checking on positions
@@ -345,7 +346,10 @@ Foam::spray::spray
         }
     }
 
-    parcel::readFields(*this);
+    if (readFields)
+    {
+        parcel::readFields(*this);
+    }
 }
 
 
@@ -355,12 +359,4 @@ Foam::spray::~spray()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::spray::writeFields() const
-{
-    parcel::writeFields(*this);
-}
-
-
 // ************************************************************************* //
diff --git a/src/lagrangian/dieselSpray/spray/spray.H b/src/lagrangian/dieselSpray/spray/spray.H
index 60127656e7..5c641ec84e 100644
--- a/src/lagrangian/dieselSpray/spray/spray.H
+++ b/src/lagrangian/dieselSpray/spray/spray.H
@@ -194,7 +194,8 @@ public:
             const basicMultiComponentMixture& composition,
             const PtrList<gasThermoPhysics>& gasProperties,
             const dictionary& thermophysicalProperties,
-            const dimensionedVector& g
+            const dimensionedVector& g,
+            bool readFields = true
         );
 
 
@@ -339,12 +340,6 @@ public:
 
             //- Return Ambient Temperature
             void calculateAmbientTemperature();
-
-
-    // I/O
-
-            //- Write fields
-            virtual void writeFields() const;
 };
 
 
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
index 23aa3dc478..5ac29101e5 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
@@ -28,11 +28,9 @@ License
 #include "BinaryCollisionModel.H"
 #include "WallInteractionModel.H"
 #include "InflowBoundaryModel.H"
+#include "constants.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-template<class ParcelType>
-Foam::scalar Foam::DsmcCloud<ParcelType>::kb = 1.380650277e-23;
+using namespace Foam::constant;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -526,7 +524,8 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
 (
     const word& cloudName,
     const volScalarField& T,
-    const volVectorField& U
+    const volVectorField& U,
+    bool readFields
 )
 :
     Cloud<ParcelType>(T.mesh(), cloudName, false),
@@ -630,6 +629,11 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
     {
         collisionSelectionRemainder_[i] = rndGen_.scalar01();
     }
+
+    if (readFields)
+    {
+        ParcelType::readFields(*this);
+    }
 }
 
 
@@ -853,7 +857,7 @@ Foam::vector Foam::DsmcCloud<ParcelType>::equipartitionLinearVelocity
 )
 {
     return
-        sqrt(kb*temperature/mass)
+        sqrt(physicoChemical::k.value()*temperature/mass)
        *vector
         (
             rndGen_.GaussNormal(),
@@ -879,7 +883,7 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
     else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
     {
         // Special case for iDof = 2, i.e. diatomics;
-        Ei = -log(rndGen_.scalar01())*kb*temperature;
+        Ei = -log(rndGen_.scalar01())*physicoChemical::k.value()*temperature;
     }
     else
     {
@@ -897,7 +901,7 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
 
         } while (P < rndGen_.scalar01());
 
-        Ei = energyRatio*kb*temperature;
+        Ei = energyRatio*physicoChemical::k.value()*temperature;
     }
 
     return Ei;
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
index 76b1dd7e79..c3765ba0df 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
@@ -171,12 +171,6 @@ class DsmcCloud
 
 public:
 
-    // Static data members
-
-        //- Boltzmann constant
-        static scalar kb;
-
-
     // Constructors
 
         //- Construct given name and mesh, will read Parcels from file
@@ -184,7 +178,8 @@ public:
         (
             const word& cloudName,
             const volScalarField& T,
-            const volVectorField& U
+            const volVectorField& U,
+            bool readFields = true
         );
 
         //- Construct given name and mesh.  Used to initialise.
@@ -338,6 +333,7 @@ public:
                     scalar mass
                 ) const;
 
+
             // Sub-models
 
                 //- Return reference to binary elastic collision model
@@ -429,11 +425,8 @@ public:
             //- Evolve the cloud (move, collide)
             void evolve();
 
-
             //- Clear the Cloud
             inline void clear();
-
-
 };
 
 
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
index 492bc98e18..1c58b87c24 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
@@ -24,7 +24,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "mathConstants.H"
+#include "constants.H"
+
+using namespace Foam::constant;
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -314,7 +316,8 @@ inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianAverageSpeed
     scalar mass
 ) const
 {
-    return 2.0*sqrt(2.0*kb*temperature/(constant::math::pi*mass));
+    return
+        2.0*sqrt(2.0*physicoChemical::k.value()*temperature/(math::pi*mass));
 }
 
 
@@ -325,7 +328,8 @@ inline Foam::scalarField Foam::DsmcCloud<ParcelType>::maxwellianAverageSpeed
     scalar mass
 ) const
 {
-    return 2.0*sqrt(2.0*kb*temperature/(constant::math::pi*mass));
+    return
+        2.0*sqrt(2.0*physicoChemical::k.value()*temperature/(math::pi*mass));
 }
 
 
@@ -336,7 +340,7 @@ inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianRMSSpeed
     scalar mass
 ) const
 {
-    return sqrt(3.0*kb*temperature/mass);
+    return sqrt(3.0*physicoChemical::k.value()*temperature/mass);
 }
 
 
@@ -347,7 +351,7 @@ inline Foam::scalarField Foam::DsmcCloud<ParcelType>::maxwellianRMSSpeed
     scalar mass
 ) const
 {
-    return sqrt(3.0*kb*temperature/mass);
+    return sqrt(3.0*physicoChemical::k.value()*temperature/mass);
 }
 
 
@@ -359,7 +363,7 @@ Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
     scalar mass
 ) const
 {
-    return sqrt(2.0*kb*temperature/mass);
+    return sqrt(2.0*physicoChemical::k.value()*temperature/mass);
 }
 
 
@@ -371,7 +375,7 @@ Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
     scalar mass
 ) const
 {
-    return sqrt(2.0*kb*temperature/mass);
+    return sqrt(2.0*physicoChemical::k.value()*temperature/mass);
 }
 
 
diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C
index 425e47a705..7f94f3e739 100644
--- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C
@@ -44,9 +44,7 @@ Foam::dsmcCloud::dsmcCloud
 )
 :
     DsmcCloud<dsmcParcel>(cloudName, T, U)
-{
-    dsmcParcel::readFields(*this);
-}
+{}
 
 
 Foam::dsmcCloud::dsmcCloud
@@ -65,11 +63,4 @@ Foam::dsmcCloud::~dsmcCloud()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::dsmcCloud::writeFields() const
-{
-    dsmcParcel::writeFields(*this);
-}
-
 // ************************************************************************* //
diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H
index ad63aca9fc..99343c2b3e 100644
--- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H
+++ b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H
@@ -86,12 +86,6 @@ public:
 
     //- Destructor
     ~dsmcCloud();
-
-
-    // Member functions
-
-        //- Write fields
-        virtual void writeFields() const;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H
index d96d515843..918767b1af 100644
--- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H
@@ -308,9 +308,9 @@ public:
 
         // I-O
 
-            static void readFields(DsmcCloud<ParcelType>& c);
+            static void readFields(Cloud<ParcelType>& c);
 
-            static void writeFields(const DsmcCloud<ParcelType>& c);
+            static void writeFields(const Cloud<ParcelType>& c);
 
 
     // Ostream Operator
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C
index 6677591009..6e43e47508 100644
--- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C
@@ -74,10 +74,7 @@ Foam::DsmcParcel<ParcelType>::DsmcParcel
 
 
 template<class ParcelType>
-void Foam::DsmcParcel<ParcelType>::readFields
-(
-    DsmcCloud<ParcelType>& c
-)
+void Foam::DsmcParcel<ParcelType>::readFields(Cloud<ParcelType>& c)
 {
     if (!c.size())
     {
@@ -107,10 +104,7 @@ void Foam::DsmcParcel<ParcelType>::readFields
 
 
 template<class ParcelType>
-void Foam::DsmcParcel<ParcelType>::writeFields
-(
-    const DsmcCloud<ParcelType>& c
-)
+void Foam::DsmcParcel<ParcelType>::writeFields(const Cloud<ParcelType>& c)
 {
     Particle<ParcelType>::writeFields(c);
 
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
index 3ab8c5511e..e445f7cde5 100644
--- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
@@ -25,7 +25,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "LarsenBorgnakkeVariableHardSphere.H"
-#include "mathConstants.H"
+#include "constants.H"
+
+using namespace Foam::constant;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -157,8 +159,8 @@ Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::sigmaTcR
 
     // calculating cross section = pi*dPQ^2, where dPQ is from Bird, eq. 4.79
     scalar sigmaTPQ =
-        constant::math::pi*dPQ*dPQ
-       *pow(2.0*CloudType::kb*Tref_/(mR*cR*cR), omegaPQ - 0.5)
+        math::pi*dPQ*dPQ
+       *pow(2.0*physicoChemical::k.value()*Tref_/(mR*cR*cR), omegaPQ - 0.5)
        /exp(Foam::lgamma(2.5 - omegaPQ));
 
     return sigmaTPQ*cR;
@@ -253,7 +255,7 @@ void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
 
     scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
 
-    scalar phi = constant::math::twoPi*rndGen.scalar01();
+    scalar phi = math::twoPi*rndGen.scalar01();
 
     vector postCollisionRelU =
         cR
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
index 8afc892c84..3174445ac6 100644
--- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
@@ -25,7 +25,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "VariableHardSphere.H"
-#include "mathConstants.H"
+#include "constants.H"
+
+using namespace Foam::constant;
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -91,8 +93,8 @@ Foam::scalar Foam::VariableHardSphere<CloudType>::sigmaTcR
 
     // calculating cross section = pi*dPQ^2, where dPQ is from Bird, eq. 4.79
     scalar sigmaTPQ =
-        constant::math::pi*dPQ*dPQ
-       *pow(2.0*CloudType::kb*Tref_/(mR*cR*cR), omegaPQ - 0.5)
+        math::pi*dPQ*dPQ
+       *pow(2.0*physicoChemical::k.value()*Tref_/(mR*cR*cR), omegaPQ - 0.5)
        /exp(Foam::lgamma(2.5 - omegaPQ));
 
     return sigmaTPQ*cR;
@@ -126,7 +128,7 @@ void Foam::VariableHardSphere<CloudType>::collide
 
     scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
 
-    scalar phi = constant::math::twoPi*rndGen.scalar01();
+    scalar phi = math::twoPi*rndGen.scalar01();
 
     vector postCollisionRelU =
         cR
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index 5472e6c237..dded033d73 100644
--- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -25,7 +25,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "FreeStream.H"
-#include "mathConstants.H"
+#include "constants.H"
+
+using namespace Foam::constant;
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -131,7 +133,7 @@ void Foam::FreeStream<CloudType>::inflow()
 
     Random& rndGen(cloud.rndGen());
 
-    scalar sqrtPi = sqrt(constant::math::pi);
+    scalar sqrtPi = sqrt(math::pi);
 
     label particlesInserted = 0;
 
@@ -370,7 +372,7 @@ void Foam::FreeStream<CloudType>::inflow()
                     } while (P < rndGen.scalar01());
 
                     vector U =
-                        sqrt(CloudType::kb*faceTemperature/mass)
+                        sqrt(physicoChemical::k.value()*faceTemperature/mass)
                        *(
                             rndGen.GaussNormal()*t1
                           + rndGen.GaussNormal()*t2
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
index e7213561c0..569e9b21e6 100644
--- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
@@ -25,6 +25,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "MaxwellianThermal.H"
+#include "constants.H"
+
+using namespace Foam::constant;
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -108,7 +111,7 @@ void Foam::MaxwellianThermal<CloudType>::correct
     scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
 
     U =
-        sqrt(CloudType::kb*T/mass)
+        sqrt(physicoChemical::k.value()*T/mass)
        *(
             rndGen.GaussNormal()*tw1
           + rndGen.GaussNormal()*tw2
diff --git a/src/lagrangian/intermediate/Make/files b/src/lagrangian/intermediate/Make/files
index 7911383146..270bc40e2e 100644
--- a/src/lagrangian/intermediate/Make/files
+++ b/src/lagrangian/intermediate/Make/files
@@ -18,13 +18,6 @@ $(BASECLOUDS)/reactingCloud/reactingCloud.C
 $(BASECLOUDS)/reactingMultiphaseCloud/reactingMultiphaseCloud.C
 
 
-/* Cloud container/injection mechanisms */
-$(DERIVEDCLOUDS)/basicKinematicCloud/basicKinematicCloud.C
-$(DERIVEDCLOUDS)/basicThermoCloud/basicThermoCloud.C
-$(DERIVEDCLOUDS)/BasicReactingCloud/defineBasicReactingCloud.C
-$(DERIVEDCLOUDS)/BasicReactingMultiphaseCloud/defineBasicReactingMultiphaseCloud.C
-
-
 /* kinematic parcel sub-models */
 KINEMATICPARCEL=$(DERIVEDPARCELS)/basicKinematicParcel
 $(KINEMATICPARCEL)/basicKinematicParcel.C
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 650cb7ddcb..f034fbbaa3 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -43,7 +43,8 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     const volScalarField& rho,
     const volVectorField& U,
     const volScalarField& mu,
-    const dimensionedVector& g
+    const dimensionedVector& g,
+    bool readFields
 )
 :
     Cloud<ParcelType>(rho.mesh(), cloudName, false),
@@ -136,7 +137,12 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
         mesh_,
         dimensionedVector("zero", dimMass*dimVelocity, vector::zero)
     )
-{}
+{
+    if (readFields)
+    {
+        ParcelType::readFields(*this);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index ee2b3bafe2..4d5295b08d 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -204,7 +204,8 @@ public:
             const volScalarField& rho,
             const volVectorField& U,
             const volScalarField& mu,
-            const dimensionedVector& g
+            const dimensionedVector& g,
+            bool readFields = true
         );
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index 1e58d2876c..6f4b4bc763 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -67,10 +67,11 @@ Foam::ReactingCloud<ParcelType>::ReactingCloud
     const volScalarField& rho,
     const volVectorField& U,
     const dimensionedVector& g,
-    basicThermo& thermo
+    basicThermo& thermo,
+    bool readFields
 )
 :
-    ThermoCloud<ParcelType>(cloudName, rho, U, g, thermo),
+    ThermoCloud<ParcelType>(cloudName, rho, U, g, thermo, false),
     reactingCloud(),
     constProps_(this->particleProperties()),
     mcCarrierThermo_
@@ -118,6 +119,11 @@ Foam::ReactingCloud<ParcelType>::ReactingCloud
             )
         );
     }
+
+    if (readFields)
+    {
+        ParcelType::readFields(*this);
+    }
 }
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
index 74b75ea416..23f579649a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
@@ -152,7 +152,8 @@ public:
             const volScalarField& rho,
             const volVectorField& U,
             const dimensionedVector& g,
-            basicThermo& thermo
+            basicThermo& thermo,
+            bool readFields = true
         );
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
index 5e7dee64bc..e81aa2ee8d 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
@@ -38,10 +38,11 @@ Foam::ReactingMultiphaseCloud<ParcelType>::ReactingMultiphaseCloud
     const volScalarField& rho,
     const volVectorField& U,
     const dimensionedVector& g,
-    basicThermo& thermo
+    basicThermo& thermo,
+    bool readFields
 )
 :
-    ReactingCloud<ParcelType>(cloudName, rho, U, g, thermo),
+    ReactingCloud<ParcelType>(cloudName, rho, U, g, thermo, false),
     reactingMultiphaseCloud(),
     constProps_(this->particleProperties()),
     devolatilisationModel_
@@ -61,7 +62,12 @@ Foam::ReactingMultiphaseCloud<ParcelType>::ReactingMultiphaseCloud
         )
     ),
     dMassDevolatilisation_(0.0)
-{}
+{
+    if (readFields)
+    {
+        ParcelType::readFields(*this);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.H b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.H
index b073308407..61dcdf6679 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.H
@@ -134,7 +134,8 @@ public:
             const volScalarField& rho,
             const volVectorField& U,
             const dimensionedVector& g,
-            basicThermo& thermo
+            basicThermo& thermo,
+            bool readFields = true
         );
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index 0585629246..48ef62ee19 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -39,7 +39,8 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
     const volScalarField& rho,
     const volVectorField& U,
     const dimensionedVector& g,
-    basicThermo& thermo
+    basicThermo& thermo,
+    bool readFields
 )
 :
     KinematicCloud<ParcelType>
@@ -48,7 +49,8 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
         rho,
         U,
         thermo.mu(),
-        g
+        g,
+        false
     ),
     thermoCloud(),
     constProps_(this->particleProperties()),
@@ -98,7 +100,12 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
         this->mesh(),
         dimensionedScalar("zero", dimensionSet(1, 2, -2, 0, 0), 0.0)
     )
-{}
+{
+    if (readFields)
+    {
+        ParcelType::readFields(*this);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
index 151b076411..793a563f5b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
@@ -140,7 +140,8 @@ public:
             const volScalarField& rho,
             const volVectorField& U,
             const dimensionedVector& g,
-            basicThermo& thermo
+            basicThermo& thermo,
+            bool readFields = true
         );
 
 
diff --git a/src/lagrangian/intermediate/clouds/derived/BasicReactingCloud/BasicReactingCloud.C b/src/lagrangian/intermediate/clouds/derived/BasicReactingCloud/BasicReactingCloud.C
deleted file mode 100644
index e59ca18ecf..0000000000
--- a/src/lagrangian/intermediate/clouds/derived/BasicReactingCloud/BasicReactingCloud.C
+++ /dev/null
@@ -1,70 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "BasicReactingCloud.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-Foam::BasicReactingCloud<ThermoType>::BasicReactingCloud
-(
-    const word& cloudName,
-    const volScalarField& rho,
-    const volVectorField& U,
-    const dimensionedVector& g,
-    basicThermo& thermo
-)
-:
-    ReactingCloud<BasicReactingParcel<ThermoType> >
-    (
-        cloudName,
-        rho,
-        U,
-        g,
-        thermo
-    )
-{
-    BasicReactingParcel<ThermoType>::readFields(*this);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-Foam::BasicReactingCloud<ThermoType>::~BasicReactingCloud()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class ThermoType>
-void Foam::BasicReactingCloud<ThermoType>::writeFields() const
-{
-    BasicReactingParcel<ThermoType>::writeFields(*this);
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/derived/BasicReactingCloud/BasicReactingCloud.H b/src/lagrangian/intermediate/clouds/derived/BasicReactingCloud/BasicReactingCloud.H
index 0dea7344a6..7c737d1c62 100644
--- a/src/lagrangian/intermediate/clouds/derived/BasicReactingCloud/BasicReactingCloud.H
+++ b/src/lagrangian/intermediate/clouds/derived/BasicReactingCloud/BasicReactingCloud.H
@@ -38,73 +38,21 @@ SourceFiles
 
 #include "ReactingCloud.H"
 #include "BasicReactingParcel.H"
+#include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+    typedef ReactingCloud<BasicReactingParcel<constGasThermoPhysics> >
+        constThermoReactingCloud;
 
-// Forward declaration of classes
-template<class ThermoType>
-class BasicReactingCloud;
+    typedef ReactingCloud<BasicReactingParcel<gasThermoPhysics> >
+        thermoReactingCloud;
 
-/*---------------------------------------------------------------------------*\
-                        Class BasicReactingCloud Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class ThermoType>
-class BasicReactingCloud
-:
-    public ReactingCloud<BasicReactingParcel<ThermoType> >
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        BasicReactingCloud(const BasicReactingCloud&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const BasicReactingCloud&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("BasicReactingCloud");
-
-
-    // Constructors
-
-        //- Construct given carrier gas fields
-        BasicReactingCloud
-        (
-            const word& cloudName,
-            const volScalarField& rho,
-            const volVectorField& U,
-            const dimensionedVector& g,
-            basicThermo& thermo
-        );
-
-
-    //- Destructor
-    ~BasicReactingCloud();
-
-
-    // Member Functions
-
-        //- Write fields
-        virtual void writeFields() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "BasicReactingCloud.C"
-#endif
+    typedef ReactingCloud<BasicReactingParcel<icoPoly8ThermoPhysics> >
+        icoPoly8ThermoReactingCloud;
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/intermediate/clouds/derived/BasicReactingCloud/defineBasicReactingCloud.C b/src/lagrangian/intermediate/clouds/derived/BasicReactingCloud/defineBasicReactingCloud.C
deleted file mode 100644
index ca5b9e8090..0000000000
--- a/src/lagrangian/intermediate/clouds/derived/BasicReactingCloud/defineBasicReactingCloud.C
+++ /dev/null
@@ -1,38 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "createReactingCloudTypes.H"
-#include "BasicReactingCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    createReactingCloudType(BasicReactingCloud);
-};
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.C b/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.C
deleted file mode 100644
index 6f08db2016..0000000000
--- a/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.C
+++ /dev/null
@@ -1,70 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "BasicReactingMultiphaseCloud.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-Foam::BasicReactingMultiphaseCloud<ThermoType>::BasicReactingMultiphaseCloud
-(
-    const word& cloudName,
-    const volScalarField& rho,
-    const volVectorField& U,
-    const dimensionedVector& g,
-    basicThermo& thermo
-)
-:
-    ReactingMultiphaseCloud<BasicReactingMultiphaseParcel<ThermoType> >
-    (
-        cloudName,
-        rho,
-        U,
-        g,
-        thermo
-    )
-{
-    BasicReactingMultiphaseParcel<ThermoType>::readFields(*this);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-Foam::BasicReactingMultiphaseCloud<ThermoType>::~BasicReactingMultiphaseCloud()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class ThermoType>
-void Foam::BasicReactingMultiphaseCloud<ThermoType>::writeFields() const
-{
-    BasicReactingMultiphaseParcel<ThermoType>::writeFields(*this);
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.H b/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.H
index 11b6fb2378..f02475a316 100644
--- a/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.H
+++ b/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.H
@@ -26,10 +26,8 @@ Class
     Foam::BasicReactingMultiphaseCloud
 
 Description
-    Reacting multiphase cloud templated on the type of carrier phase thermodynamics
-
-SourceFiles
-    BasicReactingMultiphaseCloud.C
+    Reacting multiphase cloud templated on the type of carrier phase
+    thermodynamics
 
 \*---------------------------------------------------------------------------*/
 
@@ -38,73 +36,21 @@ SourceFiles
 
 #include "ReactingMultiphaseCloud.H"
 #include "BasicReactingMultiphaseParcel.H"
+#include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+    typedef ReactingMultiphaseCloud<BasicReactingParcel<constGasThermoPhysics> >
+        constThermoReactingMultiphaseCloud;
 
-// Forward declaration of classes
-template<class ThermoType>
-class BasicReactingMultiphaseCloud;
+    typedef ReactingMultiphaseCloud<BasicReactingParcel<gasThermoPhysics> >
+        thermoReactingMultiphaseCloud;
 
-/*---------------------------------------------------------------------------*\
-                Class BasicReactingMultiphaseCloud Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class ThermoType>
-class BasicReactingMultiphaseCloud
-:
-    public ReactingMultiphaseCloud<BasicReactingMultiphaseParcel<ThermoType> >
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        BasicReactingMultiphaseCloud(const BasicReactingMultiphaseCloud&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const BasicReactingMultiphaseCloud&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("BasicReactingMultiphaseCloud");
-
-
-    // Constructors
-
-        //- Construct given carrier gas fields
-        BasicReactingMultiphaseCloud
-        (
-            const word& cloudName,
-            const volScalarField& rho,
-            const volVectorField& U,
-            const dimensionedVector& g,
-            basicThermo& thermo
-        );
-
-
-    //- Destructor
-    ~BasicReactingMultiphaseCloud();
-
-
-    // Member Functions
-
-        //- Write fields
-        virtual void writeFields() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "BasicReactingMultiphaseCloud.C"
-#endif
+    typedef ReactingMultiphaseCloud<BasicReactingParcel<icoPoly8ThermoPhysics> >
+        icoPoly8ThermoReactingMultiphaseCloud;
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/defineBasicReactingMultiphaseCloud.C b/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/defineBasicReactingMultiphaseCloud.C
deleted file mode 100644
index d23115997d..0000000000
--- a/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/defineBasicReactingMultiphaseCloud.C
+++ /dev/null
@@ -1,38 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "createReactingCloudTypes.H"
-#include "BasicReactingMultiphaseCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    createReactingCloudType(BasicReactingMultiphaseCloud);
-};
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/derived/basicKinematicCloud/basicKinematicCloud.C b/src/lagrangian/intermediate/clouds/derived/basicKinematicCloud/basicKinematicCloud.C
deleted file mode 100644
index 15ab7afdfc..0000000000
--- a/src/lagrangian/intermediate/clouds/derived/basicKinematicCloud/basicKinematicCloud.C
+++ /dev/null
@@ -1,68 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "basicKinematicCloud.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(basicKinematicCloud, 0);
-};
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::basicKinematicCloud::basicKinematicCloud
-(
-    const word& cloudName,
-    const volScalarField& rho,
-    const volVectorField& U,
-    const volScalarField& mu,
-    const dimensionedVector& g
-)
-:
-    KinematicCloud<basicKinematicParcel>(cloudName, rho, U, mu, g)
-{
-    basicKinematicParcel::readFields(*this);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::basicKinematicCloud::~basicKinematicCloud()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::basicKinematicCloud::writeFields() const
-{
-    basicKinematicParcel::writeFields(*this);
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/derived/basicKinematicCloud/basicKinematicCloud.H b/src/lagrangian/intermediate/clouds/derived/basicKinematicCloud/basicKinematicCloud.H
index 9b2c6003d3..3a277e415a 100644
--- a/src/lagrangian/intermediate/clouds/derived/basicKinematicCloud/basicKinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/derived/basicKinematicCloud/basicKinematicCloud.H
@@ -43,56 +43,8 @@ SourceFiles
 
 namespace Foam
 {
-
-/*---------------------------------------------------------------------------*\
-                      Class basicKinematicCloud Declaration
-\*---------------------------------------------------------------------------*/
-
-class basicKinematicCloud
-:
-    public KinematicCloud<basicKinematicParcel>
-{
-    // Private member functions
-
-        //- Disallow default bitwise copy construct
-        basicKinematicCloud(const basicKinematicCloud&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const basicKinematicCloud&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("basicKinematicCloud");
-
-
-    // Constructors
-
-        //- Construct from components
-        basicKinematicCloud
-        (
-            const word& cloudName,
-            const volScalarField& rho,
-            const volVectorField& U,
-            const volScalarField& mu,
-            const dimensionedVector& g
-        );
-
-
-    //- Destructor
-    ~basicKinematicCloud();
-
-
-    // Member functions
-
-        //- Write fields
-        virtual void writeFields() const;
-};
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
+    typedef KinematicCloud<basicKinematicParcel> basicKinematicCloud;
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.C b/src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.C
deleted file mode 100644
index fc7ee7dc26..0000000000
--- a/src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.C
+++ /dev/null
@@ -1,68 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "basicThermoCloud.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(basicThermoCloud, 0);
-};
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::basicThermoCloud::basicThermoCloud
-(
-    const word& cloudName,
-    const volScalarField& rho,
-    const volVectorField& U,
-    const dimensionedVector& g,
-    basicThermo& thermo
-)
-:
-    ThermoCloud<basicThermoParcel>(cloudName, rho, U, g, thermo)
-{
-    basicThermoParcel::readFields(*this);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::basicThermoCloud::~basicThermoCloud()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::basicThermoCloud::writeFields() const
-{
-    basicThermoParcel::writeFields(*this);
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.H b/src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.H
index cbe3695d5d..c8ce18494f 100644
--- a/src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.H
@@ -43,60 +43,8 @@ SourceFiles
 
 namespace Foam
 {
-
-/*---------------------------------------------------------------------------*\
-                        Class basicThermoCloud Declaration
-\*---------------------------------------------------------------------------*/
-
-class basicThermoCloud
-:
-    public ThermoCloud<basicThermoParcel>
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        basicThermoCloud
-        (
-            const basicThermoCloud&
-        );
-
-        //- Disallow default bitwise assignment
-        void operator=(const basicThermoCloud&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("basicThermoCloud");
-
-
-    // Constructors
-
-        //- Construct given carrier gas fields
-        basicThermoCloud
-        (
-            const word& cloudName,
-            const volScalarField& rho,
-            const volVectorField& U,
-            const dimensionedVector& g,
-            basicThermo& thermo
-        );
-
-
-    //- Destructor
-    ~basicThermoCloud();
-
-
-    // Member Functions
-
-        //- Write fields
-        virtual void writeFields() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
+    typedef ThermoCloud<basicThermoParcel> basicThermoCloud;
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/intermediate/clouds/include/createReactingCloudTypes.H b/src/lagrangian/intermediate/clouds/include/createReactingCloudTypes.H
deleted file mode 100644
index 14f439e285..0000000000
--- a/src/lagrangian/intermediate/clouds/include/createReactingCloudTypes.H
+++ /dev/null
@@ -1,50 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef createReactingCloudTypes_H
-#define createReactingCloudTypes_H
-
-#include "thermoPhysicsTypes.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#define createReactingCloudType(CloudType)                                    \
-                                                                              \
-    createReactingCloudThermoType(CloudType, constGasThermoPhysics);          \
-    createReactingCloudThermoType(CloudType, gasThermoPhysics);               \
-    createReactingCloudThermoType(CloudType, icoPoly8ThermoPhysics);
-
-
-#define createReactingCloudThermoType(CloudType, ThermoType)                  \
-                                                                              \
-    defineTemplateTypeNameAndDebug(CloudType<ThermoType>, 0);
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index 973c97a2d8..b09cac4790 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -520,10 +520,10 @@ public:
         // I-O
 
             //- Read
-            static void readFields(KinematicCloud<ParcelType>& c);
+            static void readFields(Cloud<ParcelType>& c);
 
             //- Write
-            static void writeFields(const KinematicCloud<ParcelType>& c);
+            static void writeFields(const Cloud<ParcelType>& c);
 
 
     // Ostream Operator
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
index c69199cfb7..637de03109 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
@@ -103,10 +103,7 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
 
 
 template<class ParcelType>
-void Foam::KinematicParcel<ParcelType>::readFields
-(
-    KinematicCloud<ParcelType>& c
-)
+void Foam::KinematicParcel<ParcelType>::readFields(Cloud<ParcelType>& c)
 {
     if (!c.size())
     {
@@ -153,10 +150,7 @@ void Foam::KinematicParcel<ParcelType>::readFields
 
 
 template<class ParcelType>
-void Foam::KinematicParcel<ParcelType>::writeFields
-(
-    const KinematicCloud<ParcelType>& c
-)
+void Foam::KinematicParcel<ParcelType>::writeFields(const Cloud<ParcelType>& c)
 {
     Particle<ParcelType>::writeFields(c);
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index fabe2f6136..1bc1ee99d7 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -399,13 +399,10 @@ public:
         // I-O
 
             //- Read
-            static void readFields(ReactingMultiphaseCloud<ParcelType>& c);
+            static void readFields(Cloud<ParcelType>& c);
 
             //- Write
-            static void writeFields
-            (
-                const ReactingMultiphaseCloud<ParcelType>& c
-            );
+            static void writeFields(const Cloud<ParcelType>& c);
 
 
     // Ostream Operator
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
index 512f21e953..04747d9420 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
@@ -94,14 +94,17 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
 template<class ParcelType>
 void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
 (
-    ReactingMultiphaseCloud<ParcelType>& c
+    Cloud<ParcelType>& cIn
 )
 {
-    if (!c.size())
+    if (!cIn.size())
     {
         return;
     }
 
+    ReactingMultiphaseCloud<ParcelType>& c =
+        dynamic_cast<ReactingMultiphaseCloud<ParcelType>&>(cIn);
+
     ReactingParcel<ParcelType>::readFields(c);
 
     // Get names and sizes for each Y...
@@ -185,9 +188,12 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
 template<class ParcelType>
 void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
 (
-    const ReactingMultiphaseCloud<ParcelType>& c
+    const Cloud<ParcelType>& cIn
 )
 {
+    const ReactingMultiphaseCloud<ParcelType>& c =
+        dynamic_cast<const ReactingMultiphaseCloud<ParcelType>&>(cIn);
+
     ReactingParcel<ParcelType>::writeFields(c);
 
     label np =  c.size();
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index c3910a46ab..8219ec2ac7 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -348,10 +348,10 @@ public:
         // I-O
 
             //- Read
-            static void readFields(ReactingCloud<ParcelType>& c);
+            static void readFields(Cloud<ParcelType>& c);
 
             //- Write
-            static void writeFields(const ReactingCloud<ParcelType>& c);
+            static void writeFields(const Cloud<ParcelType>& c);
 
 
     // Ostream Operator
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
index b72d2af25a..f8bfbf24f7 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
@@ -88,16 +88,16 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
 
 
 template<class ParcelType>
-void Foam::ReactingParcel<ParcelType>::readFields
-(
-    ReactingCloud<ParcelType>& c
-)
+void Foam::ReactingParcel<ParcelType>::readFields(Cloud<ParcelType>& cIn)
 {
-    if (!c.size())
+    if (!cIn.size())
     {
         return;
     }
 
+    ReactingCloud<ParcelType>& c =
+        dynamic_cast<ReactingCloud<ParcelType>&>(cIn);
+
     ThermoParcel<ParcelType>::readFields(c);
 
     IOField<scalar> mass0(c.fieldIOobject("mass0", IOobject::MUST_READ));
@@ -152,12 +152,15 @@ void Foam::ReactingParcel<ParcelType>::readFields
 template<class ParcelType>
 void Foam::ReactingParcel<ParcelType>::writeFields
 (
-    const ReactingCloud<ParcelType>& c
+    const Cloud<ParcelType>& cIn
 )
 {
+    const ReactingCloud<ParcelType>& c =
+        dynamic_cast<const ReactingCloud<ParcelType>&>(cIn);
+
     ThermoParcel<ParcelType>::writeFields(c);
 
-    const label np =  c.size();
+    const label np = c.size();
 
     if (np > 0)
     {
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index f0e728bc1f..4bcdc43707 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -359,10 +359,10 @@ public:
         // I-O
 
             //- Read
-            static void readFields(ThermoCloud<ParcelType>& c);
+            static void readFields(Cloud<ParcelType>& c);
 
             //- Write
-            static void writeFields(const ThermoCloud<ParcelType>& c);
+            static void writeFields(const Cloud<ParcelType>& c);
 
 
     // Ostream Operator
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
index 827b2f36f4..b0a075af06 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
@@ -79,10 +79,7 @@ Foam::ThermoParcel<ParcelType>::ThermoParcel
 
 
 template<class ParcelType>
-void Foam::ThermoParcel<ParcelType>::readFields
-(
-    ThermoCloud<ParcelType>& c
-)
+void Foam::ThermoParcel<ParcelType>::readFields(Cloud<ParcelType>& c)
 {
     if (!c.size())
     {
@@ -111,10 +108,7 @@ void Foam::ThermoParcel<ParcelType>::readFields
 
 
 template<class ParcelType>
-void Foam::ThermoParcel<ParcelType>::writeFields
-(
-    const ThermoCloud<ParcelType>& c
-)
+void Foam::ThermoParcel<ParcelType>::writeFields(const Cloud<ParcelType>& c)
 {
     KinematicParcel<ParcelType>::writeFields(c);
 
diff --git a/src/lagrangian/intermediate/particleForces/particleForces.H b/src/lagrangian/intermediate/particleForces/particleForces.H
index fbe88a6dc8..f0be486381 100644
--- a/src/lagrangian/intermediate/particleForces/particleForces.H
+++ b/src/lagrangian/intermediate/particleForces/particleForces.H
@@ -26,7 +26,8 @@ Class
     Foam::particleForces
 
 Description
-    Particle forces
+    Provides a mechanism to calculate particle forces
+    Note: forces are force per unit mass (accelerations)
 
 SourceFiles
     particleForces.C
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
index 322ca30284..2589813289 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
@@ -58,7 +58,7 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::Sh
     const scalar Sc
 ) const
 {
-    return 2.0 + 0.6*Foam::sqrt(Re)*pow(Sc, 0.333);
+    return 2.0 + 0.6*Foam::sqrt(Re)*cbrt(Sc);
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.H b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.H
index 1f9015800c..72453a2124 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.H
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.H
@@ -115,42 +115,42 @@ public:
 
     // Member Functions
 
-    // Access
+        // Access
 
-        //- Return the cloud dictionary
-        const dictionary& dict() const;
+            //- Return the cloud dictionary
+            const dictionary& dict() const;
 
-        //- Return the coefficients dictionary
-        const dictionary& coeffDict() const;
+            //- Return the coefficients dictionary
+            const dictionary& coeffDict() const;
 
-        //- Return the owner cloud object
-        const CloudType& owner() const;
+            //- Return the owner cloud object
+            const CloudType& owner() const;
 
-        //- Return the Bird htc correction flag
-        const Switch& BirdCorrection() const;
+            //- Return the Bird htc correction flag
+            const Switch& BirdCorrection() const;
 
-        //- Flag to indicate whether model activates heat transfer model
-        virtual bool active() const = 0;
+            //- Flag to indicate whether model activates heat transfer model
+            virtual bool active() const = 0;
 
 
-    // Evaluation
+        // Evaluation
 
-        //- Nusselt number
-        virtual scalar Nu
-        (
-            const scalar Re,
-            const scalar Pr
-        ) const = 0;
+            //- Nusselt number
+            virtual scalar Nu
+            (
+                const scalar Re,
+                const scalar Pr
+            ) const = 0;
 
-        //- Return heat transfer coefficient
-        virtual scalar htc
-        (
-            const scalar dp,
-            const scalar Re,
-            const scalar Pr,
-            const scalar kappa,
-            const scalar NCpW
-        ) const;
+            //- Return heat transfer coefficient
+            virtual scalar htc
+            (
+                const scalar dp,
+                const scalar Re,
+                const scalar Pr,
+                const scalar kappa,
+                const scalar NCpW
+            ) const;
 };
 
 
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/NoHeatTransfer/NoHeatTransfer.H b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/NoHeatTransfer/NoHeatTransfer.H
index 2941b4c2b7..7a5ea0a86a 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/NoHeatTransfer/NoHeatTransfer.H
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/NoHeatTransfer/NoHeatTransfer.H
@@ -39,6 +39,7 @@ Description
 
 namespace Foam
 {
+
 /*---------------------------------------------------------------------------*\
                        Class NoHeatTransfer Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/RanzMarshall/RanzMarshall.H b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/RanzMarshall/RanzMarshall.H
index e30c51cef4..3df31f1f8d 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/RanzMarshall/RanzMarshall.H
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/RanzMarshall/RanzMarshall.H
@@ -39,6 +39,7 @@ Description
 
 namespace Foam
 {
+
 /*---------------------------------------------------------------------------*\
                          Class RanzMarshall Declaration
 \*---------------------------------------------------------------------------*/
@@ -76,9 +77,6 @@ public:
             //- Flag to indicate whether model activates heat transfer model
             virtual bool active() const;
 
-            //- Return the Bird correction flag
-            const Switch& BirdCorrection() const;
-
 
         // Evaluation
 
diff --git a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
index 788763b163..34dd65156a 100644
--- a/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
+++ b/src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.H
@@ -57,7 +57,7 @@ if (runTime.outputTime())
 
         averageTemperature =
         (
-            2.0/(moleculeCloud::kb * accumulatedDOFs)
+            2.0/(physicoChemical::k.value()*accumulatedDOFs)
             *
             (
                 accumulatedTotalLinearKE + accumulatedTotalAngularKE
@@ -70,10 +70,8 @@ if (runTime.outputTime())
         (
             (
                 (accumulatedNMols/nAveragingSteps)
-                *
-                moleculeCloud::kb * averageTemperature
-                +
-                accumulatedTotalrDotfSum/(6.0 * nAveragingSteps)
+               *physicoChemical::k.value()*averageTemperature
+              + accumulatedTotalrDotfSum/(6.0*nAveragingSteps)
             )
             /
             meshVolume
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
index ef2de145c4..990c4af6be 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
@@ -382,9 +382,12 @@ public:
             int&
         );
 
-        static void readFields(moleculeCloud& mC);
 
-        static void writeFields(const moleculeCloud& mC);
+    // I-O
+
+        static void readFields(Cloud<molecule>& mC);
+
+        static void writeFields(const Cloud<molecule>& mC);
 
 
     // IOstream Operators
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
index fa2e239304..fb732c77f1 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
@@ -98,7 +98,7 @@ Foam::molecule::molecule
 }
 
 
-void Foam::molecule::readFields(moleculeCloud& mC)
+void Foam::molecule::readFields(Cloud<molecule>& mC)
 {
     if (!mC.size())
     {
@@ -150,7 +150,7 @@ void Foam::molecule::readFields(moleculeCloud& mC)
 }
 
 
-void Foam::molecule::writeFields(const moleculeCloud& mC)
+void Foam::molecule::writeFields(const Cloud<molecule>& mC)
 {
     Particle<molecule>::writeFields(mC);
 
@@ -241,9 +241,10 @@ void Foam::molecule::writeFields(const moleculeCloud& mC)
     orientation2.write();
     orientation3.write();
 
-    mC.writeXYZ
+    const moleculeCloud& m = dynamic_cast<const moleculeCloud&>(mC);
+    m.writeXYZ
     (
-        mC.mesh().time().timePath() + "/lagrangian" + "/moleculeCloud.xmol"
+        m.mesh().time().timePath() + "/lagrangian" + "/moleculeCloud.xmol"
     );
 }
 
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
index 40ebfba487..d51e92f453 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -28,6 +28,8 @@ License
 #include "fvMesh.H"
 #include "mathConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
@@ -36,13 +38,6 @@ namespace Foam
     defineTemplateTypeNameAndDebug(Cloud<molecule>, 0);
 };
 
-Foam::scalar Foam::moleculeCloud::kb = 1.380650277e-23;
-
-Foam::scalar Foam::moleculeCloud::elementaryCharge = 1.602176487e-19;
-
-Foam::scalar Foam::moleculeCloud::vacuumPermittivity = 8.854187817e-12;
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::moleculeCloud::buildConstProps()
@@ -103,9 +98,7 @@ void Foam::moleculeCloud::buildConstProps()
 
 void Foam::moleculeCloud::setSiteSizesAndPositions()
 {
-    iterator mol(this->begin());
-
-    for (mol = this->begin(); mol != this->end(); ++mol)
+    forAllIter(moleculeCloud, *this, mol)
     {
         const molecule::constantProperties& cP = constProps(mol().id());
 
@@ -123,14 +116,7 @@ void Foam::moleculeCloud::buildCellOccupancy()
         cellOccupancy_[cO].clear();
     }
 
-    iterator mol(this->begin());
-
-    for
-    (
-        mol = this->begin();
-        mol != this->end();
-        ++mol
-    )
+    forAllIter(moleculeCloud, *this, mol)
     {
         cellOccupancy_[mol().cell()].append(&mol());
     }
@@ -225,9 +211,7 @@ void Foam::moleculeCloud::calculateTetherForce()
 {
     const tetherPotentialList& tetherPot(pot_.tetherPotentials());
 
-    iterator mol(this->begin());
-
-    for (mol = this->begin(); mol != this->end(); ++mol)
+    forAllIter(moleculeCloud, *this, mol)
     {
         if (mol().tethered())
         {
@@ -252,9 +236,7 @@ void Foam::moleculeCloud::calculateTetherForce()
 
 void Foam::moleculeCloud::calculateExternalForce()
 {
-    iterator mol(this->begin());
-
-    for (mol = this->begin(); mol != this->end(); ++mol)
+    forAllIter(moleculeCloud, *this, mol)
     {
         mol().a() += pot_.gravity();
     }
@@ -624,11 +606,11 @@ void Foam::moleculeCloud::initialiseMolecules
                     zoneDict.lookup("orientationAngles")
                 );
 
-                scalar phi(orientationAngles.x()*constant::math::pi/180.0);
+                scalar phi(orientationAngles.x()*math::pi/180.0);
 
-                scalar theta(orientationAngles.y()*constant::math::pi/180.0);
+                scalar theta(orientationAngles.y()*math::pi/180.0);
 
-                scalar psi(orientationAngles.z()*constant::math::pi/180.0);
+                scalar psi(orientationAngles.z()*math::pi/180.0);
 
                 const tensor R
                 (
@@ -994,11 +976,11 @@ void Foam::moleculeCloud::createMolecule
     {
         pi = equipartitionAngularMomentum(temperature, cP);
 
-        scalar phi(rndGen_.scalar01()*constant::math::twoPi);
+        scalar phi(rndGen_.scalar01()*math::twoPi);
 
-        scalar theta(rndGen_.scalar01()*constant::math::twoPi);
+        scalar theta(rndGen_.scalar01()*math::twoPi);
 
-        scalar psi(rndGen_.scalar01()*constant::math::twoPi);
+        scalar psi(rndGen_.scalar01()*math::twoPi);
 
         Q = tensor
         (
@@ -1039,9 +1021,7 @@ Foam::label Foam::moleculeCloud::nSites() const
 {
     label n = 0;
 
-    const_iterator mol(this->begin());
-
-    for (mol = this->begin(); mol != this->end(); ++mol)
+    forAllConstIter(moleculeCloud, *this, mol)
     {
         n += constProps(mol().id()).nSites();
     }
@@ -1055,7 +1035,8 @@ Foam::label Foam::moleculeCloud::nSites() const
 Foam::moleculeCloud::moleculeCloud
 (
     const polyMesh& mesh,
-    const potential& pot
+    const potential& pot,
+    bool readFields
 )
 :
     Cloud<molecule>(mesh, "moleculeCloud", false),
@@ -1066,7 +1047,10 @@ Foam::moleculeCloud::moleculeCloud
     constPropList_(),
     rndGen_(clock::getTime())
 {
-    molecule::readFields(*this);
+    if (readFields)
+    {
+        molecule::readFields(*this);
+    }
 
     buildConstProps();
 
@@ -1082,9 +1066,10 @@ Foam::moleculeCloud::moleculeCloud
 (
     const polyMesh& mesh,
     const potential& pot,
-    const IOdictionary& mdInitialiseDict
+    const IOdictionary& mdInitialiseDict,
+    bool readFields
 )
-    :
+:
     Cloud<molecule>(mesh, "moleculeCloud", false),
     mesh_(mesh),
     pot_(pot),
@@ -1092,7 +1077,10 @@ Foam::moleculeCloud::moleculeCloud
     constPropList_(),
     rndGen_(clock::getTime())
 {
-    molecule::readFields(*this);
+    if (readFields)
+    {
+        molecule::readFields(*this);
+    }
 
     clear();
 
@@ -1126,10 +1114,8 @@ void Foam::moleculeCloud::calculateForce()
 {
     buildCellOccupancy();
 
-    iterator mol(this->begin());
-
     // Set accumulated quantities to zero
-    for (mol = this->begin(); mol != this->end(); ++mol)
+    forAllIter(moleculeCloud, *this, mol)
     {
         mol().siteForces() = vector::zero;
 
@@ -1166,9 +1152,7 @@ void Foam::moleculeCloud::applyConstraintsAndThermostats
         << "----------------------------------------"
         << endl;
 
-    iterator mol(this->begin());
-
-    for (mol = this->begin(); mol != this->end(); ++mol)
+    forAllIter(moleculeCloud, *this, mol)
     {
         mol().v() *= temperatureCorrectionFactor;
 
@@ -1177,21 +1161,13 @@ void Foam::moleculeCloud::applyConstraintsAndThermostats
 }
 
 
-void Foam::moleculeCloud::writeFields() const
-{
-    molecule::writeFields(*this);
-}
-
-
 void Foam::moleculeCloud::writeXYZ(const fileName& fName) const
 {
     OFstream str(fName);
 
     str << nSites() << nl << "moleculeCloud site positions in angstroms" << nl;
 
-    const_iterator mol(this->begin());
-
-    for (mol = this->begin(); mol != this->end(); ++mol)
+    forAllConstIter(moleculeCloud, *this, mol)
     {
         const molecule::constantProperties& cP = constProps(mol().id());
 
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
index 6863227c12..81623e715f 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
@@ -156,22 +156,14 @@ private:
 
 public:
 
-    // Static data members
-
-        static scalar kb;
-
-        static scalar elementaryCharge;
-
-        static scalar vacuumPermittivity;
-
-
     // Constructors
 
         //- Construct given mesh and potential references
         moleculeCloud
         (
             const polyMesh& mesh,
-            const potential& pot
+            const potential& pot,
+            bool readFields = true
         );
 
         //- Construct given mesh, potential and mdInitialiseDict
@@ -179,7 +171,8 @@ public:
         (
             const polyMesh& mesh,
             const potential& pot,
-            const IOdictionary& mdInitialiseDict
+            const IOdictionary& mdInitialiseDict,
+            bool readFields = true
         );
 
 
@@ -217,9 +210,6 @@ public:
 
     // Member Operators
 
-        //- Write fields
-        void writeFields() const;
-
         //- Write molecule sites in XYZ format
         void writeXYZ(const fileName& fName) const;
 };
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
index 94c9f49286..fee67cfced 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
@@ -24,6 +24,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "constants.H"
+
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 inline void Foam::moleculeCloud::evaluatePair
@@ -576,7 +580,7 @@ inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
     scalar mass
 )
 {
-    return sqrt(kb*temperature/mass)*vector
+    return sqrt(physicoChemical::k.value()*temperature/mass)*vector
     (
         rndGen_.GaussNormal(),
         rndGen_.GaussNormal(),
@@ -591,7 +595,7 @@ inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
     const molecule::constantProperties& cP
 )
 {
-    scalar sqrtKbT = sqrt(kb*temperature);
+    scalar sqrtKbT = sqrt(physicoChemical::k.value()*temperature);
 
     if (cP.linearMolecule())
     {
diff --git a/src/lagrangian/solidParticle/solidParticleCloud.C b/src/lagrangian/solidParticle/solidParticleCloud.C
index b317f75d94..e17f4efc60 100644
--- a/src/lagrangian/solidParticle/solidParticleCloud.C
+++ b/src/lagrangian/solidParticle/solidParticleCloud.C
@@ -42,7 +42,8 @@ namespace Foam
 Foam::solidParticleCloud::solidParticleCloud
 (
     const fvMesh& mesh,
-    const word& cloudName
+    const word& cloudName,
+    bool readFields
 )
 :
     Cloud<solidParticle>(mesh, cloudName, false),
@@ -62,7 +63,10 @@ Foam::solidParticleCloud::solidParticleCloud
     e_(dimensionedScalar(particleProperties_.lookup("e")).value()),
     mu_(dimensionedScalar(particleProperties_.lookup("mu")).value())
 {
-    solidParticle::readFields(*this);
+    if (readFields)
+    {
+        solidParticle::readFields(*this);
+    }
 }
 
 
@@ -84,10 +88,4 @@ void Foam::solidParticleCloud::move(const dimensionedVector& g)
 }
 
 
-void Foam::solidParticleCloud::writeFields() const
-{
-    solidParticle::writeFields(*this);
-}
-
-
 // ************************************************************************* //
diff --git a/src/lagrangian/solidParticle/solidParticleCloud.H b/src/lagrangian/solidParticle/solidParticleCloud.H
index f171007ce3..81675b34cd 100644
--- a/src/lagrangian/solidParticle/solidParticleCloud.H
+++ b/src/lagrangian/solidParticle/solidParticleCloud.H
@@ -86,7 +86,8 @@ public:
         solidParticleCloud
         (
             const fvMesh&,
-            const word& cloudName = "defaultCloud"
+            const word& cloudName = "defaultCloud",
+            bool readFields = true
         );
 
 
@@ -106,12 +107,6 @@ public:
             //- Move the particles under the influence of the given
             //  gravitational acceleration
             void move(const dimensionedVector& g);
-
-
-        // Write
-
-            //- Write fields
-            virtual void writeFields() const;
 };
 
 
diff --git a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
index dab51b4714..e3f4983ffe 100644
--- a/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
+++ b/src/meshTools/triSurface/triangleFuncs/triangleFuncs.C
@@ -115,7 +115,7 @@ bool Foam::triangleFuncs::intersectAxesBundle
     //          V10:(-1.285715 8.99165e-16 -1.142858)
     //          V20:(0 0 -1.678573)
     //          i0:0
-    if (Foam::mag(det)/localScale < SMALL)
+    if (localScale < VSMALL || Foam::mag(det)/localScale < SMALL)
     {
         // Triangle parallel to dir
         return false;
diff --git a/src/postProcessing/functionObjects/Allwmake b/src/postProcessing/functionObjects/Allwmake
index b068fe51cb..6c642dfe7b 100755
--- a/src/postProcessing/functionObjects/Allwmake
+++ b/src/postProcessing/functionObjects/Allwmake
@@ -7,6 +7,5 @@ wmake libso forces
 wmake libso IO
 wmake libso utilities
 wmake libso systemCall
-wmake libso zones
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
index 03b5a80288..8754ed2df4 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
@@ -387,10 +387,7 @@ void Foam::streamLineParticle::hitPatch
 {}
 
 
-void Foam::streamLineParticle::readFields
-(
-    Cloud<streamLineParticle>& c
-)
+void Foam::streamLineParticle::readFields(Cloud<streamLineParticle>& c)
 {
     if (!c.size())
     {
@@ -426,10 +423,7 @@ void Foam::streamLineParticle::readFields
 }
 
 
-void Foam::streamLineParticle::writeFields
-(
-    const Cloud<streamLineParticle>& c
-)
+void Foam::streamLineParticle::writeFields(const Cloud<streamLineParticle>& c)
 {
     Particle<streamLineParticle>::writeFields(c);
 
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticleCloud.C b/src/postProcessing/functionObjects/field/streamLine/streamLineParticleCloud.C
index 95a92d024a..51f16c91f9 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticleCloud.C
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticleCloud.C
@@ -26,15 +26,11 @@ License
 
 #include "streamLineParticleCloud.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTemplateTypeNameAndDebug(Cloud<streamLineParticle>, 0);
-
+    defineTemplateTypeNameAndDebug(Cloud<streamLineParticle>, 0);
 };
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -42,12 +38,16 @@ defineTemplateTypeNameAndDebug(Cloud<streamLineParticle>, 0);
 Foam::streamLineParticleCloud::streamLineParticleCloud
 (
     const polyMesh& mesh,
-    const word& cloudName
+    const word& cloudName,
+    bool readFields
 )
 :
     Cloud<streamLineParticle>(mesh, cloudName, false)
 {
-    readFields();
+    if (readFields)
+    {
+        streamLineParticle::readFields(*this);
+    }
 }
 
 
@@ -62,18 +62,4 @@ Foam::streamLineParticleCloud::streamLineParticleCloud
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::streamLineParticleCloud::readFields()
-{
-    streamLineParticle::readFields(*this);
-}
-
-
-void Foam::streamLineParticleCloud::writeFields() const
-{
-    streamLineParticle::writeFields(*this);
-}
-
-
 // ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticleCloud.H b/src/postProcessing/functionObjects/field/streamLine/streamLineParticleCloud.H
index 40171fca9d..5894124277 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticleCloud.H
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticleCloud.H
@@ -45,7 +45,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class streamLineCloud Declaration
+                      Class streamLineCloud Declaration
 \*---------------------------------------------------------------------------*/
 
 class streamLineParticleCloud
@@ -72,7 +72,8 @@ public:
         streamLineParticleCloud
         (
             const polyMesh&,
-            const word& cloudName = "defaultCloud"
+            const word& cloudName = "defaultCloud",
+            bool readFields = true
         );
 
         //- Construct from mesh, cloud name, and a list of particles
@@ -82,15 +83,6 @@ public:
             const word& cloudName,
             const IDLList<streamLineParticle>& particles
         );
-
-
-    // Member Functions
-
-        //- Read fields
-        virtual void readFields();
-
-        //- Write fields
-        virtual void writeFields() const;
 };
 
 
diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
index 8957dee17c..eac6499043 100644
--- a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
+++ b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
@@ -29,6 +29,10 @@ License
 #include "dictionary.H"
 #include "dsmcCloud.H"
 
+#include "constants.H"
+
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
@@ -164,7 +168,7 @@ void Foam::dsmcFields::write()
                     obr_,
                     IOobject::NO_READ
                 ),
-                2.0/(3.0*dsmcCloud::kb*rhoNMean)
+                2.0/(3.0*physicoChemical::k.value()*rhoNMean)
                 *(linearKEMean - 0.5*rhoMMean*(UMean & UMean))
             );
 
@@ -178,7 +182,7 @@ void Foam::dsmcFields::write()
                     obr_,
                     IOobject::NO_READ
                 ),
-                2.0/(dsmcCloud::kb*iDofMean)*internalEMean
+                2.0/(physicoChemical::k.value()*iDofMean)*internalEMean
             );
 
             Info<< "    Calculating overallT field." << endl;
@@ -191,7 +195,7 @@ void Foam::dsmcFields::write()
                     obr_,
                     IOobject::NO_READ
                 ),
-                2.0/(dsmcCloud::kb*(3.0*rhoNMean + iDofMean))
+                2.0/(physicoChemical::k.value()*(3.0*rhoNMean + iDofMean))
                 *(linearKEMean - 0.5*rhoMMean*(UMean & UMean) + internalEMean)
             );
 
diff --git a/src/postProcessing/functionObjects/zones/Make/files b/src/postProcessing/functionObjects/zones/Make/files
deleted file mode 100644
index 90e33f7524..0000000000
--- a/src/postProcessing/functionObjects/zones/Make/files
+++ /dev/null
@@ -1,4 +0,0 @@
-faceZoneIntegration/faceZonesIntegration.C
-faceZoneIntegration/faceZonesIntegrationFunctionObject.C
-
-LIB = $(FOAM_LIBBIN)/libzoneFunctionObjects
diff --git a/src/postProcessing/functionObjects/zones/Make/options b/src/postProcessing/functionObjects/zones/Make/options
deleted file mode 100644
index 5166bcc9e3..0000000000
--- a/src/postProcessing/functionObjects/zones/Make/options
+++ /dev/null
@@ -1,9 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
-
-LIB_LIBS = \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lsampling
diff --git a/src/postProcessing/functionObjects/zones/faceZoneIntegration/IOfaceZonesIntegration.H b/src/postProcessing/functionObjects/zones/faceZoneIntegration/IOfaceZonesIntegration.H
deleted file mode 100644
index ef1ecbf004..0000000000
--- a/src/postProcessing/functionObjects/zones/faceZoneIntegration/IOfaceZonesIntegration.H
+++ /dev/null
@@ -1,50 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Typedef
-    Foam::IOfaceZonesIntegration
-
-Description
-    Instance of the generic IOOutputFilter for faceZonesIntegration.
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef IOfaceZonesIntegration_H
-#define IOfaceZonesIntegration_H
-
-#include "faceZonesIntegration.H"
-#include "IOOutputFilter.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    typedef IOOutputFilter<faceZonesIntegration> IOfaceZonesIntegration;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.C b/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.C
deleted file mode 100644
index 2899224d29..0000000000
--- a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.C
+++ /dev/null
@@ -1,330 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "faceZonesIntegration.H"
-#include "volFields.H"
-#include "dictionary.H"
-#include "Time.H"
-#include "IOmanip.H"
-#include "ListListOps.H"
-#include "processorPolyPatch.H"
-#include "cyclicPolyPatch.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(faceZonesIntegration, 0);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::faceZonesIntegration::faceZonesIntegration
-(
-    const word& name,
-    const objectRegistry& obr,
-    const dictionary& dict,
-    const bool loadFromFiles
-)
-:
-    name_(name),
-    obr_(obr),
-    active_(true),
-    log_(false),
-    zoneNames_(),
-    fieldNames_(),
-    filePtr_(NULL)
-{
-    // Check if the available mesh is an fvMesh otherise deactivate
-    if (!isA<fvMesh>(obr_))
-    {
-        active_ = false;
-        WarningIn
-        (
-            "Foam::faceZonesIntegration::faceZonesIntegration"
-            "("
-                "const word&, "
-                "const objectRegistry&, "
-                "const dictionary&, "
-                "const bool"
-            ")"
-        )   << "No fvMesh available, deactivating."
-            << endl;
-    }
-
-    read(dict);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::faceZonesIntegration::~faceZonesIntegration()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::faceZonesIntegration::read(const dictionary& dict)
-{
-    if (active_)
-    {
-        log_ = dict.lookupOrDefault<Switch>("log", false);
-
-        dict.lookup("fields") >> fieldNames_;
-
-        dict.lookup("faceZones") >> zoneNames_;
-    }
-}
-
-
-void Foam::faceZonesIntegration::makeFile()
-{
-    // Create the face Zone file if not already created
-    if (filePtr_.empty())
-    {
-        if (debug)
-        {
-            Info<< "Creating faceZonesIntegration file." << endl;
-        }
-
-        // File update
-        if (Pstream::master())
-        {
-            fileName faceZonesIntegrationDir;
-            if (Pstream::parRun())
-            {
-                // Put in undecomposed case (Note: gives problems for
-                // distributed data running)
-                faceZonesIntegrationDir =
-                    obr_.time().path()/".."/name_/obr_.time().timeName();
-            }
-            else
-            {
-                faceZonesIntegrationDir =
-                    obr_.time().path()/name_/obr_.time().timeName();
-            }
-
-            // Create directory if does not exist.
-            mkDir(faceZonesIntegrationDir);
-
-            // Open new file at start up
-            filePtr_.resize(fieldNames_.size());
-
-            forAll(fieldNames_, fieldI)
-            {
-                const word& fieldName = fieldNames_[fieldI];
-
-                OFstream* sPtr = new OFstream
-                    (
-                        faceZonesIntegrationDir/fieldName
-                    );
-
-                filePtr_.insert(fieldName, sPtr);
-            }
-
-            // Add headers to output data
-            writeFileHeader();
-        }
-    }
-}
-
-
-void Foam::faceZonesIntegration::writeFileHeader()
-{
-    forAllIter(HashPtrTable<OFstream>, filePtr_, iter)
-    {
-        unsigned int w = IOstream::defaultPrecision() + 7;
-
-        OFstream& os = *filePtr_[iter.key()];
-
-        os  << "#Time " << setw(w);
-
-        forAll (zoneNames_, zoneI)
-        {
-            os  << zoneNames_[zoneI] << setw(w);
-        }
-
-        os  << nl << endl;
-    }
-}
-
-
-void Foam::faceZonesIntegration::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::faceZonesIntegration::end()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::faceZonesIntegration::write()
-{
-    if (active_)
-    {
-        makeFile();
-
-        forAll(fieldNames_, fieldI)
-        {
-            const word& fieldName = fieldNames_[fieldI];
-
-            const surfaceScalarField& sField =
-                obr_.lookupObject<surfaceScalarField>(fieldName);
-
-            const fvMesh& mesh = sField.mesh();
-
-            // 1. integrate over all face zones
-
-            scalarField integralVals(zoneNames_.size());
-
-            forAll(integralVals, zoneI)
-            {
-                const word& name = zoneNames_[zoneI];
-
-                label zoneID = mesh.faceZones().findZoneID(name);
-
-                const faceZone& fZone = mesh.faceZones()[zoneID];
-
-                integralVals[zoneI] = returnReduce
-                (
-                    calcIntegral(sField, fZone),
-                    sumOp<scalar>()
-                );
-            }
-
-
-            unsigned int w = IOstream::defaultPrecision() + 7;
-
-            // 2. Write only on master
-
-            if (Pstream::master() && filePtr_.found(fieldName))
-            {
-                OFstream& os = *filePtr_(fieldName);
-
-                os  << obr_.time().value();
-
-                forAll(integralVals, zoneI)
-                {
-                    os  << ' ' << setw(w) << integralVals[zoneI];
-
-                    if (log_)
-                    {
-                        Info<< "faceZonesIntegration output:" << nl
-                            << "    Integration[" << zoneI << "] "
-                            << integralVals[zoneI] << endl;
-                    }
-                }
-
-                os  << endl;
-            }
-        }
-    }
-}
-
-
-Foam::scalar Foam::faceZonesIntegration::calcIntegral
-(
-    const surfaceScalarField& sField,
-    const faceZone& fZone
-) const
-{
-    scalar sum = 0.0;
-    const fvMesh& mesh = sField.mesh();
-
-    forAll (fZone, i)
-    {
-        label faceI = fZone[i];
-
-        if (mesh.isInternalFace(faceI))
-        {
-            if (fZone.flipMap()[i])
-            {
-                sum -= sField[faceI];
-            }
-            else
-            {
-                sum += sField[faceI];
-            }
-        }
-        else
-        {
-            label patchI = mesh.boundaryMesh().whichPatch(faceI);
-            const polyPatch& pp = mesh.boundaryMesh()[patchI];
-            const fvsPatchScalarField& bField = sField.boundaryField()[patchI];
-            if (isA<processorPolyPatch>(pp))
-            {
-                if (refCast<const processorPolyPatch>(pp).owner())
-                {
-                    label patchFaceI = pp.whichFace(faceI);
-                    if (fZone.flipMap()[i])
-                    {
-                        sum -= bField[patchFaceI];
-                    }
-                    else
-                    {
-                        sum += bField[patchFaceI];
-                    }
-                }
-            }
-            else if (isA<cyclicPolyPatch>(pp))
-            {
-                label patchFaceI = faceI - pp.start();
-                if (patchFaceI < pp.size()/2)
-                {
-                    if (fZone.flipMap()[i])
-                    {
-                        sum -= bField[patchFaceI];
-                    }
-                    else
-                    {
-                        sum += bField[patchFaceI];
-                    }
-                }
-            }
-            else if (!isA<emptyPolyPatch>(pp))
-            {
-                label patchFaceI = faceI - pp.start();
-                if (fZone.flipMap()[i])
-                {
-                    sum -= bField[patchFaceI];
-                }
-                else
-                {
-                    sum += bField[patchFaceI];
-                }
-            }
-        }
-    }
-
-    return sum;
-}
-
-
-// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.H b/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.H
deleted file mode 100644
index d1a1bacfd2..0000000000
--- a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.H
+++ /dev/null
@@ -1,180 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::faceZonesIntegration
-
-Description
-    Integrates surfaceScalarFields on faceZones
-
-SourceFiles
-    faceZonesIntegration.C
-    IOfaceZonesIntegration.H
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef faceZonesIntegration_H
-#define faceZonesIntegration_H
-
-#include "fvCFD.H"
-#include "primitiveFieldsFwd.H"
-#include "volFieldsFwd.H"
-#include "HashPtrTable.H"
-#include "OFstream.H"
-#include "Switch.H"
-#include "pointFieldFwd.H"
-#include "polyMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class objectRegistry;
-class dictionary;
-class mapPolyMesh;
-
-/*---------------------------------------------------------------------------*\
-                    Class faceZonesIntegration Declaration
-\*---------------------------------------------------------------------------*/
-
-class faceZonesIntegration
-{
-
-protected:
-
-    // Protected data
-
-        //- Name of this set of face zone integration,
-        //  Also used as the name of the output directory.
-        word name_;
-
-        const objectRegistry& obr_;
-
-
-        // Read from dictionary
-
-            //- On/off switch
-            bool active_;
-
-            //- Switch to send output to Info as well as to file
-            Switch log_;
-
-            //- List of face zone names to integrate over
-            wordList zoneNames_;
-
-            //- Names of the surface fields
-            wordList fieldNames_;
-
-
-        //- Current open files
-        HashPtrTable<OFstream> filePtr_;
-
-
-
-    // Protected Member Functions
-
-        //- If the integration file has not been created create it
-        void makeFile();
-
-        scalar calcIntegral
-        (
-            const surfaceScalarField& sField,
-            const faceZone& fZone
-        ) const;
-
-
-        //- Disallow default bitwise copy construct
-        faceZonesIntegration(const faceZonesIntegration&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const faceZonesIntegration&);
-
-        //- Output file header information
-        virtual void writeFileHeader();
-
-
-public:
-
-    //- Runtime type information
-    TypeName("faceZonesIntegration");
-
-
-    // Constructors
-
-        //- Construct for given objectRegistry and dictionary.
-        //  Allow the possibility to load fields from files
-        faceZonesIntegration
-        (
-            const word& name,
-            const objectRegistry&,
-            const dictionary&,
-            const bool loadFromFiles = false
-        );
-
-
-    //- Destructor
-    virtual ~faceZonesIntegration();
-
-
-    // Member Functions
-
-        //- Return name of the set of zones
-        virtual const word& name() const
-        {
-            return name_;
-        };
-
-        //- Read the zone integration data
-        virtual void read(const dictionary&);
-
-        //- Execute, currently does nothing
-        virtual void execute();
-
-        //- Execute at the final time-loop, currently does nothing
-        virtual void end();
-
-        //- Write the integration
-        virtual void write();
-
-        //- Update for changes of mesh
-        virtual void updateMesh(const mapPolyMesh&)
-        {}
-
-        //- Update for changes of mesh
-        virtual void movePoints(const pointField&)
-        {}
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegrationFunctionObject.C b/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegrationFunctionObject.C
deleted file mode 100644
index 43665119b1..0000000000
--- a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegrationFunctionObject.C
+++ /dev/null
@@ -1,43 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "faceZonesIntegrationFunctionObject.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineNamedTemplateTypeNameAndDebug(faceZonesIntegrationFunctionObject, 0);
-
-    addToRunTimeSelectionTable
-    (
-        functionObject,
-        faceZonesIntegrationFunctionObject,
-        dictionary
-    );
-}
-
-// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegrationFunctionObject.H b/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegrationFunctionObject.H
deleted file mode 100644
index 2cc95e2af2..0000000000
--- a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegrationFunctionObject.H
+++ /dev/null
@@ -1,55 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Typedef
-    Foam::faceZonesIntegrationFunctionObject
-
-Description
-    FunctionObject wrapper around faceZonesIntegration to allow them to be
-    created via the functions list within controlDict.
-
-SourceFiles
-    faceZonesIntegrationFunctionObject.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef faceZonesIntegrationFunctionObject_H
-#define faceZonesIntegrationFunctionObject_H
-
-#include "faceZonesIntegration.H"
-#include "OutputFilterFunctionObject.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    typedef OutputFilterFunctionObject<faceZonesIntegration>
-        faceZonesIntegrationFunctionObject;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/thermophysicalModels/pdfs/RosinRammler/RosinRammler.C b/src/thermophysicalModels/pdfs/RosinRammler/RosinRammler.C
index 36b84679ed..f157c2645f 100644
--- a/src/thermophysicalModels/pdfs/RosinRammler/RosinRammler.C
+++ b/src/thermophysicalModels/pdfs/RosinRammler/RosinRammler.C
@@ -47,14 +47,15 @@ Foam::RosinRammler::RosinRammler(const dictionary& dict, Random& rndGen)
     d_(pdfDict_.lookup("d")),
     n_(pdfDict_.lookup("n")),
     ls_(d_),
-    range_(maxValue_-minValue_)
+    range_(maxValue_ - minValue_)
 {
     if (minValue_<0)
     {
         FatalErrorIn
         (
             "RosinRammler::RosinRammler(const dictionary& dict)"
-        ) << " minValue = " << minValue_ << ", it must be >0." << abort(FatalError);
+        )   << " minValue = " << minValue_ << ", it must be >0." << nl
+            << abort(FatalError);
     }
 
     if (maxValue_<minValue_)
@@ -62,7 +63,8 @@ Foam::RosinRammler::RosinRammler(const dictionary& dict, Random& rndGen)
         FatalErrorIn
         (
             "RosinRammler::RosinRammler(const dictionary& dict)"
-        ) << " maxValue is smaller than minValue." << abort(FatalError);
+        )   << " maxValue is smaller than minValue." << nl
+            << abort(FatalError);
     }
 
     // find max value so that it can be normalized to 1.0
@@ -103,26 +105,20 @@ Foam::scalar Foam::RosinRammler::sample() const
 {
     scalar y = 0;
     scalar x = 0;
-    label n = d_.size();
-    bool success = false;
+    scalar p = 0.0;
 
-    while (!success)
+    do
     {
         x = minValue_ + range_*rndGen_.scalar01();
         y = rndGen_.scalar01();
-        scalar p = 0.0;
+        p = 0.0;
 
-        for (label i=0; i<n; i++)
+        forAll(d_, i)
         {
             scalar xx = pow(x/d_[i], n_[i]);
             p += ls_[i]*xx*exp(-xx);
         }
-
-        if (y<p)
-        {
-            success = true;
-        }
-    }
+    } while(y>p);
 
     return x;
 }
diff --git a/src/thermophysicalModels/pdfs/exponential/exponential.C b/src/thermophysicalModels/pdfs/exponential/exponential.C
index 973161537f..6ebe918fe8 100644
--- a/src/thermophysicalModels/pdfs/exponential/exponential.C
+++ b/src/thermophysicalModels/pdfs/exponential/exponential.C
@@ -46,14 +46,15 @@ Foam::exponential::exponential(const dictionary& dict, Random& rndGen)
     maxValue_(readScalar(pdfDict_.lookup("maxValue"))),
     lambda_(pdfDict_.lookup("lambda")),
     ls_(lambda_),
-    range_(maxValue_-minValue_)
+    range_(maxValue_ - minValue_)
 {
     if (minValue_<0)
     {
         FatalErrorIn
         (
             "exponential::exponential(const dictionary& dict)"
-        ) << " minValue = " << minValue_ << ", it must be >0." << abort(FatalError);
+        )   << " minValue = " << minValue_ << ", it must be >0." << nl
+            << abort(FatalError);
     }
 
     scalar sMax = 0;
@@ -73,7 +74,7 @@ Foam::exponential::exponential(const dictionary& dict, Random& rndGen)
         sMax = max(sMax, s);
     }
 
-    for(label i=0; i<n; i++)
+    for (label i=0; i<n; i++)
     {
         ls_[i] /= sMax;
     }
@@ -90,27 +91,21 @@ Foam::exponential::~exponential()
 
 Foam::scalar Foam::exponential::sample() const
 {
-    scalar y = 0;
-    scalar x = 0;
-    label n = lambda_.size();
-    bool success = false;
+    scalar y = 0.0;
+    scalar x = 0.0;
+    scalar p = 0.0;
 
-    while (!success)
+    do
     {
         x = minValue_ + range_*rndGen_.scalar01();
         y = rndGen_.scalar01();
-        scalar p = 0.0;
+        p = 0.0;
 
-        for(label i=0; i<n; i++)
+        forAll(lambda_, i)
         {
             p += ls_[i]*exp(-lambda_[i]*x);
         }
-
-        if (y<p)
-        {
-            success = true;
-        }
-    }
+    } while(p>y);
 
     return x;
 }
diff --git a/src/thermophysicalModels/pdfs/general/general.C b/src/thermophysicalModels/pdfs/general/general.C
index ec5d2b7abf..59b1a75953 100644
--- a/src/thermophysicalModels/pdfs/general/general.C
+++ b/src/thermophysicalModels/pdfs/general/general.C
@@ -45,7 +45,7 @@ Foam::general::general(const dictionary& dict, Random& rndGen)
     nEntries_(xy_.size()),
     minValue_(xy_[0][0]),
     maxValue_(xy_[nEntries_-1][0]),
-    range_(maxValue_-minValue_)
+    range_(maxValue_ - minValue_)
 {
     // normalize the pdf
     scalar yMax = 0;
@@ -72,12 +72,11 @@ Foam::general::~general()
 
 Foam::scalar Foam::general::sample() const
 {
-    scalar y = 0;
-    scalar x = 0;
+    scalar x = 0.0;
+    scalar y = 0.0;
+    scalar p = 0.0;
 
-    bool success = false;
-
-    while (!success)
+    do
     {
         x = minValue_ + range_*rndGen_.scalar01();
         y = rndGen_.scalar01();
@@ -87,23 +86,17 @@ Foam::scalar Foam::general::sample() const
         while (!intervalFound)
         {
             i++;
-            if ( (x>xy_[i][0]) && (x<xy_[i+1][0]) )
+            if ((x>xy_[i][0]) && (x<xy_[i+1][0]))
             {
                 intervalFound = true;
             }
         }
 
-        scalar p =
-            xy_[i][1]
-          + (x-xy_[i][0])
-           *(xy_[i+1][1]-xy_[i][1])
-           /(xy_[i+1][0]-xy_[i][0]);
-
-        if (y<p)
-        {
-            success = true;
-        }
-    }
+        p = xy_[i][1]
+          + (x - xy_[i][0])
+           *(xy_[i+1][1] - xy_[i][1])
+           /(xy_[i+1][0] - xy_[i][0]);
+    } while(p>y);
 
     return x;
 }
diff --git a/src/thermophysicalModels/pdfs/normal/normal.C b/src/thermophysicalModels/pdfs/normal/normal.C
index f8633c270b..4a66fc9883 100644
--- a/src/thermophysicalModels/pdfs/normal/normal.C
+++ b/src/thermophysicalModels/pdfs/normal/normal.C
@@ -46,7 +46,7 @@ Foam::normal::normal(const dictionary& dict, Random& rndGen)
     expectation_(pdfDict_.lookup("expectation")),
     variance_(pdfDict_.lookup("variance")),
     strength_(pdfDict_.lookup("strength")),
-    range_(maxValue_-minValue_)
+    range_(maxValue_ - minValue_)
 {
     scalar sMax = 0;
     label n = strength_.size();
@@ -84,31 +84,25 @@ Foam::normal::~normal()
 
 Foam::scalar Foam::normal::sample() const
 {
-    scalar y = 0;
-    scalar x = 0;
-    label n = expectation_.size();
-    bool success = false;
+    scalar x = 0.0;
+    scalar y = 0.0;
+    scalar p = 0.0;
 
-    while (!success)
+    do
     {
         x = minValue_ + range_*rndGen_.scalar01();
         y = rndGen_.scalar01();
-        scalar p = 0.0;
+        p = 0.0;
 
-        for (label i=0; i<n; i++)
+        forAll(expectation_, i)
         {
             scalar nu = expectation_[i];
             scalar sigma = variance_[i];
             scalar s = strength_[i];
-            scalar v = (x-nu)/sigma;
+            scalar v = (x - nu)/sigma;
             p += s*exp(-0.5*v*v);
         }
-
-        if (y<p)
-        {
-            success = true;
-        }
-    }
+    } while (p>y);
 
     return x;
 }
diff --git a/src/thermophysicalModels/pdfs/uniform/uniform.C b/src/thermophysicalModels/pdfs/uniform/uniform.C
index 988c337071..09c2c4b05a 100644
--- a/src/thermophysicalModels/pdfs/uniform/uniform.C
+++ b/src/thermophysicalModels/pdfs/uniform/uniform.C
@@ -43,7 +43,7 @@ Foam::uniform::uniform(const dictionary& dict, Random& rndGen)
     pdfDict_(dict.subDict(typeName + "PDF")),
     minValue_(readScalar(pdfDict_.lookup("minValue"))),
     maxValue_(readScalar(pdfDict_.lookup("maxValue"))),
-    range_(maxValue_-minValue_)
+    range_(maxValue_ - minValue_)
 {}
 
 
diff --git a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C
index 1d4054fb04..f0b77bb1f2 100644
--- a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C
+++ b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C
@@ -176,7 +176,7 @@ tmp<scalarField> RASModel::yPlus(const label patchNo, const scalar Cmu) const
 
     if (isType<wallFvPatch>(curPatch))
     {
-        Yp = pow(Cmu, 0.25)
+        Yp = pow025(Cmu)
             *y_[patchNo]
             *sqrt(k()().boundaryField()[patchNo].patchInternalField())
            /(
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
index 5893f2072e..bacfe9615e 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
@@ -118,7 +118,7 @@ void turbulentMixingLengthFrequencyInletFvPatchScalarField::updateCoeffs()
     const scalar Cmu =
         rasModel.coeffDict().lookupOrDefault<scalar>("Cmu", 0.09);
 
-    const scalar Cmu25 = pow(Cmu, 0.25);
+    const scalar Cmu25 = pow025(Cmu);
 
     const fvPatchField<scalar>& kp =
         patch().lookupPatchField<volScalarField, scalar>(kName_);
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C
index e8fdf3edd5..b5022e2756 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C
@@ -204,7 +204,7 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
 
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
 
-    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu25 = pow025(Cmu_);
 
     const scalarField& y = rasModel.y()[patchI];
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
index 07382bad52..a129ec8bff 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
@@ -188,7 +188,7 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
 
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
 
-    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu25 = pow025(Cmu_);
     const scalar Cmu75 = pow(Cmu_, 0.75);
 
     const scalarField& y = rasModel.y()[patchI];
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutkRoughWallFunction/mutkRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutkRoughWallFunction/mutkRoughWallFunctionFvPatchScalarField.C
index 6af9687f7b..67353e4319 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutkRoughWallFunction/mutkRoughWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutkRoughWallFunction/mutkRoughWallFunctionFvPatchScalarField.C
@@ -75,7 +75,7 @@ tmp<scalarField> mutkRoughWallFunctionFvPatchScalarField::calcMut() const
     const volScalarField& k = tk();
     const scalarField& muw = rasModel.mu().boundaryField()[patchI];
 
-    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu25 = pow025(Cmu_);
 
     tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
     scalarField& mutw = tmutw();
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutkWallFunction/mutkWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutkWallFunction/mutkWallFunctionFvPatchScalarField.C
index 5ab8324cda..674097af09 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutkWallFunction/mutkWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutkWallFunction/mutkWallFunctionFvPatchScalarField.C
@@ -83,7 +83,7 @@ tmp<scalarField> mutkWallFunctionFvPatchScalarField::calcMut() const
     const volScalarField& k = tk();
     const scalarField& muw = rasModel.mu().boundaryField()[patchI];
 
-    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu25 = pow025(Cmu_);
 
     tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
     scalarField& mutw = tmutw();
@@ -215,7 +215,7 @@ tmp<scalarField> mutkWallFunctionFvPatchScalarField::yPlus() const
     const scalarField& muw = rasModel.mu().boundaryField()[patchI];
     const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
 
-    return pow(Cmu_, 0.25)*y*sqrt(kwc)/(muw/rhow);
+    return pow025(Cmu_)*y*sqrt(kwc)/(muw/rhow);
 }
 
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
index 6e475094ec..c60d6628e5 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
@@ -187,7 +187,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
     const scalarField& y = rasModel.y()[patch().index()];
 
-    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu25 = pow025(Cmu_);
 
     volScalarField& G = const_cast<volScalarField&>
         (db().lookupObject<volScalarField>(GName_));
diff --git a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
index b8c1f48ee4..a755d4e56e 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
@@ -346,7 +346,7 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
 
     volVectorField gradK = fvc::grad(k_);
     volVectorField gradOmega = fvc::grad(omega_);
-    volScalarField L = sqrt(k_)/(pow(Cmu_, 0.25)*(omega_ + omegaSmall_));
+    volScalarField L = sqrt(k_)/(pow025(Cmu_)*(omega_ + omegaSmall_));
     volScalarField CDkOmega =
         (2.0*alphaOmega2_)*(gradK & gradOmega)/(omega_ + omegaSmall_);
     volScalarField F1 = this->F1(CDkOmega);
diff --git a/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C b/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C
index c2546e659e..3e8412689b 100644
--- a/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C
+++ b/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C
@@ -171,7 +171,7 @@ tmp<scalarField> RASModel::yPlus(const label patchNo, const scalar Cmu) const
 
     if (isType<wallFvPatch>(curPatch))
     {
-        Yp = pow(Cmu, 0.25)
+        Yp = pow025(Cmu)
             *y_[patchNo]
             *sqrt(k()().boundaryField()[patchNo].patchInternalField())
             /nu().boundaryField()[patchNo];
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
index dac5ccb0f1..179ea640c5 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
@@ -118,7 +118,7 @@ void turbulentMixingLengthFrequencyInletFvPatchScalarField::updateCoeffs()
     const scalar Cmu =
         rasModel.coeffDict().lookupOrDefault<scalar>("Cmu", 0.09);
 
-    const scalar Cmu25 = pow(Cmu, 0.25);
+    const scalar Cmu25 = pow025(Cmu);
 
     const fvPatchField<scalar>& kp =
         patch().lookupPatchField<volScalarField, scalar>(kName_);
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
index 5f771ec663..cdd6582d7a 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
@@ -185,7 +185,7 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
     const scalarField& y = rasModel.y()[patchI];
 
-    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu25 = pow025(Cmu_);
     const scalar Cmu75 = pow(Cmu_, 0.75);
 
     volScalarField& G =
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C
index f81cb76c5c..fdd5428a71 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C
@@ -74,7 +74,7 @@ tmp<scalarField> nutkRoughWallFunctionFvPatchScalarField::calcNut() const
     const volScalarField& k = tk();
     const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
 
-    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu25 = pow025(Cmu_);
 
     tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
     scalarField& nutw = tnutw();
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C
index dbb07712d8..6f786d8f5c 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C
@@ -83,7 +83,7 @@ tmp<scalarField> nutkWallFunctionFvPatchScalarField::calcNut() const
     const volScalarField& k = tk();
     const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
 
-    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu25 = pow025(Cmu_);
 
     tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
     scalarField& nutw = tnutw();
@@ -223,7 +223,7 @@ tmp<scalarField> nutkWallFunctionFvPatchScalarField::yPlus() const
     const scalarField kwc = k.boundaryField()[patchI].patchInternalField();
     const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
 
-    return pow(Cmu_, 0.25)*y*sqrt(kwc)/nuw;
+    return pow025(Cmu_)*y*sqrt(kwc)/nuw;
 }
 
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
index c39b20bf0c..b494541846 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
@@ -186,7 +186,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
     const scalarField& y = rasModel.y()[patchI];
 
-    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu25 = pow025(Cmu_);
 
     volScalarField& G =
         const_cast<volScalarField&>(db().lookupObject<volScalarField>(GName_));
diff --git a/src/turbulenceModels/incompressible/RAS/include/nonLinearWallFunctionsI.H b/src/turbulenceModels/incompressible/RAS/include/nonLinearWallFunctionsI.H
index e88e0f8ee3..8b10335481 100644
--- a/src/turbulenceModels/incompressible/RAS/include/nonLinearWallFunctionsI.H
+++ b/src/turbulenceModels/incompressible/RAS/include/nonLinearWallFunctionsI.H
@@ -75,7 +75,7 @@ Description
                 label faceCelli = curPatch.faceCells()[facei];
 
                 //- using local Cmu !
-                scalar Cmu25 = pow(Cmu_[faceCelli], 0.25);
+                scalar Cmu25 = pow025(Cmu_[faceCelli]);
                 scalar Cmu75 = pow(Cmu_[faceCelli], 0.75);
 
                 scalar yPlus =
diff --git a/src/turbulenceModels/incompressible/RAS/include/wallNonlinearViscosityI.H b/src/turbulenceModels/incompressible/RAS/include/wallNonlinearViscosityI.H
index b24b23c282..5ff4e58a00 100644
--- a/src/turbulenceModels/incompressible/RAS/include/wallNonlinearViscosityI.H
+++ b/src/turbulenceModels/incompressible/RAS/include/wallNonlinearViscosityI.H
@@ -49,7 +49,7 @@ Description
                 label faceCelli = curPatch.faceCells()[facei];
 
                 //- Using local Cmu
-                scalar Cmu25 = pow(Cmu_[faceCelli], 0.25);
+                scalar Cmu25 = pow025(Cmu_[faceCelli]);
 
                 scalar yPlus =
                     Cmu25*y_[patchi][facei]*sqrt(k_[faceCelli])/nuw[facei];