diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/Make/options b/applications/solvers/compressible/rhoPorousSimpleFoam/Make/options
index 61456e486af0422ae88af5d5b11203a507794c01..4f198992b905351e4a49dcb054dd964e24008622 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermalPorousZone/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
@@ -8,6 +9,7 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lbasicThermophysicalModels \
+    -lthermalPorousZone \
     -lspecie \
     -lcompressibleRASModels \
     -lfiniteVolume \
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
index 1177cba2a4b86adcdc0081e5305b9095214976a9..09b75191db13da90f31406132a949d4d81157ad8 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
@@ -64,7 +64,7 @@
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
-    porousZones pZones(mesh);
+    thermalPorousZones pZones(mesh);
     Switch pressureImplicitPorosity(false);
 
     int nUCorr = 0;
@@ -84,4 +84,3 @@
             pressureImplicitPorosity = true;
         }
     }
-
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
index 605b8820d1816daeaa6a6b4a2f20965e388a1c8f..10e84962063c7d96dc5365318f1fc60c9b808c7b 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
@@ -9,6 +9,8 @@
       - p*fvc::div(phi/fvc::interpolate(rho))
     );
 
+    pZones.addEnthalpySource(thermo, hEqn);
+
     hEqn.relax();
 
     eqnResidual = hEqn.solve().initialResidual();
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C b/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
index 10bfeb61db53b4769d3cfd5b006a12ba587fa4fb..5b620be33502eaecbede0e123b94f0404b4a9dae 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
@@ -34,7 +34,7 @@ Description
 #include "fvCFD.H"
 #include "basicPsiThermo.H"
 #include "RASModel.H"
-#include "porousZones.H"
+#include "thermalPorousZones.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/dimensionSet/dimensionSet.C b/src/OpenFOAM/dimensionSet/dimensionSet.C
index 9fba472d9801d92233fefc7e7fe941cf7d6e3f69..e60a3998d449859a078c6a750e98d45181c0d999 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 e61d3ec4df03d416bec6e57905516119aae3caba..13b0a899a9c87507d022fd5745f3c1cbad711bdc 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 8761396a0ff0cef0c55fe96f587d58a7d11ba03a..054d1eb43b45512ec79cf8a61bb8abf814fad8c6 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 492d5064ef52376b7bf9a559ee352078b0e8d4e9..c6e46168c9ee70283884b2bf2e50689996e9439d 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 ade367e8ff019c78a2f0f128d9e1c6768df99377..4aad204f790f2367394cf093d54fe75f043c35f0 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 1b5e5d22bffd475fbf79e8825e8639077532f4c3..e2d565283e2929771cee3d868fcf18d15d03429e 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 fe9067a8f57dab085599253b534afb4cf6a7d68d..8a308adbb58238d64cefcdeee479866c4d68d6c8 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 45f1e86ea01e95a8f36496075a1d09b4595d9184..019b7319797a47c54519a86505f32fc8ce2954ef 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 4f2eccfa9ceea08bab64a6720d31c2161cf1dd6d..ed8e2da5e3ded05a93ad8d0afe011fd91a7f04b2 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 a1d0865d205664100e4826ace2cbf57aab3792b3..bd25d64a426c4f407399e5e4f77452b29685895d 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 fa6a2a58d0a2874f56b6ed7db2e7bc15af060f4c..3aaf5ddf3b48a6fabc522eb9345237bedaf5564e 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 4a2d9166b4c3bbaf0d0e028974dd11d333639e66..3b50684200b1b16d37761662e8a9502fd0c3b5dc 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 cf5d9cb31dd694a31d3b30845fc71660ff58d35b..69961ab2cc3ae4f4886729172b14269e0ecea0ed 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/finiteVolume/cfdTools/general/porousMedia/PorousZones.C b/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.C
new file mode 100644
index 0000000000000000000000000000000000000000..444cd2513107f330b2d5a071462f9e3571dda3ed
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.C
@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "PorousZones.H"
+#include "Time.H"
+#include "volFields.H"
+#include "fvm.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class ZoneType>
+template<class Type>
+void Foam::PorousZones<ZoneType>::modifyDdt(fvMatrix<Type>& m) const
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).modifyDdt(m);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ZoneType>
+Foam::PorousZones<ZoneType>::PorousZones
+(
+    const fvMesh& mesh
+)
+:
+    IOPtrList<ZoneType>
+    (
+        IOobject
+        (
+            "porousZones",
+            mesh.time().constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        typename ZoneType::iNew(mesh)
+    ),
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ZoneType>
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> >
+Foam::PorousZones<ZoneType>::ddt
+(
+    GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type> > tres = fvm::ddt(vf);
+    modifyDdt(tres());
+    return tres;
+}
+
+
+template<class ZoneType>
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> >
+Foam::PorousZones<ZoneType>::ddt
+(
+    const oneField&,
+    GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type> > tres = fvm::ddt(vf);
+    modifyDdt(tres());
+    return tres;
+}
+
+
+template<class ZoneType>
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> >
+Foam::PorousZones<ZoneType>::ddt
+(
+    const dimensionedScalar& rho,
+    GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type> > tres = fvm::ddt(rho,vf);
+    modifyDdt(tres());
+    return tres;
+}
+
+
+template<class ZoneType>
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type> >
+Foam::PorousZones<ZoneType>::ddt
+(
+    const volScalarField& rho,
+    GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type> > tres = fvm::ddt(rho,vf);
+    modifyDdt(tres());
+    return tres;
+}
+
+template<class ZoneType>
+void Foam::PorousZones<ZoneType>::addResistance(fvVectorMatrix& UEqn) const
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).addResistance(UEqn);
+    }
+}
+
+
+template<class ZoneType>
+void Foam::PorousZones<ZoneType>::addResistance
+(
+    const fvVectorMatrix& UEqn,
+    volTensorField& AU
+) const
+{
+    // addResistance for each zone, delaying the correction of the
+    // precessor BCs of AU
+    forAll(*this, i)
+    {
+        this->operator[](i).addResistance(UEqn, AU, false);
+    }
+
+    // Correct the boundary conditions of the tensorial diagonal to ensure
+    // processor bounaries are correctly handled when AU^-1 is interpolated
+    // for the pressure equation.
+    AU.correctBoundaryConditions();
+}
+
+
+template<class ZoneType>
+bool Foam::PorousZones<ZoneType>::readData(Istream& is)
+{
+    this->clear();
+
+    IOPtrList<ZoneType> newLst
+    (
+        IOobject
+        (
+            "porousZones",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false     // Don't re-register new zones with objectRegistry
+        ),
+        typename ZoneType::iNew(mesh_)
+    );
+
+    transfer(newLst);
+
+    return is.good();
+}
+
+
+template<class ZoneType>
+bool Foam::PorousZones<ZoneType>::writeData(Ostream& os, bool subDict) const
+{
+    // Write size of list
+    os << nl << this->size();
+
+    // Write beginning of contents
+    os << nl << token::BEGIN_LIST;
+
+    // Write list contents
+    forAll(*this, i)
+    {
+        os << nl;
+        this->operator[](i).writeDict(os, subDict);
+    }
+
+    // Write end of contents
+    os << token::END_LIST << nl;
+
+    // Check state of IOstream
+    return os.good();
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.H b/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.H
new file mode 100644
index 0000000000000000000000000000000000000000..cc55fc42646b30589cf4f32da1b40d16c3461ca4
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/porousMedia/PorousZones.H
@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::PorousZones<ZoneType>
+
+Description
+    A centralized ZoneType collection.
+
+    Container class for a set of ZoneType with the ZoneType member
+    functions implemented to loop over the functions for each ZoneType.
+
+SourceFiles
+    PorousZones.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PorousZones_H
+#define PorousZones_H
+
+#include "IOPtrList.H"
+
+#include "volFieldsFwd.H"
+#include "fvMatricesFwd.H"
+#include "dimensionedScalarFwd.H"
+#include "oneField.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+class fvMesh;
+
+
+/*---------------------------------------------------------------------------*\
+                           Class PorousZones Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ZoneType>
+class PorousZones
+:
+    public IOPtrList<ZoneType>
+{
+    // Private data
+
+        //- Reference to the finite volume mesh this zone is part of
+        const fvMesh& mesh_;
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        PorousZones(const PorousZones<ZoneType>&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const PorousZones<ZoneType>&);
+
+
+        //- modify time derivative elements
+        template<class Type>
+        void modifyDdt(fvMatrix<Type>&) const;
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh
+        //  with automatically constructed coordinate systems list
+        PorousZones(const fvMesh&);
+
+
+    // Member Functions
+
+        //- mirror fvm::ddt with porosity
+        template<class Type>
+        tmp<fvMatrix<Type> > ddt
+        (
+            GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- mirror fvm::ddt with porosity
+        template<class Type>
+        tmp<fvMatrix<Type> > ddt
+        (
+            const oneField&,
+            GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- mirror fvm::ddt with porosity
+        template<class Type>
+        tmp<fvMatrix<Type> > ddt
+        (
+            const dimensionedScalar&,
+            GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- mirror fvm::ddt with porosity
+        template<class Type>
+        tmp<fvMatrix<Type> > ddt
+        (
+            const volScalarField&,
+            GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- Add the viscous and inertial resistance force contribution
+        //  to the momentum equation
+        void addResistance(fvVectorMatrix& UEqn) const;
+
+        //- Add the viscous and inertial resistance force contribution
+        //  to the tensorial diagonal
+        void addResistance
+        (
+            const fvVectorMatrix& UEqn,
+            volTensorField& AU
+        ) const;
+
+        //- read modified data
+        virtual bool readData(Istream&);
+
+        //- write data
+        bool writeData(Ostream&, bool subDict = true) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "PorousZones.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H
index 5b4d16af5f6ebe4e498db526856656cc0c8ab2ae..9f35446e4dae33e7ef85765da7d1671ee3996e1c 100644
--- a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H
+++ b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H
@@ -198,7 +198,6 @@ public:
             return autoPtr<porousZone>(NULL);
         }
 
-
         //- Return pointer to new porousZone created on freestore from Istream
         class iNew
         {
@@ -222,6 +221,11 @@ public:
         };
 
 
+    //- Destructor
+    virtual ~porousZone()
+    {}
+
+
     // Member Functions
 
         // Access
@@ -232,6 +236,12 @@ public:
                 return name_;
             }
 
+            //- Return mesh
+            const fvMesh& mesh() const
+            {
+                return mesh_;
+            }
+
             //- cellZone number
             label zoneId() const
             {
@@ -275,7 +285,7 @@ public:
             }
 
 
-        //- modify time derivative elements according to porosity
+        //- Modify time derivative elements according to porosity
         template<class Type>
         void modifyDdt(fvMatrix<Type>&) const;
 
@@ -294,7 +304,7 @@ public:
         ) const;
 
         //- Write the porousZone dictionary
-        void writeDict(Ostream&, bool subDict = true) const;
+        virtual void writeDict(Ostream&, bool subDict = true) const;
 
 
     // Ostream Operator
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.C b/src/finiteVolume/cfdTools/general/porousMedia/porousZones.C
index 596398e6d620636f8ae676cb3efeebba137568f3..ffcb1ef2e8cc6637e2cc80700b8c0f095820d581 100644
--- a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.C
+++ b/src/finiteVolume/cfdTools/general/porousMedia/porousZones.C
@@ -25,8 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "porousZones.H"
-#include "Time.H"
-#include "volFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -35,105 +33,4 @@ namespace Foam
     defineTemplateTypeNameAndDebug(IOPtrList<porousZone>, 0);
 }
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::porousZones::porousZones
-(
-    const fvMesh& mesh
-)
-:
-    IOPtrList<porousZone>
-    (
-        IOobject
-        (
-            "porousZones",
-            mesh.time().constant(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
-        porousZone::iNew(mesh)
-    ),
-    mesh_(mesh)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::porousZones::addResistance(fvVectorMatrix& UEqn) const
-{
-    forAll(*this, i)
-    {
-        operator[](i).addResistance(UEqn);
-    }
-}
-
-
-void Foam::porousZones::addResistance
-(
-    const fvVectorMatrix& UEqn,
-    volTensorField& AU
-) const
-{
-    // addResistance for each zone, delaying the correction of the
-    // precessor BCs of AU
-    forAll(*this, i)
-    {
-        operator[](i).addResistance(UEqn, AU, false);
-    }
-
-    // Correct the boundary conditions of the tensorial diagonal to ensure
-    // processor bounaries are correctly handled when AU^-1 is interpolated
-    // for the pressure equation.
-    AU.correctBoundaryConditions();
-}
-
-
-bool Foam::porousZones::readData(Istream& is)
-{
-    clear();
-
-    IOPtrList<porousZone> newLst
-    (
-        IOobject
-        (
-            "porousZones",
-            mesh_.time().constant(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false     // Don't re-register new zones with objectRegistry
-        ),
-        porousZone::iNew(mesh_)
-    );
-
-    transfer(newLst);
-
-    return is.good();
-}
-
-
-bool Foam::porousZones::writeData(Ostream& os, bool subDict) const
-{
-    // Write size of list
-    os << nl << size();
-
-    // Write beginning of contents
-    os << nl << token::BEGIN_LIST;
-
-    // Write list contents
-    forAll(*this, i)
-    {
-        os << nl;
-        operator[](i).writeDict(os, subDict);
-    }
-
-    // Write end of contents
-    os << token::END_LIST << nl;
-
-    // Check state of IOstream
-    return os.good();
-}
-
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H b/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H
index 0d9dd2bd9906eb62f7f6e582bad92eba08dea910..abc0ea633d81f4a65d9bd208f684e8a71486f896 100644
--- a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H
+++ b/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H
@@ -22,152 +22,23 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Class
+Typedef
     Foam::porousZones
 
-Description
-    A centralized porousZone collection.
-
-    Container class for a set of porousZones with the porousZone member
-    functions implemented to loop over the functions for each porousZone.
-
-    The input file @c constant/porousZone is implemented as
-    IOPtrList\<porousZone\> and contains the following type of data:
-
-    @verbatim
-    1
-    (
-    cat1
-    {
-        coordinateSystem    system_10;
-        porosity    0.781;
-        Darcy
-        {
-            d   d [0 -2 0 0 0]  (-1000 -1000 0.50753e+08);
-            f   f [0 -1 0 0 0]  (-1000 -1000 12.83);
-        }
-    }
-    )
-    @endverbatim
-
-SourceFiles
-    porousZones.C
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef porousZones_H
 #define porousZones_H
 
+#include "PorousZones.H"
 #include "porousZone.H"
-#include "IOPtrList.H"
-
-#include "volFieldsFwd.H"
-#include "fvMatrix.H"
-#include "oneField.H"
-
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-/*---------------------------------------------------------------------------*\
-                           Class porousZones Declaration
-\*---------------------------------------------------------------------------*/
-
-class porousZones
-:
-    public IOPtrList<porousZone>
-{
-    // Private data
-
-        //- Reference to the finite volume mesh this zone is part of
-        const fvMesh& mesh_;
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        porousZones(const porousZones&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const porousZones&);
-
-
-        //- modify time derivative elements
-        template<class Type>
-        void modifyDdt(fvMatrix<Type>&) const;
-
-public:
-
-    // Constructors
-
-        //- Construct from fvMesh
-        //  with automatically constructed coordinate systems list
-        porousZones(const fvMesh&);
-
-
-    // Member Functions
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            const oneField&,
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- mirror fvm::ddt with porosity
-        template<class Type>
-        tmp<fvMatrix<Type> > ddt
-        (
-            const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
-        );
-
-        //- Add the viscous and inertial resistance force contribution
-        //  to the momentum equation
-        void addResistance(fvVectorMatrix& UEqn) const;
-
-        //- Add the viscous and inertial resistance force contribution
-        //  to the tensorial diagonal
-        void addResistance
-        (
-            const fvVectorMatrix& UEqn,
-            volTensorField& AU
-        ) const;
-
-        //- read modified data
-        virtual bool readData(Istream&);
-
-        //- write data
-        bool writeData(Ostream&, bool subDict = true) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "porousZonesTemplates.C"
-#endif
+    typedef PorousZones<porousZone> porousZones;
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/basic/Particle/ParticleI.H b/src/lagrangian/basic/Particle/ParticleI.H
index 65ac64e5251c25a460e79eb70ae5aeadc3d523bb..2e6d8f868f0ec8526f457bbe75b28a376dec1c5d 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/intermediate/particleForces/particleForces.H b/src/lagrangian/intermediate/particleForces/particleForces.H
index fbe88a6dc832e7fd545a1a673aaf2dce90178a23..f0be486381613268ee513fca876b3b2581ecd30f 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/thermophysicalModels/Allwmake b/src/thermophysicalModels/Allwmake
index 7fbee69034eae6d4ec0cff99068f1a44b1cfee24..38f4d54b4767c2475ea80c205aa84b824ae0c098 100755
--- a/src/thermophysicalModels/Allwmake
+++ b/src/thermophysicalModels/Allwmake
@@ -16,5 +16,6 @@ wmake libso chemistryModel
 wmake libso pdfs
 wmake libso radiation
 wmake libso barotropicCompressibilityModel
+wmake libso thermalPorousZone
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/src/thermophysicalModels/thermalPorousZone/Make/files b/src/thermophysicalModels/thermalPorousZone/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..046355feaf96fe50a6283172613cf0cf9cac070f
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/Make/files
@@ -0,0 +1,4 @@
+thermalPorousZone/thermalPorousZone.C
+thermalPorousZone/thermalPorousZones.C
+
+LIB = $(FOAM_LIBBIN)/libthermalPorousZone
diff --git a/src/thermophysicalModels/thermalPorousZone/Make/options b/src/thermophysicalModels/thermalPorousZone/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..ca2d6183be35fd2d39863ff61a2c74c01666ee4c
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/Make/options
@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
+
+LIB_LIBS = \
+    -lbasicThermophysicalModels \
+    -lmeshTools \
+    -lfiniteVolume
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C
new file mode 100644
index 0000000000000000000000000000000000000000..699d029e0836026222f8d348439943f4c0034756
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C
@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "thermalPorousZone.H"
+#include "fvMesh.H"
+#include "fvMatrices.H"
+#include "basicThermo.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::thermalPorousZone::thermalPorousZone
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    porousZone(name, mesh, dict),
+    T_("T", dimTemperature, -GREAT)
+{
+    if (const dictionary* dictPtr = dict.subDictPtr("thermalModel"))
+    {
+        word thermalModel(dictPtr->lookup("type"));
+
+        if (thermalModel == "fixedTemperature")
+        {
+            dictPtr->lookup("T") >> T_;
+        }
+        else
+        {
+            FatalIOErrorIn
+            (
+                "thermalPorousZone::thermalPorousZone"
+                "("
+                "const word& name, "
+                "const fvMesh& mesh, "
+                "const dictionary& dict"
+                ")",
+                *dictPtr
+            )   << "thermalModel " << thermalModel << " is not supported" << nl
+                << "    Supported thermalModels are: fixedTemperature"
+                << exit(FatalIOError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::thermalPorousZone::addEnthalpySource
+(
+    const basicThermo& thermo,
+    fvScalarMatrix& hEqn
+) const
+{
+    if (zoneId() == -1 || T_.value() < 0.0)
+    {
+        return;
+    }
+
+    const labelList& cells = mesh().cellZones()[zoneId()];
+    const scalarField& V = mesh().V();
+    scalarField& hDiag = hEqn.diag();
+    scalarField& hSource = hEqn.source();
+    const scalarField& rho = thermo.rho();
+
+    scalarField hZone = thermo.h(scalarField(cells.size(), T_.value()), cells);
+    scalar rate = 1e6;
+
+    forAll (cells, i)
+    {
+        hDiag[cells[i]] += rate*V[cells[i]]*rho[cells[i]];
+        hSource[cells[i]] += rate*V[cells[i]]*rho[cells[i]]*hZone[i];
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H
new file mode 100644
index 0000000000000000000000000000000000000000..ab4c982818af637a65d2b720302660cca880d581
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H
@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::thermalPorousZone
+
+Description
+    Porous zone definition based on cell zones including terms for energy
+    equations.
+
+See Also
+    porousZone, thermalPorousZones and coordinateSystems
+
+SourceFiles
+    thermalPorousZone.C
+    thermalPorousZoneTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef thermalPorousZone_H
+#define thermalPorousZone_H
+
+#include "porousZone.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class fvMesh;
+class basicThermo;
+
+/*---------------------------------------------------------------------------*\
+                        Class thermalPorousZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class thermalPorousZone
+:
+    public porousZone
+{
+    // Private data
+
+        //- Temperature in the porous-zone
+        dimensionedScalar T_;
+
+
+        //- Disallow default bitwise copy construct
+        thermalPorousZone(const thermalPorousZone&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const thermalPorousZone&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        thermalPorousZone(const word& name, const fvMesh&, const dictionary&);
+
+        //- Return clone
+        autoPtr<thermalPorousZone> clone() const
+        {
+            notImplemented("autoPtr<thermalPorousZone> clone() const");
+            return autoPtr<thermalPorousZone>(NULL);
+        }
+
+        //- Return pointer to new thermalPorousZone
+        //  created on freestore from Istream
+        class iNew
+        {
+            //- Reference to the finite volume mesh this zone is part of
+            const fvMesh& mesh_;
+
+        public:
+
+            iNew(const fvMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<thermalPorousZone> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+
+                return autoPtr<thermalPorousZone>
+                (
+                    new thermalPorousZone(name, mesh_, dict)
+                );
+            }
+        };
+
+
+    //- Destructor
+    virtual ~thermalPorousZone()
+    {}
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the temperature in the porous-zone
+            const dimensionedScalar& T() const
+            {
+                return T_;
+            }
+
+            //- Edit access to the temperature in the porous-zone
+            dimensionedScalar& T()
+            {
+                return T_;
+            }
+
+        //- Add the thermal source to the enthalpy equation
+        void addEnthalpySource
+        (
+            const basicThermo& thermo,
+            fvScalarMatrix& hEqn
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+//#   include "thermalPorousZoneTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZonesTemplates.C b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZoneTemplates.C
similarity index 52%
rename from src/finiteVolume/cfdTools/general/porousMedia/porousZonesTemplates.C
rename to src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZoneTemplates.C
index f2b7f4b2ad5ea3ece8c2c2aaa1cb930d0479290d..9b2e734c549fb9b82f9afe80dd1d0eea36360828 100644
--- a/src/finiteVolume/cfdTools/general/porousMedia/porousZonesTemplates.C
+++ b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZoneTemplates.C
@@ -24,77 +24,57 @@ License
 
 \*----------------------------------------------------------------------------*/
 
-#include "porousZones.H"
-#include "volFields.H"
-#include "fvMatrix.H"
-#include "fvm.H"
+#include "porousZone.H"
+#include "fvMesh.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class Type>
-void Foam::porousZones::modifyDdt(fvMatrix<Type>& m) const
+template<class RhoFieldType>
+void Foam::porousZone::addPowerLawResistance
+(
+    scalarField& Udiag,
+    const labelList& cells,
+    const scalarField& V,
+    const RhoFieldType& rho,
+    const vectorField& U
+) const
 {
-    forAll(*this, i)
+    const scalar C0 = C0_;
+    const scalar C1m1b2 = (C1_ - 1.0)/2.0;
+
+    forAll (cells, i)
     {
-        operator[](i).modifyDdt(m);
+        Udiag[cells[i]] +=
+            V[cells[i]]*rho[cells[i]]*C0*pow(magSqr(U[cells[i]]), C1m1b2);
     }
 }
 
 
-// * * * * * * * * * * * * * * *  Member Functions * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type> >
-Foam::porousZones::ddt
+template<class RhoFieldType>
+void Foam::porousZone::addViscousInertialResistance
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
-)
+    scalarField& Udiag,
+    vectorField& Usource,
+    const labelList& cells,
+    const scalarField& V,
+    const RhoFieldType& rho,
+    const scalarField& mu,
+    const vectorField& U
+) const
 {
-    tmp<fvMatrix<Type> > tres = fvm::ddt(vf);
-    modifyDdt(tres());
-    return tres;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type> >
-Foam::porousZones::ddt
-(
-    const oneField&,
-    GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    tmp<fvMatrix<Type> > tres = fvm::ddt(vf);
-    modifyDdt(tres());
-    return tres;
-}
+    const tensor& D = D_.value();
+    const tensor& F = F_.value();
 
+    forAll (cells, i)
+    {
+        tensor dragCoeff = mu[cells[i]]*D + (rho[cells[i]]*mag(U[cells[i]]))*F;
+        scalar isoDragCoeff = tr(dragCoeff);
 
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type> >
-Foam::porousZones::ddt
-(
-    const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    tmp<fvMatrix<Type> > tres = fvm::ddt(rho,vf);
-    modifyDdt(tres());
-    return tres;
+        Udiag[cells[i]] += V[cells[i]]*isoDragCoeff;
+        Usource[cells[i]] -=
+            V[cells[i]]*((dragCoeff - I*isoDragCoeff) & U[cells[i]]);
+    }
 }
 
 
-template<class Type>
-Foam::tmp<Foam::fvMatrix<Type> >
-Foam::porousZones::ddt
-(
-    const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    tmp<fvMatrix<Type> > tres = fvm::ddt(rho,vf);
-    modifyDdt(tres());
-    return tres;
-}
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C
new file mode 100644
index 0000000000000000000000000000000000000000..640aebc618e7e5f5c635a1b277dde3c10a2b142d
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C
@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "thermalPorousZones.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug(IOPtrList<thermalPorousZone>, 0);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::thermalPorousZones::thermalPorousZones
+(
+    const fvMesh& mesh
+)
+:
+    PorousZones<thermalPorousZone>(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::thermalPorousZones::addEnthalpySource
+(
+    const basicThermo& thermo,
+    fvScalarMatrix& hEqn
+) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).addEnthalpySource(thermo, hEqn);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H
new file mode 100644
index 0000000000000000000000000000000000000000..e53698b3bf8923ebcff47cc1167a66e42c5f9fba
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Class
+    Foam::thermalPorousZones
+
+Description
+    A centralized thermalPorousZone collection.
+
+    Container class for a set of thermalPorousZones with the thermalPorousZone
+    member functions implemented to loop over the functions for each
+    thermalPorousZone.
+
+    The input file @c constant/thermalPorousZone is implemented as
+    IOPtrList\<thermalPorousZone\> and contains the following type of data:
+
+    @verbatim
+    1
+    (
+    cat1
+    {
+        coordinateSystem    system_10;
+        porosity    0.781;
+        Darcy
+        {
+            d   d [0 -2 0 0 0]  (-1000 -1000 0.50753e+08);
+            f   f [0 -1 0 0 0]  (-1000 -1000 12.83);
+        }
+        Temperature [0 0 1 0 0] 600;
+    }
+    )
+    @endverbatim
+
+SourceFiles
+    thermalPorousZones.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef thermalPorousZones_H
+#define thermalPorousZones_H
+
+#include "PorousZones.H"
+#include "thermalPorousZone.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class thermalPorousZones Declaration
+\*---------------------------------------------------------------------------*/
+
+class thermalPorousZones
+:
+    public PorousZones<thermalPorousZone>
+{
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh
+        thermalPorousZones(const fvMesh&);
+
+
+    // Member Functions
+
+        //- Add the thermal source to the enthalpy equation
+        void addEnthalpySource
+        (
+            const basicThermo& thermo,
+            fvScalarMatrix& hEqn
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C
index 1d4054fb04fc37f1cdd8c8bcc7aae8acd3512cf8..f0b77bb1f26666be385bd073550a775db0bb317d 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 5893f2072ebff0ddc2a2fb4476d45ba8897029de..bacfe9615e079b000086aa87f5e49fee5fc66cad 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 e8fdf3edd5eb70e6e8e03646de0d51d0c536d422..b5022e2756469fcb0b70eea08ae1be15633cb0fd 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 07382bad5236b4d841ae82a7a88d33cff391931d..a129ec8bffbdcade0afde63a5015b02b804589c7 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 6af9687f7b8563407a3f7b8391fc7fca487d6631..67353e43198075d5c553e33f66fbd1ca96d2fd84 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 5ab8324cda85beb987ea456744ebee6ae8e4356c..674097af098bdd9929aaab63560d53f8ad343bd9 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 6e475094ec7cfad5ad5c401508c3a65b78426ab0..c60d6628e5efea62d1e5c94cb821222104a65586 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 b8c1f48ee4ab853fb65603df611d3ff7b4e5109e..a755d4e56e8a02b6fccea56adf22b75043471c6a 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 c2546e659eba7a99df77c333eab929d5bae50096..3e8412689bdeab0af918d859641d79d06070368f 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 dac5ccb0f16498e57b9541d9149d6c0310a48104..179ea640c50354a9608f01c9d0f698d9a11b826f 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 5f771ec663cedf7d7278d016af67bf1e712d3a0e..cdd6582d7a911d6cc7b8196e54fb0bb6cded8b5d 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 f81cb76c5c5d8c95f9cd59f4eab862e179e00b47..fdd5428a711884cd37865034e77ba51c266adfd2 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 dbb07712d8b0e0ee520035ed92bd67da8699ee50..6f786d8f5cd58945680ac6bb909a107d8f124718 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 c39b20bf0c68ce37346a1bae3551eb843b34de8c..b49454184672c7a8343059053502e4d5b743305c 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 e88e0f8ee317b0bcff4337a0ad7ccf83d8da4ef1..8b103354816b2ce5d3814f33042376a8c81e9e65 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 b24b23c2820c9f7d61ad7a4f02455f3282d84892..5ff4e58a0054b7ebb63e509434d143c41b659371 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];
diff --git a/tutorials/compressible/rhoPimpleFoam/angledDuct/constant/polyMesh/boundary b/tutorials/compressible/rhoPimpleFoam/angledDuct/constant/polyMesh/boundary
index d060c920afda73a3022c1c66814682862962ac9d..0abd1608aba0dcb6aa66c9488133a3c4b51c7588 100644
--- a/tutorials/compressible/rhoPimpleFoam/angledDuct/constant/polyMesh/boundary
+++ b/tutorials/compressible/rhoPimpleFoam/angledDuct/constant/polyMesh/boundary
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
index 5adb07e26c13092828a40377945ec2f4a98ba0d0..0abd1608aba0dcb6aa66c9488133a3c4b51c7588 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/polyMesh/boundary
@@ -1,14 +1,14 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
     version     2.0;
-    format      binary;
+    format      ascii;
     class       polyBoundaryMesh;
     location    "constant/polyMesh";
     object      boundary;
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/porousZones b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/porousZones
index 634799837eaf7009528df07d027deafc7c6cc1f4..43602a72a6b0fed8c3528eb0791f130bcab38fcc 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/porousZones
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/porousZones
@@ -30,6 +30,13 @@ FoamFile
             d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
             f   f [0 -1 0 0 0 0 0] (0 0 0);
         }
+
+        thermalModel
+        {
+            type fixedTemperature;
+
+            T  T [0 0 0 1 0] 350;
+        }
     }
 )
 
diff --git a/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCubeArgon/constant/polyMesh/boundary b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCubeArgon/constant/polyMesh/boundary
index 0eb90da00f8ac5dfe6fc2600c7eef27018a5b5ce..850c9baa1ec078697b17ec724c99578c6047e10d 100644
--- a/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCubeArgon/constant/polyMesh/boundary
+++ b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCubeArgon/constant/polyMesh/boundary
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
diff --git a/tutorials/incompressible/pimpleDyMFoam/movingCone/system/tetFemSolution b/tutorials/incompressible/pimpleDyMFoam/movingCone/system/tetFemSolution
deleted file mode 100644
index 7e16a5b03d6529c4af7fc10422797b31a431d842..0000000000000000000000000000000000000000
--- a/tutorials/incompressible/pimpleDyMFoam/movingCone/system/tetFemSolution
+++ /dev/null
@@ -1,24 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      binary;
-    class       dictionary;
-    location    "system";
-    object      tetFemSolution;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-solvers
-{
-    motionU         ICCG 1e-06 0;
-}
-
-
-// ************************************************************************* //