diff --git a/src/thermophysicalModels/basic/Make/files b/src/thermophysicalModels/basic/Make/files
index e88b8bf66075109dc0418c3946c347a34d606f4f..4cea8bf57d4c6f7e179760a661e095f4a862044c 100644
--- a/src/thermophysicalModels/basic/Make/files
+++ b/src/thermophysicalModels/basic/Make/files
@@ -2,14 +2,21 @@ basicMixture = mixtures/basicMixture
 basicThermo = basicThermo
 
 $(basicMixture)/basicMixture.C
+$(basicMixture)/basicMixtures.C
 $(basicThermo)/basicThermo.C
 $(basicThermo)/newBasicThermo.C
-$(basicThermo)/basicThermos.C
+
+hThermo/hThermos.C
+eThermo/eThermos.C
 
 derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C
 derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C
 derivedFvPatchFields/mixedEnthalpy/mixedEnthalpyFvPatchScalarField.C
 
+derivedFvPatchFields/fixedInternalEnergy/fixedInternalEnergyFvPatchScalarField.C
+derivedFvPatchFields/gradientInternalEnergy/gradientInternalEnergyFvPatchScalarField.C
+derivedFvPatchFields/mixedInternalEnergy/mixedInternalEnergyFvPatchScalarField.C
+
 derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C
 
 LIB = $(FOAM_LIBBIN)/libbasicThermophysicalModels
diff --git a/src/thermophysicalModels/basic/basicThermo/basicThermo.C b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
index dffda25d6923938ec7fe14fec6d0a92e34436889..e2e92ed4349bd7f121819de7106fba7acd208c55 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.C
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
@@ -31,6 +31,9 @@ License
 #include "fixedEnthalpyFvPatchScalarField.H"
 #include "gradientEnthalpyFvPatchScalarField.H"
 #include "mixedEnthalpyFvPatchScalarField.H"
+#include "fixedInternalEnergyFvPatchScalarField.H"
+#include "gradientInternalEnergyFvPatchScalarField.H"
+#include "mixedInternalEnergyFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -74,9 +77,9 @@ wordList basicThermo::hBoundaryTypes()
     return hbt;
 }
 
-void basicThermo::hBoundaryCorrection(volScalarField& h_)
+void basicThermo::hBoundaryCorrection(volScalarField& h)
 {
-    volScalarField::GeometricBoundaryField& hbf = h_.boundaryField();
+    volScalarField::GeometricBoundaryField& hbf = h.boundaryField();
 
     forAll(hbf, patchi)
     {
@@ -93,6 +96,53 @@ void basicThermo::hBoundaryCorrection(volScalarField& h_)
     }
 }
 
+wordList basicThermo::eBoundaryTypes()
+{
+    const volScalarField::GeometricBoundaryField& tbf = T_.boundaryField();
+
+    wordList ebt = tbf.types();
+
+    forAll(tbf, patchi)
+    {
+        if (isA<fixedValueFvPatchScalarField>(tbf[patchi]))
+        {
+            ebt[patchi] = fixedInternalEnergyFvPatchScalarField::typeName;
+        }
+        else if
+        (
+            isA<zeroGradientFvPatchScalarField>(tbf[patchi])
+         || isA<fixedGradientFvPatchScalarField>(tbf[patchi])
+        )
+        {
+            ebt[patchi] = gradientInternalEnergyFvPatchScalarField::typeName;
+        }
+        else if (isA<mixedFvPatchScalarField>(tbf[patchi]))
+        {
+            ebt[patchi] = mixedInternalEnergyFvPatchScalarField::typeName;
+        }
+    }
+
+    return ebt;
+}
+
+void basicThermo::eBoundaryCorrection(volScalarField& e)
+{
+    volScalarField::GeometricBoundaryField& ebf = e.boundaryField();
+
+    forAll(ebf, patchi)
+    {
+        if (isA<gradientInternalEnergyFvPatchScalarField>(ebf[patchi]))
+        {
+            refCast<gradientInternalEnergyFvPatchScalarField>(ebf[patchi])
+                .gradient() = ebf[patchi].fvPatchField::snGrad();
+        }
+        else if (isA<mixedInternalEnergyFvPatchScalarField>(ebf[patchi]))
+        {
+            refCast<mixedInternalEnergyFvPatchScalarField>(ebf[patchi])
+                .refGrad() = ebf[patchi].fvPatchField::snGrad();
+        }
+    }
+}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
diff --git a/src/thermophysicalModels/basic/basicThermo/basicThermo.H b/src/thermophysicalModels/basic/basicThermo/basicThermo.H
index 30e8aa65a218cf0860335043a96017e4f2006003..da1cf8ea30882d31b7b49e81c057c99e23746c14 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.H
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.H
@@ -74,6 +74,9 @@ protected:
         wordList hBoundaryTypes();
         void hBoundaryCorrection(volScalarField& h);
 
+        wordList eBoundaryTypes();
+        void eBoundaryCorrection(volScalarField& e);
+
         //- Construct as copy (not implemented)
         basicThermo(const basicThermo&);
 
@@ -107,9 +110,8 @@ public:
         static autoPtr<basicThermo> New(const fvMesh&);
 
 
-    // Destructor
-
-        virtual ~basicThermo();
+    //- Destructor
+    virtual ~basicThermo();
 
 
     // Member functions
@@ -122,13 +124,13 @@ public:
 
             //- Pressure [Pa]
             //  Non-const access allowed for transport equations
-            volScalarField& p()
+            virtual volScalarField& p()
             {
                 return p_;
             }
 
             //- Pressure [Pa]
-            const volScalarField& p() const
+            virtual const volScalarField& p() const
             {
                 return p_;
             }
@@ -193,23 +195,52 @@ public:
                 return volScalarField::null();
             }
 
+            //- Internal energy for cell-set [J/kg]
+            virtual tmp<scalarField> e
+            (
+                const scalarField& T,
+                const labelList& cells
+            ) const
+            {
+                notImplemented
+                (
+                    "basicThermo::e"
+                    "(const scalarField& T, const labelList& cells) const"
+                );
+                return tmp<scalarField>(NULL);
+            }
+
+            //-Internal energy for patch [J/kg]
+            virtual tmp<scalarField> e
+            (
+                const scalarField& T,
+                const label patchi
+            ) const
+            {
+                notImplemented
+                (
+                    "basicThermo::e"
+                    "(const scalarField& T, const label patchi) const"
+                );
+                return tmp<scalarField>(NULL);
+            }
 
         // Fields derived from thermodynamic state variables
 
             //- Temperature [K]
-            const volScalarField& T() const
+            virtual const volScalarField& T() const
             {
                 return T_;
             }
 
             //- Density [kg/m^3]
-            tmp<volScalarField> rho() const
+            virtual tmp<volScalarField> rho() const
             {
                 return p_*psi();
             }
 
             //- Compressibility [s^2/m^2]
-            const volScalarField& psi() const
+            virtual const volScalarField& psi() const
             {
                 return psi_;
             }
@@ -236,6 +267,21 @@ public:
                 return volScalarField::null();
             }
 
+            //- Heat capacity at constant volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cv
+            (
+                const scalarField& T,
+                const label patchi
+            ) const
+            {
+                notImplemented
+                (
+                    "basicThermo::Cv"
+                    "(const scalarField& T, const label patchi) const"
+                );
+                return tmp<scalarField>(NULL);
+            }
+
             //- Heat capacity at constant volume [J/kg/K]
             virtual tmp<volScalarField> Cv() const
             {
@@ -247,13 +293,13 @@ public:
         // Access to transport state variables
 
             //- Dynamic viscosity of mixture [kg/ms]
-            const volScalarField& mu() const
+            virtual const volScalarField& mu() const
             {
                 return mu_;
             }
 
             //- Thermal diffusivity for enthalpy of mixture [kg/ms]
-            const volScalarField& alpha() const
+            virtual const volScalarField& alpha() const
             {
                 return alpha_;
             }
diff --git a/src/thermophysicalModels/basic/basicThermo/makeBasicThermo.H b/src/thermophysicalModels/basic/basicThermo/makeBasicThermo.H
index e37bfc7b05130c75c84a132838bd83f5e104748a..13a1bb853c9f82a03fd86b28ccc21a2b9a6ba78b 100644
--- a/src/thermophysicalModels/basic/basicThermo/makeBasicThermo.H
+++ b/src/thermophysicalModels/basic/basicThermo/makeBasicThermo.H
@@ -38,13 +38,6 @@ Description
 
 #define makeBasicThermo(Cthermo,Mixture,Transport,Thermo,EqnOfState)          \
                                                                               \
-typedef Mixture<Transport<specieThermo<Thermo<EqnOfState> > > >               \
-    Mixture##Transport##Thermo##EqnOfState;                                   \
-                                                                              \
-defineTemplateTypeNameAndDebugWithName                                        \
-    (Mixture##Transport##Thermo##EqnOfState,                                  \
-    #Mixture"<"#Transport"<specieThermo<"#Thermo"<"#EqnOfState">>>>", 0)      \
-                                                                              \
 typedef Cthermo<Mixture<Transport<specieThermo<Thermo<EqnOfState> > > > >     \
     Cthermo##Mixture##Transport##Thermo##EqnOfState;                          \
                                                                               \
@@ -60,7 +53,6 @@ addToRunTimeSelectionTable                                                    \
     fvMesh                                                                    \
 )
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/fixedInternalEnergy/fixedInternalEnergyFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/fixedInternalEnergy/fixedInternalEnergyFvPatchScalarField.C
new file mode 100644
index 0000000000000000000000000000000000000000..1ad0a1c0d5aba0217424d04a45fbbb0b780218fc
--- /dev/null
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/fixedInternalEnergy/fixedInternalEnergyFvPatchScalarField.C
@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 "fixedInternalEnergyFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "basicThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+fixedInternalEnergyFvPatchScalarField::fixedInternalEnergyFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF)
+{}
+
+
+fixedInternalEnergyFvPatchScalarField::fixedInternalEnergyFvPatchScalarField
+(
+    const fixedInternalEnergyFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper)
+{}
+
+
+fixedInternalEnergyFvPatchScalarField::fixedInternalEnergyFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict)
+{}
+
+
+fixedInternalEnergyFvPatchScalarField::fixedInternalEnergyFvPatchScalarField
+(
+    const fixedInternalEnergyFvPatchScalarField& tppsf
+)
+:
+    fixedValueFvPatchScalarField(tppsf)
+{}
+
+
+fixedInternalEnergyFvPatchScalarField::fixedInternalEnergyFvPatchScalarField
+(
+    const fixedInternalEnergyFvPatchScalarField& tppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(tppsf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void fixedInternalEnergyFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const basicThermo& thermo = db().lookupObject<basicThermo>
+    (
+        "thermophysicalProperties"
+    );
+    
+    const label patchi = patch().index();
+
+    fvPatchScalarField& Tw = 
+        const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);
+    Tw.evaluate();
+    operator==(thermo.e(Tw, patchi));
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, fixedInternalEnergyFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/fixedInternalEnergy/fixedInternalEnergyFvPatchScalarField.H b/src/thermophysicalModels/basic/derivedFvPatchFields/fixedInternalEnergy/fixedInternalEnergyFvPatchScalarField.H
new file mode 100644
index 0000000000000000000000000000000000000000..ae335ed2efd60aa374a5b762e2663dea5dbfe34f
--- /dev/null
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/fixedInternalEnergy/fixedInternalEnergyFvPatchScalarField.H
@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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::fixedInternalEnergyFvPatchScalarField
+
+Description
+    A fixed boundary condition for internal energy
+
+SourceFiles
+    fixedInternalEnergyFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedInternalEnergyFvPatchScalarField_H
+#define fixedInternalEnergyFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                Class fixedInternalEnergyFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class fixedInternalEnergyFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("fixedInternalEnergy");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        fixedInternalEnergyFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        fixedInternalEnergyFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given fixedInternalEnergyFvPatchScalarField
+        // onto a new patch
+        fixedInternalEnergyFvPatchScalarField
+        (
+            const fixedInternalEnergyFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        fixedInternalEnergyFvPatchScalarField
+        (
+            const fixedInternalEnergyFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new fixedInternalEnergyFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        fixedInternalEnergyFvPatchScalarField
+        (
+            const fixedInternalEnergyFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new fixedInternalEnergyFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/gradientInternalEnergy/gradientInternalEnergyFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientInternalEnergy/gradientInternalEnergyFvPatchScalarField.C
new file mode 100644
index 0000000000000000000000000000000000000000..0479eed1e839f62665d7aad0d538ab6f5b7b21b2
--- /dev/null
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientInternalEnergy/gradientInternalEnergyFvPatchScalarField.C
@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 "gradientInternalEnergyFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "basicThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+gradientInternalEnergyFvPatchScalarField::gradientInternalEnergyFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedGradientFvPatchScalarField(p, iF)
+{}
+
+
+gradientInternalEnergyFvPatchScalarField::gradientInternalEnergyFvPatchScalarField
+(
+    const gradientInternalEnergyFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedGradientFvPatchScalarField(ptf, p, iF, mapper)
+{}
+
+
+gradientInternalEnergyFvPatchScalarField::gradientInternalEnergyFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedGradientFvPatchScalarField(p, iF, dict)
+{}
+
+
+gradientInternalEnergyFvPatchScalarField::gradientInternalEnergyFvPatchScalarField
+(
+    const gradientInternalEnergyFvPatchScalarField& tppsf
+)
+:
+    fixedGradientFvPatchScalarField(tppsf)
+{}
+
+
+gradientInternalEnergyFvPatchScalarField::gradientInternalEnergyFvPatchScalarField
+(
+    const gradientInternalEnergyFvPatchScalarField& tppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedGradientFvPatchScalarField(tppsf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void gradientInternalEnergyFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const basicThermo& thermo = db().lookupObject<basicThermo>
+    (
+        "thermophysicalProperties"
+    );
+    
+    const label patchi = patch().index();
+
+    fvPatchScalarField& Tw = 
+        const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);
+
+    Tw.evaluate();
+
+    gradient() = thermo.Cv(Tw, patchi)*Tw.snGrad()
+      + patch().deltaCoeffs()*
+        (
+            thermo.e(Tw, patchi)
+          - thermo.e(Tw, patch().faceCells())
+        );
+
+    fixedGradientFvPatchScalarField::updateCoeffs();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, gradientInternalEnergyFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/gradientInternalEnergy/gradientInternalEnergyFvPatchScalarField.H b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientInternalEnergy/gradientInternalEnergyFvPatchScalarField.H
new file mode 100644
index 0000000000000000000000000000000000000000..562a7aa0c358cf730e1f91bd97dfc84c8887ff68
--- /dev/null
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientInternalEnergy/gradientInternalEnergyFvPatchScalarField.H
@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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::gradientInternalEnergyFvPatchScalarField
+
+Description
+    Gradient boundary condition for internal energy
+
+SourceFiles
+    gradientInternalEnergyFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef gradientInternalEnergyFvPatchScalarField_H
+#define gradientInternalEnergyFvPatchScalarField_H
+
+#include "fixedGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+               Class gradientInternalEnergyFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class gradientInternalEnergyFvPatchScalarField
+:
+    public fixedGradientFvPatchScalarField
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("gradientInternalEnergy");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        gradientInternalEnergyFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        gradientInternalEnergyFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given gradientInternalEnergyFvPatchScalarField
+        // onto a new patch
+        gradientInternalEnergyFvPatchScalarField
+        (
+            const gradientInternalEnergyFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        gradientInternalEnergyFvPatchScalarField
+        (
+            const gradientInternalEnergyFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new gradientInternalEnergyFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        gradientInternalEnergyFvPatchScalarField
+        (
+            const gradientInternalEnergyFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new gradientInternalEnergyFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/mixedInternalEnergy/mixedInternalEnergyFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/mixedInternalEnergy/mixedInternalEnergyFvPatchScalarField.C
new file mode 100644
index 0000000000000000000000000000000000000000..83f80b6890bd4b5f49741a9b669d13eb00ea8502
--- /dev/null
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/mixedInternalEnergy/mixedInternalEnergyFvPatchScalarField.C
@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 "mixedInternalEnergyFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "basicThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+mixedInternalEnergyFvPatchScalarField::mixedInternalEnergyFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(p, iF)
+{
+    valueFraction() = 0.0;
+    refValue() = 0.0;
+    refGrad() = 0.0;
+}
+
+
+mixedInternalEnergyFvPatchScalarField::mixedInternalEnergyFvPatchScalarField
+(
+    const mixedInternalEnergyFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchScalarField(ptf, p, iF, mapper)
+{}
+
+
+mixedInternalEnergyFvPatchScalarField::mixedInternalEnergyFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchScalarField(p, iF, dict)
+{}
+
+
+mixedInternalEnergyFvPatchScalarField::mixedInternalEnergyFvPatchScalarField
+(
+    const mixedInternalEnergyFvPatchScalarField& tppsf
+)
+:
+    mixedFvPatchScalarField(tppsf)
+{}
+
+
+mixedInternalEnergyFvPatchScalarField::mixedInternalEnergyFvPatchScalarField
+(
+    const mixedInternalEnergyFvPatchScalarField& tppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(tppsf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void mixedInternalEnergyFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const basicThermo& thermo = db().lookupObject<basicThermo>
+    (
+        "thermophysicalProperties"
+    );
+    
+    const label patchi = patch().index();
+
+    mixedFvPatchScalarField& Tw = refCast<mixedFvPatchScalarField>
+    (
+        const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi])
+    );
+
+    Tw.evaluate();
+
+    valueFraction() = Tw.valueFraction();
+    refValue() = thermo.e(Tw.refValue(), patchi);
+    refGrad() = thermo.Cv(Tw, patchi)*Tw.refGrad()
+      + patch().deltaCoeffs()*
+        (
+            thermo.e(Tw, patchi)
+          - thermo.e(Tw, patch().faceCells())
+        );
+
+    mixedFvPatchScalarField::updateCoeffs();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, mixedInternalEnergyFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/mixedInternalEnergy/mixedInternalEnergyFvPatchScalarField.H b/src/thermophysicalModels/basic/derivedFvPatchFields/mixedInternalEnergy/mixedInternalEnergyFvPatchScalarField.H
new file mode 100644
index 0000000000000000000000000000000000000000..ac6661d609268bf9cfabf1da0cce19ce6791492b
--- /dev/null
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/mixedInternalEnergy/mixedInternalEnergyFvPatchScalarField.H
@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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::mixedInternalEnergyFvPatchScalarField
+
+Description
+    Mixed boundary conditions for internal energy
+
+SourceFiles
+    mixedInternalEnergyFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mixedInternalEnergyFvPatchScalarField_H
+#define mixedInternalEnergyFvPatchScalarField_H
+
+#include "mixedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+               Class mixedInternalEnergyFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class mixedInternalEnergyFvPatchScalarField
+:
+    public mixedFvPatchScalarField
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("mixedInternalEnergy");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        mixedInternalEnergyFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        mixedInternalEnergyFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given mixedInternalEnergyFvPatchScalarField
+        // onto a new patch
+        mixedInternalEnergyFvPatchScalarField
+        (
+            const mixedInternalEnergyFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        mixedInternalEnergyFvPatchScalarField
+        (
+            const mixedInternalEnergyFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mixedInternalEnergyFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        mixedInternalEnergyFvPatchScalarField
+        (
+            const mixedInternalEnergyFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mixedInternalEnergyFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/eThermo/eThermo.C b/src/thermophysicalModels/basic/eThermo/eThermo.C
new file mode 100644
index 0000000000000000000000000000000000000000..b4629a1d15a4bd473b27d44c7d21453bdb690382
--- /dev/null
+++ b/src/thermophysicalModels/basic/eThermo/eThermo.C
@@ -0,0 +1,343 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 "eThermo.H"
+#include "fvMesh.H"
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::eThermo<MixtureType>::eThermo(const fvMesh& mesh)
+:
+    basicThermo(mesh),
+    MixtureType(*this, mesh),
+
+    e_
+    (
+        IOobject
+        (
+            "e",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionSet(0, 2, -2, 0, 0),
+        eBoundaryTypes()
+    )
+{
+    scalarField& eCells = e_.internalField();
+    const scalarField& TCells = T_.internalField();
+
+    forAll(eCells, celli)
+    {
+        eCells[celli] = this->cellMixture(celli).E(TCells[celli]);
+    }
+
+    forAll(e_.boundaryField(), patchi)
+    {
+        e_.boundaryField()[patchi] == e(T_.boundaryField()[patchi], patchi);
+    }
+
+    eBoundaryCorrection(e_);
+
+    calculate();
+    psi_.oldTime();   // Switch on saving old time
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::eThermo<MixtureType>::~eThermo()
+{}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::eThermo<MixtureType>::calculate()
+{
+    const scalarField& eCells = e_.internalField();
+    const scalarField& pCells = p_.internalField();
+
+    scalarField& TCells = T_.internalField();
+    scalarField& psiCells = psi_.internalField();
+    scalarField& muCells = mu_.internalField();
+    scalarField& alphaCells = alpha_.internalField();
+
+    forAll(TCells, celli)
+    {
+        const typename MixtureType::thermoType& mixture_ =
+            this->cellMixture(celli);
+
+        TCells[celli] = mixture_.TE(eCells[celli], TCells[celli]);
+        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
+
+        muCells[celli] = mixture_.mu(TCells[celli]);
+        alphaCells[celli] = mixture_.alpha(TCells[celli]);
+    }
+
+    forAll(T_.boundaryField(), patchi)
+    {
+        fvPatchScalarField& pp = p_.boundaryField()[patchi];
+        fvPatchScalarField& pT = T_.boundaryField()[patchi];
+        fvPatchScalarField& ppsi = psi_.boundaryField()[patchi];
+
+        fvPatchScalarField& pe = e_.boundaryField()[patchi];
+
+        fvPatchScalarField& pmu = mu_.boundaryField()[patchi];
+        fvPatchScalarField& palpha = alpha_.boundaryField()[patchi];
+
+        if (pT.fixesValue())
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+
+                pe[facei] = mixture_.E(pT[facei]);
+
+                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                palpha[facei] = mixture_.alpha(pT[facei]);
+            }
+        }
+        else
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+
+                pT[facei] = mixture_.TE(pe[facei], pT[facei]);
+
+                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                palpha[facei] = mixture_.alpha(pT[facei]);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::eThermo<MixtureType>::correct()
+{
+    if (debug)
+    {
+        Info<< "entering eThermo<MixtureType>::correct()" << endl;
+    }
+
+    // force the saving of the old-time values
+    psi_.oldTime();
+
+    calculate();
+
+    if (debug)
+    {
+        Info<< "exiting eThermo<MixtureType>::correct()" << endl;
+    }
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::eThermo<MixtureType>::e
+(
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    tmp<scalarField> te(new scalarField(T.size()));
+    scalarField& h = te();
+
+    forAll(T, celli)
+    {
+        h[celli] = this->cellMixture(cells[celli]).E(T[celli]);
+    }
+
+    return te;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::eThermo<MixtureType>::e
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> te(new scalarField(T.size()));
+    scalarField& h = te();
+
+    forAll(T, facei)
+    {
+        h[facei] = this->patchFaceMixture(patchi, facei).E(T[facei]);
+    }
+
+    return te;
+}
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::eThermo<MixtureType>::Cp
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> tCp(new scalarField(T.size()));
+    scalarField& cp = tCp();
+
+    forAll(T, facei)
+    {
+        cp[facei] = this->patchFaceMixture(patchi, facei).Cp(T[facei]);
+    }
+
+    return tCp;
+}
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::eThermo<MixtureType>::Cp() const
+{
+    const fvMesh& mesh = T_.mesh();
+
+    tmp<volScalarField> tCp
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cp",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionSet(0, 2, -2, -1, 0),
+            T_.boundaryField().types()
+        )
+    );
+
+    volScalarField& cp = tCp();
+
+    forAll(T_, celli)
+    {
+        cp[celli] = this->cellMixture(celli).Cp(T_[celli]);
+    }
+
+    forAll(T_.boundaryField(), patchi)
+    {
+        const fvPatchScalarField& pT = T_.boundaryField()[patchi];
+        fvPatchScalarField& pCp = cp.boundaryField()[patchi];
+
+        forAll(pT, facei)
+        {
+            pCp[facei] = this->patchFaceMixture(patchi, facei).Cp(pT[facei]);
+        }
+    }
+
+    return tCp;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::eThermo<MixtureType>::Cv
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> tCv(new scalarField(T.size()));
+    scalarField& cv = tCv();
+
+    forAll(T, facei)
+    {
+        cv[facei] = this->patchFaceMixture(patchi, facei).Cv(T[facei]);
+    }
+
+    return tCv;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::eThermo<MixtureType>::Cv() const
+{
+    const fvMesh& mesh = T_.mesh();
+
+    tmp<volScalarField> tCv
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cv",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionSet(0, 2, -2, -1, 0)
+        )
+    );
+
+    volScalarField& cv = tCv();
+
+    forAll(T_, celli)
+    {
+        cv[celli] = this->cellMixture(celli).Cv(T_[celli]);
+    }
+
+    forAll(T_.boundaryField(), patchi)
+    {
+        cv.boundaryField()[patchi] = Cv(T_.boundaryField()[patchi], patchi);
+    }
+
+    return tCv;
+}
+
+
+template<class MixtureType>
+bool Foam::eThermo<MixtureType>::read()
+{
+    if (basicThermo::read())
+    {
+        MixtureType::read(*this);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/eThermo/eThermo.H b/src/thermophysicalModels/basic/eThermo/eThermo.H
new file mode 100644
index 0000000000000000000000000000000000000000..3071d6150d7cefb646941e3b31ba75f52aa2df6f
--- /dev/null
+++ b/src/thermophysicalModels/basic/eThermo/eThermo.H
@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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::eThermo
+
+Description
+    Internal energy for a mixture
+
+SourceFiles
+    eThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef eThermo_H
+#define eThermo_H
+
+#include "basicThermo.H"
+#include "basicMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class eThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MixtureType>
+class eThermo
+:
+    public basicThermo,
+    public MixtureType
+{
+    // Private data
+
+        volScalarField e_;
+
+
+    // Private member functions
+
+        void calculate();
+
+        //- Construct as copy (not implemented)
+        eThermo(const eThermo<MixtureType>&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("eThermo");
+
+
+    // Constructors
+
+        //- Construct from mese
+        eThermo(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~eThermo();
+
+
+    // Member functions
+
+        //- Return the compostion of the combustion mixture
+        virtual basicMixture& composition()
+        {
+            return *this;
+        }
+
+        //- Return the compostion of the combustion mixture
+        virtual const basicMixture& composition() const
+        {
+            return *this;
+        }
+
+        //- Update properties
+        virtual void correct();
+
+
+        // Access to thermodynamic state variables
+
+            //- Internal energy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& e()
+            {
+                return e_;
+            }
+
+            //- Internal energy [J/kg]
+            virtual const volScalarField& e() const
+            {
+                return e_;
+            }
+
+
+        // Fields derived from thermodynamic state variables
+
+            //- Internal energy for cell-set [J/kg]
+            virtual tmp<scalarField> e
+            (
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+            //- Internal energy for patch [J/kg]
+            virtual tmp<scalarField> e
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure for patch [J/kg/K]
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure [J/kg/K]
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Heat capacity at constant volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cv
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant volume [J/kg/K]
+            virtual tmp<volScalarField> Cv() const;
+
+
+        //- Read thermophysicalProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+#ifdef NoRepository
+#   include "eThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/eThermo/eThermos.C b/src/thermophysicalModels/basic/eThermo/eThermos.C
new file mode 100644
index 0000000000000000000000000000000000000000..fa4f0d220a217cab504baf59931a87093c459316
--- /dev/null
+++ b/src/thermophysicalModels/basic/eThermo/eThermos.C
@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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
+
+Description
+    Selection function for internal energy based thermodynamics package.
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "basicThermo.H"
+#include "makeBasicThermo.H"
+
+#include "perfectGas.H"
+
+#include "hConstThermo.H"
+#include "janafThermo.H"
+#include "specieThermo.H"
+
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+
+#include "eThermo.H"
+#include "pureMixture.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+
+makeBasicThermo
+(
+    eThermo,
+    pureMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicThermo
+(
+    eThermo,
+    pureMixture,
+    sutherlandTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicThermo
+(
+    eThermo,
+    pureMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/hThermo/hThermo.C b/src/thermophysicalModels/basic/hThermo/hThermo.C
index 7e32f7a0636288169e7bc36f02a3b9140c20efba..988a790ce40b2468406fa07594370c43d41eb49b 100644
--- a/src/thermophysicalModels/basic/hThermo/hThermo.C
+++ b/src/thermophysicalModels/basic/hThermo/hThermo.C
@@ -205,7 +205,6 @@ Foam::tmp<Foam::scalarField> Foam::hThermo<MixtureType>::h
     return th;
 }
 
-
 template<class MixtureType>
 Foam::tmp<Foam::scalarField> Foam::hThermo<MixtureType>::Cp
 (
@@ -224,7 +223,6 @@ Foam::tmp<Foam::scalarField> Foam::hThermo<MixtureType>::Cp
     return tCp;
 }
 
-
 template<class MixtureType>
 Foam::tmp<Foam::volScalarField> Foam::hThermo<MixtureType>::Cp() const
 {
@@ -243,7 +241,8 @@ Foam::tmp<Foam::volScalarField> Foam::hThermo<MixtureType>::Cp() const
                 IOobject::NO_WRITE
             ),
             mesh,
-            dimensionSet(0, 2, -2, -1, 0)
+            dimensionSet(0, 2, -2, -1, 0),
+            T_.boundaryField().types()
         )
     );
 
@@ -256,13 +255,38 @@ Foam::tmp<Foam::volScalarField> Foam::hThermo<MixtureType>::Cp() const
 
     forAll(T_.boundaryField(), patchi)
     {
-        cp.boundaryField()[patchi] = Cp(T_.boundaryField()[patchi], patchi);
+        const fvPatchScalarField& pT = T_.boundaryField()[patchi];
+        fvPatchScalarField& pCp = cp.boundaryField()[patchi];
+
+        forAll(pT, facei)
+        {
+            pCp[facei] = this->patchFaceMixture(patchi, facei).Cp(pT[facei]);
+        }
     }
 
     return tCp;
 }
 
 
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hThermo<MixtureType>::Cv
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> tCv(new scalarField(T.size()));
+    scalarField& cv = tCv();
+
+    forAll(T, facei)
+    {
+        cv[facei] = this->patchFaceMixture(patchi, facei).Cv(T[facei]);
+    }
+
+    return tCv;
+}
+
+
 template<class MixtureType>
 Foam::tmp<Foam::volScalarField> Foam::hThermo<MixtureType>::Cv() const
 {
@@ -281,8 +305,7 @@ Foam::tmp<Foam::volScalarField> Foam::hThermo<MixtureType>::Cv() const
                 IOobject::NO_WRITE
             ),
             mesh,
-            dimensionSet(0, 2, -2, -1, 0),
-            T_.boundaryField().types()
+            dimensionSet(0, 2, -2, -1, 0)
         )
     );
 
@@ -295,19 +318,12 @@ Foam::tmp<Foam::volScalarField> Foam::hThermo<MixtureType>::Cv() const
 
     forAll(T_.boundaryField(), patchi)
     {
-        const fvPatchScalarField& pT = T_.boundaryField()[patchi];
-        fvPatchScalarField& pCv = cv.boundaryField()[patchi];
-
-        forAll(pT, facei)
-        {
-            pCv[facei] = this->patchFaceMixture(patchi, facei).Cv(pT[facei]);
-        }
+        cv.boundaryField()[patchi] = Cv(T_.boundaryField()[patchi], patchi);
     }
 
     return tCv;
 }
 
-
 template<class MixtureType>
 bool Foam::hThermo<MixtureType>::read()
 {
diff --git a/src/thermophysicalModels/basic/hThermo/hThermo.H b/src/thermophysicalModels/basic/hThermo/hThermo.H
index b7effc7f702619e0ae2adf763d65938ef6796194..6ed34db0b8046d69e2b06d5cb919f535d0498b6a 100644
--- a/src/thermophysicalModels/basic/hThermo/hThermo.H
+++ b/src/thermophysicalModels/basic/hThermo/hThermo.H
@@ -79,40 +79,39 @@ public:
         hThermo(const fvMesh&);
 
 
-    // Destructor
-
-        ~hThermo();
+    //- Destructor
+    virtual ~hThermo();
 
 
     // Member functions
 
         //- Return the compostion of the combustion mixture
-        basicMixture& composition()
+        virtual basicMixture& composition()
         {
             return *this;
         }
 
         //- Return the compostion of the combustion mixture
-        const basicMixture& composition() const
+        virtual const basicMixture& composition() const
         {
             return *this;
         }
 
         //- Update properties
-        void correct();
+        virtual void correct();
 
 
         // Access to thermodynamic state variables
 
             //- Enthalpy [J/kg]
             //  Non-const access allowed for transport equations
-            volScalarField& h()
+            virtual volScalarField& h()
             {
                 return h_;
             }
 
             //- Enthalpy [J/kg]
-            const volScalarField& h() const
+            virtual const volScalarField& h() const
             {
                 return h_;
             }
@@ -121,31 +120,42 @@ public:
         // Fields derived from thermodynamic state variables
 
             //- Enthalpy for cell-set [J/kg]
-            tmp<scalarField> h
+            virtual tmp<scalarField> h
             (
                 const scalarField& T,
                 const labelList& cells
             ) const;
 
             //- Enthalpy for patch [J/kg]
-            tmp<scalarField> h
+            virtual tmp<scalarField> h
             (
                 const scalarField& T,
                 const label patchi
             ) const;
 
             //- Heat capacity at constant pressure for patch [J/kg/K]
-            tmp<scalarField> Cp(const scalarField& T, const label patchi) const;
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
 
             //- Heat capacity at constant pressure [J/kg/K]
-            tmp<volScalarField> Cp() const;
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Heat capacity at constant volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cv
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
 
             //- Heat capacity at constant volume [J/kg/K]
-            tmp<volScalarField> Cv() const;
+            virtual tmp<volScalarField> Cv() const;
 
 
         //- Read thermophysicalProperties dictionary
-        bool read();
+        virtual bool read();
 };
 
 
diff --git a/src/thermophysicalModels/basic/hThermo/hThermos.C b/src/thermophysicalModels/basic/hThermo/hThermos.C
new file mode 100644
index 0000000000000000000000000000000000000000..6b144725c0c8494ba51f78914ac5667198a0d373
--- /dev/null
+++ b/src/thermophysicalModels/basic/hThermo/hThermos.C
@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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
+
+Description
+    Selection function for internal energy based thermodynamics package.
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "basicThermo.H"
+#include "makeBasicThermo.H"
+
+#include "perfectGas.H"
+
+#include "hConstThermo.H"
+#include "janafThermo.H"
+#include "specieThermo.H"
+
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+
+#include "hThermo.H"
+#include "pureMixture.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+
+makeBasicThermo
+(
+    hThermo,
+    pureMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicThermo
+(
+    hThermo,
+    pureMixture,
+    sutherlandTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicThermo
+(
+    hThermo,
+    pureMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C b/src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C
new file mode 100644
index 0000000000000000000000000000000000000000..8efab94648ef57acc383eb758efecca90ca371b8
--- /dev/null
+++ b/src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C
@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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
+
+Description
+    Selection function for internal energy based thermodynamics package.
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "basicMixture.H"
+#include "makeBasicMixture.H"
+
+#include "perfectGas.H"
+
+#include "hConstThermo.H"
+#include "janafThermo.H"
+#include "specieThermo.H"
+
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+
+#include "pureMixture.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+makeBasicMixture
+(
+    pureMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/mixtures/basicMixture/makeBasicMixture.H b/src/thermophysicalModels/basic/mixtures/basicMixture/makeBasicMixture.H
new file mode 100644
index 0000000000000000000000000000000000000000..6214aa0814c960dbd4372bc4cc38e1923e4ea66e
--- /dev/null
+++ b/src/thermophysicalModels/basic/mixtures/basicMixture/makeBasicMixture.H
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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
+
+InClass
+    Foam::basicMixture
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeBasicMixture_H
+#define makeBasicMixture_H
+
+#include "basicMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeBasicMixture(Mixture,Transport,Thermo,EqnOfState)                 \
+                                                                              \
+typedef Mixture<Transport<specieThermo<Thermo<EqnOfState> > > >               \
+    Mixture##Transport##Thermo##EqnOfState;                                   \
+                                                                              \
+defineTemplateTypeNameAndDebugWithName                                        \
+    (Mixture##Transport##Thermo##EqnOfState,                                  \
+    #Mixture"<"#Transport"<specieThermo<"#Thermo"<"#EqnOfState">>>>", 0)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/combustion/hCombustionThermo/hCombustionThermo.H b/src/thermophysicalModels/combustion/hCombustionThermo/hCombustionThermo.H
index c1dd274813e321b242a0dec1b956f51c34a1f744..1ba880983890fffbb516350df0a43e2d26d4d836 100644
--- a/src/thermophysicalModels/combustion/hCombustionThermo/hCombustionThermo.H
+++ b/src/thermophysicalModels/combustion/hCombustionThermo/hCombustionThermo.H
@@ -92,9 +92,8 @@ public:
         static autoPtr<hCombustionThermo> New(const fvMesh&);
 
 
-    // Destructor
-
-        virtual ~hCombustionThermo();
+    //- Destructor
+    virtual ~hCombustionThermo();
 
 
     // Member functions
@@ -110,18 +109,24 @@ public:
 
             //- Enthalpy [J/kg]
             //  Non-const access allowed for transport equations
-            volScalarField& h()
+            virtual volScalarField& h()
             {
                 return h_;
             }
 
             //- Enthalpy [J/kg]
-            const volScalarField& h() const
+            virtual const volScalarField& h() const
             {
                 return h_;
             }
 
 
+        //- Sensible enthalpy [J/kg]
+        virtual tmp<volScalarField> hs() const = 0;
+
+        //- Chemical enthalpy [J/kg]
+        virtual tmp<volScalarField> hc() const = 0;
+
         //- Update properties
         virtual void correct() = 0;
 };
diff --git a/src/thermophysicalModels/combustion/hhuCombustionThermo/hhuCombustionThermo.H b/src/thermophysicalModels/combustion/hhuCombustionThermo/hhuCombustionThermo.H
index ffc9f27f5a45e543d51c1517943c9c698e04dae9..0f235d5dd7ac84dd1b8d7ac895d901bf3c3ca76c 100644
--- a/src/thermophysicalModels/combustion/hhuCombustionThermo/hhuCombustionThermo.H
+++ b/src/thermophysicalModels/combustion/hhuCombustionThermo/hhuCombustionThermo.H
@@ -96,9 +96,8 @@ public:
         static autoPtr<hhuCombustionThermo> New(const fvMesh&);
 
 
-    // Destructor
-
-        virtual ~hhuCombustionThermo();
+    //- Destructor
+    virtual ~hhuCombustionThermo();
 
 
     // Member functions
@@ -111,13 +110,13 @@ public:
 
             //- Unburnt gas enthalpy [J/kg]
             //  Non-const access allowed for transport equations
-            volScalarField& hu()
+            virtual volScalarField& hu()
             {
                 return hu_;
             }
 
             //- Unburnt gas enthalpy [J/kg]
-            const volScalarField& hu() const
+            virtual const volScalarField& hu() const
             {
                 return hu_;
             }
@@ -140,7 +139,7 @@ public:
             ) const = 0;
 
             //- Unburnt gas temperature [K]
-            const volScalarField& Tu() const
+            virtual const volScalarField& Tu() const
             {
                 return Tu_;
             }
diff --git a/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C b/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C
index 064e7c3ebafeb0f67be03f7d547df9e1d2c31f55..936973abda7cb9ce54bde765f691cfe97d42c8a9 100644
--- a/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C
+++ b/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C
@@ -134,7 +134,120 @@ void Foam::hMixtureThermo<MixtureType>::calculate()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::h
+void Foam::hMixtureThermo<MixtureType>::correct()
+{
+    if (debug)
+    {
+        Info<< "entering hMixtureThermo<MixtureType>::correct()" << endl;
+    }
+
+    // force the saving of the old-time values
+    psi_.oldTime();
+
+    calculate();
+
+    if (debug)
+    {
+        Info<< "exiting hMixtureThermo<MixtureType>::correct()" << endl;
+    }
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField>
+Foam::hMixtureThermo<MixtureType>::hs() const
+{
+    const fvMesh& mesh = T_.mesh();
+
+    tmp<volScalarField> ths
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "hs",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            h_.dimensions()
+        )
+    );
+
+    volScalarField& hsf = ths();
+    scalarField& hsCells = hsf.internalField();
+    const scalarField& TCells = T_.internalField();
+
+    forAll(TCells, celli)
+    {
+        hsCells[celli] = this->cellMixture(celli).Hs(TCells[celli]);
+    }
+
+    forAll(T_.boundaryField(), patchi)
+    {
+        scalarField& hsp = hsf.boundaryField()[patchi];
+        const scalarField& Tp = T_.boundaryField()[patchi];
+
+        forAll(Tp, facei)
+        {
+            hsp[facei] = this->patchFaceMixture(patchi, facei).Hs(Tp[facei]);
+        }
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField>
+Foam::hMixtureThermo<MixtureType>::hc() const
+{
+    const fvMesh& mesh = T_.mesh();
+
+    tmp<volScalarField> thc
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "hc",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            h_.dimensions()
+        )
+    );
+
+    volScalarField& hcf = thc();
+    scalarField& hcCells = hcf.internalField();
+
+    forAll(hcCells, celli)
+    {
+        hcCells[celli] = this->cellMixture(celli).Hc();
+    }
+
+    forAll(hcf.boundaryField(), patchi)
+    {
+        scalarField& hcp = hcf.boundaryField()[patchi];
+
+        forAll(hcp, facei)
+        {
+            hcp[facei] = this->patchFaceMixture(patchi, facei).Hc();
+        }
+    }
+
+    return thc;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField>
+Foam::hMixtureThermo<MixtureType>::h
 (
     const scalarField& T,
     const labelList& cells
@@ -153,7 +266,8 @@ Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::h
 
 
 template<class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField>
+Foam::hMixtureThermo<MixtureType>::h
 (
     const scalarField& T,
     const label patchi
@@ -172,7 +286,8 @@ Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::h
 
 
 template<class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::Cp
+Foam::tmp<Foam::scalarField>
+Foam::hMixtureThermo<MixtureType>::Cp
 (
     const scalarField& T,
     const label patchi
@@ -192,7 +307,8 @@ Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::Cp
 
 
 template<class MixtureType>
-Foam::tmp<Foam::volScalarField> Foam::hMixtureThermo<MixtureType>::Cp() const
+Foam::tmp<Foam::volScalarField>
+Foam::hMixtureThermo<MixtureType>::Cp() const
 {
     const fvMesh& mesh = T_.mesh();
 
@@ -215,9 +331,12 @@ Foam::tmp<Foam::volScalarField> Foam::hMixtureThermo<MixtureType>::Cp() const
 
     volScalarField& cp = tCp();
 
-    forAll(T_, celli)
+    scalarField& cpCells = cp.internalField();
+    const scalarField& TCells = T_.internalField();
+
+    forAll(TCells, celli)
     {
-        cp[celli] = this->cellMixture(celli).Cp(T_[celli]);
+        cpCells[celli] = this->cellMixture(celli).Cp(TCells[celli]);
     }
 
     forAll(T_.boundaryField(), patchi)
@@ -229,26 +348,6 @@ Foam::tmp<Foam::volScalarField> Foam::hMixtureThermo<MixtureType>::Cp() const
 }
 
 
-template<class MixtureType>
-void Foam::hMixtureThermo<MixtureType>::correct()
-{
-    if (debug)
-    {
-        Info<< "entering hMixtureThermo<MixtureType>::correct()" << endl;
-    }
-
-    // force the saving of the old-time values
-    psi_.oldTime();
-
-    calculate();
-
-    if (debug)
-    {
-        Info<< "exiting hMixtureThermo<MixtureType>::correct()" << endl;
-    }
-}
-
-
 template<class MixtureType>
 bool Foam::hMixtureThermo<MixtureType>::read()
 {
diff --git a/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.H b/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.H
index 2436ff833d53186c7dfb1bd85eb0d4058a39d8e5..e1b231bac24286be9d1361017dc7e8f5fcd00da8 100644
--- a/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.H
+++ b/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.H
@@ -73,21 +73,20 @@ public:
         hMixtureThermo(const fvMesh&);
 
 
-    // Destructor
-
-        virtual ~hMixtureThermo();
+    //- Destructor
+    virtual ~hMixtureThermo();
 
 
     // Member functions
 
         //- Return the compostion of the combustion mixture
-        combustionMixture& composition()
+        virtual combustionMixture& composition()
         {
             return *this;
         }
 
         //- Return the compostion of the combustion mixture
-        const combustionMixture& composition() const
+        virtual const combustionMixture& composition() const
         {
             return *this;
         }
@@ -95,28 +94,37 @@ public:
         //- Update properties
         virtual void correct();
 
+        //- Sensible enthalpy [J/kg]
+        virtual tmp<volScalarField> hs() const;
+
+        //- Chemical enthalpy [J/kg]
+        virtual tmp<volScalarField> hc() const;
 
         // Fields derived from thermodynamic state variables
 
             //- Enthalpy for cell-set [J/kg]
-            tmp<scalarField> h
+            virtual tmp<scalarField> h
             (
                 const scalarField& T,
                 const labelList& cells
             ) const;
 
             //- Enthalpy for patch [J/kg]
-            tmp<scalarField> h
+            virtual tmp<scalarField> h
             (
                 const scalarField& T,
                 const label patchi
             ) const;
 
             //- Heat capacity at constant pressure for patch [J/kg/K]
-            tmp<scalarField> Cp(const scalarField& T, const label patchi) const;
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
 
             //- Heat capacity at constant pressure [J/kg/K]
-            tmp<volScalarField> Cp() const;
+            virtual tmp<volScalarField> Cp() const;
 
 
         //- Read thermophysicalProperties dictionary
diff --git a/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C b/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
index 0086e8cd147c2824097599d087c7c6537e0787b3..694d2e3b60df0a6b7a189494e5b4bc87fdeb4da1 100644
--- a/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
+++ b/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
@@ -173,6 +173,100 @@ void Foam::hhuMixtureThermo<MixtureType>::correct()
     }
 }
 
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField>
+Foam::hhuMixtureThermo<MixtureType>::hs() const
+{
+    const fvMesh& mesh = T_.mesh();
+
+    tmp<volScalarField> ths
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "hs",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            h_.dimensions()
+        )
+    );
+
+    volScalarField& hsf = ths();
+
+    scalarField& hsCells = hsf.internalField();
+    const scalarField& TCells = T_.internalField();
+
+    forAll(TCells, celli)
+    {
+        hsCells[celli] = this->cellMixture(celli).Hs(TCells[celli]);
+    }
+
+    forAll(T_.boundaryField(), patchi)
+    {
+        scalarField& hsp = hsf.boundaryField()[patchi];
+        const scalarField& Tp = T_.boundaryField()[patchi];
+
+        forAll(Tp, facei)
+        {
+            hsp[facei] = this->patchFaceMixture(patchi, facei).Hs(Tp[facei]);
+        }
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField>
+Foam::hhuMixtureThermo<MixtureType>::hc() const
+{
+    const fvMesh& mesh = T_.mesh();
+
+    tmp<volScalarField> thc
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "hc",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            h_.dimensions()
+        )
+    );
+
+    volScalarField& hcf = thc();
+    scalarField& hcCells = hcf.internalField();
+
+    forAll(hcCells, celli)
+    {
+        hcCells[celli] = this->cellMixture(celli).Hc();
+    }
+
+    forAll(hcf.boundaryField(), patchi)
+    {
+        scalarField& hcp = hcf.boundaryField()[patchi];
+
+        forAll(hcp, facei)
+        {
+            hcp[facei] = this->patchFaceMixture(patchi, facei).Hc();
+        }
+    }
+
+    return thc;
+}
+
+
 template<class MixtureType>
 Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::h
 (
@@ -232,7 +326,8 @@ Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::Cp
 
 
 template<class MixtureType>
-Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::Cp() const
+Foam::tmp<Foam::volScalarField>
+Foam::hhuMixtureThermo<MixtureType>::Cp() const
 {
     const fvMesh& mesh = T_.mesh();
 
@@ -254,10 +349,12 @@ Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::Cp() const
     );
 
     volScalarField& cp = tCp();
+    scalarField& cpCells = cp.internalField();
+    const scalarField& TCells = T_.internalField();
 
-    forAll(T_, celli)
+    forAll(TCells, celli)
     {
-        cp[celli] = this->cellMixture(celli).Cp(T_[celli]);
+        cpCells[celli] = this->cellMixture(celli).Cp(TCells[celli]);
     }
 
     forAll(T_.boundaryField(), patchi)
@@ -269,9 +366,9 @@ Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::Cp() const
 }
 
 
-
 template<class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::hu
+Foam::tmp<Foam::scalarField>
+Foam::hhuMixtureThermo<MixtureType>::hu
 (
     const scalarField& Tu,
     const labelList& cells
@@ -290,7 +387,8 @@ Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::hu
 
 
 template<class MixtureType>
-Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::hu
+Foam::tmp<Foam::scalarField>
+Foam::hhuMixtureThermo<MixtureType>::hu
 (
     const scalarField& Tu,
     const label patchi
@@ -309,7 +407,8 @@ Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::hu
 
 
 template<class MixtureType>
-Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::Tb() const
+Foam::tmp<Foam::volScalarField>
+Foam::hhuMixtureThermo<MixtureType>::Tb() const
 {
     tmp<volScalarField> tTb
     (
@@ -328,10 +427,14 @@ Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::Tb() const
     );
 
     volScalarField& Tb_ = tTb();
+    scalarField& TbCells = Tb_.internalField();
+    const scalarField& TCells = T_.internalField();
+    const scalarField& hCells = h_.internalField();
 
-    forAll(Tb_, celli)
+    forAll(TbCells, celli)
     {
-        Tb_[celli] = this->cellProducts(celli).TH(h_[celli], T_[celli]);
+        TbCells[celli] =
+            this->cellProducts(celli).TH(hCells[celli], TCells[celli]);
     }
 
     forAll(Tb_.boundaryField(), patchi)
@@ -344,7 +447,8 @@ Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::Tb() const
         forAll(pTb, facei)
         {
             pTb[facei] =
-                this->patchFaceProducts(patchi, facei).TH(ph[facei], pT[facei]);
+                this->patchFaceProducts(patchi, facei)
+               .TH(ph[facei], pT[facei]);
         }
     }
 
@@ -374,10 +478,14 @@ Foam::hhuMixtureThermo<MixtureType>::psiu() const
     );
 
     volScalarField& psiu = tpsiu();
+    scalarField& psiuCells = psiu.internalField();
+    const scalarField& TuCells = Tu_.internalField();
+    const scalarField& pCells = p_.internalField();
 
-    forAll(psiu, celli)
+    forAll(psiuCells, celli)
     {
-        psiu[celli] = this->cellReactants(celli).psi(p_[celli], Tu_[celli]);
+        psiuCells[celli] =
+            this->cellReactants(celli).psi(pCells[celli], TuCells[celli]);
     }
 
     forAll(psiu.boundaryField(), patchi)
@@ -421,12 +529,15 @@ Foam::hhuMixtureThermo<MixtureType>::psib() const
     );
 
     volScalarField& psib = tpsib();
-
+    scalarField& psibCells = psib.internalField();
     volScalarField Tb_ = Tb();
+    const scalarField& TbCells = Tb_.internalField();
+    const scalarField& pCells = p_.internalField();
 
-    forAll(psib, celli)
+    forAll(psibCells, celli)
     {
-        psib[celli] = this->cellReactants(celli).psi(p_[celli], Tb_[celli]);
+        psibCells[celli] =
+            this->cellReactants(celli).psi(pCells[celli], TbCells[celli]);
     }
 
     forAll(psib.boundaryField(), patchi)
@@ -449,7 +560,8 @@ Foam::hhuMixtureThermo<MixtureType>::psib() const
 
 
 template<class MixtureType>
-Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::muu() const
+Foam::tmp<Foam::volScalarField>
+Foam::hhuMixtureThermo<MixtureType>::muu() const
 {
     tmp<volScalarField> tmuu
     (
@@ -469,10 +581,12 @@ Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::muu() const
     );
 
     volScalarField& muu_ = tmuu();
+    scalarField& muuCells = muu_.internalField();
+    const scalarField& TuCells = Tu_.internalField();
 
-    forAll(muu_, celli)
+    forAll(muuCells, celli)
     {
-        muu_[celli] = this->cellReactants(celli).mu(Tu_[celli]);
+        muuCells[celli] = this->cellReactants(celli).mu(TuCells[celli]);
     }
 
     forAll(muu_.boundaryField(), patchi)
@@ -492,7 +606,8 @@ Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::muu() const
 
 
 template<class MixtureType>
-Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::mub() const
+Foam::tmp<Foam::volScalarField>
+Foam::hhuMixtureThermo<MixtureType>::mub() const
 {
     tmp<volScalarField> tmub
     (
@@ -512,11 +627,13 @@ Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::mub() const
     );
 
     volScalarField& mub_ = tmub();
+    scalarField& mubCells = mub_.internalField();
     volScalarField Tb_ = Tb();
+    const scalarField& TbCells = Tb_.internalField();
 
-    forAll(mub_, celli)
+    forAll(mubCells, celli)
     {
-        mub_[celli] = this->cellProducts(celli).mu(Tb_[celli]);
+        mubCells[celli] = this->cellProducts(celli).mu(TbCells[celli]);
     }
 
     forAll(mub_.boundaryField(), patchi)
@@ -526,7 +643,8 @@ Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::mub() const
 
         forAll(pMub, facei)
         {
-            pMub[facei] = this->patchFaceProducts(patchi, facei).mu(pTb[facei]);
+            pMub[facei] =
+                this->patchFaceProducts(patchi, facei).mu(pTb[facei]);
         }
     }
 
diff --git a/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.H b/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.H
index a5b3855599e497a5bfb1fe89f77cfc128f147971..7f3997ceee5eaf0e0fc4f78e7e72a6a56fc1705f 100644
--- a/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.H
+++ b/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.H
@@ -73,21 +73,20 @@ public:
         hhuMixtureThermo(const fvMesh&);
 
 
-    // Destructor
-
-        virtual ~hhuMixtureThermo();
+    //- Destructor
+    virtual ~hhuMixtureThermo();
 
 
     // Member functions
 
         //- Return the compostion of the combustion mixture
-        combustionMixture& composition()
+        virtual combustionMixture& composition()
         {
             return *this;
         }
 
         //- Return the compostion of the combustion mixture
-        const combustionMixture& composition() const
+        virtual const combustionMixture& composition() const
         {
             return *this;
         }
@@ -95,38 +94,48 @@ public:
         //- Update properties
         virtual void correct();
 
+        //- Sensible enthalpy [J/kg]
+        virtual tmp<volScalarField> hs() const;
+
+        //- Chemical enthalpy [J/kg]
+        virtual tmp<volScalarField> hc() const;
+
 
         // Fields derived from thermodynamic state variables
 
             //- Enthalpy for cell-set [J/kg]
-            tmp<scalarField> h
+            virtual tmp<scalarField> h
             (
                 const scalarField& T,
                 const labelList& cells
             ) const;
 
             //- Enthalpy for patch [J/kg]
-            tmp<scalarField> h
+            virtual tmp<scalarField> h
             (
                 const scalarField& T,
                 const label patchi
             ) const;
 
             //- Heat capacity at constant pressure for patch [J/kg/K]
-            tmp<scalarField> Cp(const scalarField& T, const label patchi) const;
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
 
             //- Heat capacity at constant pressure [J/kg/K]
-            tmp<volScalarField> Cp() const;
+            virtual tmp<volScalarField> Cp() const;
 
             //- Unburnt gas enthalpy for cell-set [J/kg]
-            tmp<scalarField> hu
+            virtual tmp<scalarField> hu
             (
                 const scalarField& T,
                 const labelList& cells
             ) const;
 
             //- Unburnt gas enthalpy for patch [J/kg]
-            tmp<scalarField> hu
+            virtual tmp<scalarField> hu
             (
                 const scalarField& T,
                 const label patchi
@@ -134,22 +143,22 @@ public:
 
 
             //- Burnt gas temperature [K]
-            tmp<volScalarField> Tb() const;
+            virtual tmp<volScalarField> Tb() const;
 
             //- Unburnt gas compressibility [s^2/m^2]
-            tmp<volScalarField> psiu() const;
+            virtual tmp<volScalarField> psiu() const;
 
             //- Burnt gas compressibility [s^2/m^2]
-            tmp<volScalarField> psib() const;
+            virtual tmp<volScalarField> psib() const;
 
 
         // Access to transport variables
 
             //- Dynamic viscosity of unburnt gas [kg/ms]
-            tmp<volScalarField> muu() const;
+            virtual tmp<volScalarField> muu() const;
 
             //- Dynamic viscosity of burnt gas [kg/ms]
-            tmp<volScalarField> mub() const;
+            virtual tmp<volScalarField> mub() const;
 
 
         //- Read thermophysicalProperties dictionary
diff --git a/src/thermophysicalModels/specie/specie/specie.C b/src/thermophysicalModels/specie/specie/specie.C
index ff2b2f67e1d3a63ddf8ca31527cdf96efe37c214..18941b6aea8fb4891dcef48ed77e69653eee01c6 100644
--- a/src/thermophysicalModels/specie/specie/specie.C
+++ b/src/thermophysicalModels/specie/specie/specie.C
@@ -22,35 +22,27 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    Base class of the thermophysical property types.
-
 \*---------------------------------------------------------------------------*/
 
 #include "specie.H"
 #include "IOstreams.H"
 #include "dimensionedConstants.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 /* * * * * * * * * * * * * public constants  * * * * * * * * * * * * */
 
 //- Universal gas constant (default in [J/(kmol K)])
-const scalar specie::RR = dimensionedConstant("R", 8314.51);
+const Foam::scalar Foam::specie::RR = dimensionedConstant("R", 8314.51);
 
 //- Standard pressure (default in [Pa])
-const scalar specie::Pstd = dimensionedConstant("Pstd", 1.0e5);
+const Foam::scalar Foam::specie::Pstd = dimensionedConstant("Pstd", 1.0e5);
 
 //- Standard temperature (default in [K])
-const scalar specie::Tstd = dimensionedConstant("Tstd", 298.15);
+const Foam::scalar Foam::specie::Tstd = dimensionedConstant("Tstd", 298.15);
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-specie::specie(Istream& is)
+Foam::specie::specie(Istream& is)
 :
     name_(is),
     nMoles_(readScalar(is)),
@@ -62,7 +54,7 @@ specie::specie(Istream& is)
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
-Ostream& operator<<(Ostream& os, const specie& st)
+Foam::Ostream& Foam::operator<<(Ostream& os, const specie& st)
 {
     os  << st.name_ << tab
         << st.nMoles_ << tab
@@ -73,8 +65,4 @@ Ostream& operator<<(Ostream& os, const specie& st)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.C b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.C
index a2870c86ab9368866716e2aa945d515351b116b4..4c0dcb5b5a944ed1025cf1fd4a4bcd10c86afdad 100644
--- a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.C
+++ b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.C
@@ -22,27 +22,18 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    Constant properties thermodynamics package derived from the basic
-    thermo package data type specieThermo.
-
 \*---------------------------------------------------------------------------*/
 
 #include "eConstThermo.H"
 #include "IOstreams.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-eConstThermo::eConstThermo(Istream& is)
+Foam::eConstThermo::eConstThermo(Istream& is)
 :
     specieThermo(is),
-    CV(readScalar(is)),
-    Hf(readScalar(is))
+    Cv_(readScalar(is)),
+    Hf_(readScalar(is))
 {
     is.check("eConstThermo::eConstThermo(Istream& is)");
 }
@@ -50,17 +41,13 @@ eConstThermo::eConstThermo(Istream& is)
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
-Ostream& operator<<(Ostream& os, const eConstThermo& ct)
+Foam::Ostream& Foam::operator<<(Ostream& os, const eConstThermo& ct)
 {
-    os << (const specieThermo&)ct << tab << ct.CV << tab << ct.Hf;
+    os << (const specieThermo&)ct << tab << ct.Cv_ << tab << ct.Hf_;
 
     os.check("Ostream& operator<<(Ostream& os, const eConstThermo& ct)");
     return os;
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
index fa8e6249460f11cb82527ddaa69f34b4aa53cd07..4cb1217c0e0bae1394638f1a4a97f32d9129edf9 100644
--- a/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
@@ -95,13 +95,13 @@ class eConstThermo
 {
     // Private data
 
-        scalar CV;
-        scalar Hf;
+        scalar Cv_;
+        scalar Hf_;
 
 
     // Private member functions
 
-        //- construct from components
+        //- Construct from components
         inline eConstThermo
         (
             const specieThermo& st,
@@ -109,6 +109,7 @@ class eConstThermo
             const scalar hf
         );
 
+
 public:
 
     // Constructors
@@ -136,6 +137,12 @@ public:
             //- Enthalpy [J/kmol]
             inline scalar h(const scalar T) const;
 
+            //- Sensible Enthalpy [J/kmol]
+            inline scalar hs(const scalar T) const;
+
+            //- Chemical enthalpy [J/kmol]
+            inline scalar hc() const;
+
             //- Entropy [J/(kmol K)]
             inline scalar s(const scalar T) const;
 
diff --git a/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H b/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
index 2f7d1faf4c0dba56c96cb53c167000bb1becf983..ce40ccbd94cc162a1264a468b8344365415a7527 100644
--- a/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
@@ -34,8 +34,8 @@ inline Foam::eConstThermo::eConstThermo
 )
 :
     specieThermo(st),
-    CV(cv),
-    Hf(hf)
+    Cv_(cv),
+    Hf_(hf)
 {}
 
 
@@ -48,8 +48,8 @@ inline Foam::eConstThermo::eConstThermo
 )
 :
     specieThermo(name, ct),
-    CV(ct.CV),
-    Hf(ct.Hf)
+    Cv_(ct.Cv_),
+    Hf_(ct.Hf_)
 {}
 
 
@@ -79,13 +79,26 @@ Foam::eConstThermo<equationOfState>::New(Istream& is)
 
 inline Foam::scalar Foam::eConstThermo::cp(const scalar) const
 {
-    return CV*W() + RR;
+    return Cv_*W() + RR;
 }
 
 
 inline Foam::scalar Foam::eConstThermo::h(const scalar T) const
 {
-    return cp(T)*T + Hf*W();
+    return cp(T)*T + Hf_*W();
+}
+
+
+inline Foam::scalar Foam::eConstThermo::hs(const scalar T) const
+{
+    return cp(T)*T;
+}
+
+
+template<class equationOfState>
+inline Foam::scalar Foam::eConstThermo<equationOfState>::hc() const
+{
+    return Hf_*this->W();
 }
 
 
@@ -102,7 +115,7 @@ inline Foam::scalar Foam::eConstThermo::TH
     const scalar T0
 ) const
 {
-    return (h - Hf)/Cp(T0);
+    return (h - Hf_)/Cp(T0);
 }
 
 
@@ -112,7 +125,7 @@ inline Foam::scalar Foam::eConstThermo::TE
     const scalar
 ) const
 {
-    return (e - Hf)/CV;
+    return (e - Hf_)/Cv_;
 }
 
 
@@ -125,8 +138,8 @@ inline Foam::eConstThermo& Foam::eConstThermo::operator=
 {
     specieThermo::operator=(ct);
 
-    CV = ct.CV;
-    Hf = ct.Hf;
+    Cv_ = ct.Cv_;
+    Hf_ = ct.Hf_;
 
     return *this;
 }
@@ -145,8 +158,8 @@ inline Foam::eConstThermo Foam::operator+
     return eConstThermo
     (
         st,
-        ct1.nMoles()/st.nMoles()*ct1.CV + ct2.nMoles()/st.nMoles()*ct2.CV,
-        ct1.nMoles()/st.nMoles()*ct1.Hf + ct2.nMoles()/st.nMoles()*ct2.Hf
+        ct1.nMoles()/st.nMoles()*ct1.Cv_ + ct2.nMoles()/st.nMoles()*ct2.Cv_,
+        ct1.nMoles()/st.nMoles()*ct1.Hf_ + ct2.nMoles()/st.nMoles()*ct2.Hf_
     );
 }
 
@@ -162,8 +175,8 @@ inline Foam::eConstThermo Foam::operator-
     return eConstThermo
     (
         st,
-        ct1.nMoles()/st.nMoles()*ct1.CV - ct2.nMoles()/st.nMoles()*ct2.CV,
-        ct1.nMoles()/st.nMoles()*ct1.Hf - ct2.nMoles()/st.nMoles()*ct2.Hf
+        ct1.nMoles()/st.nMoles()*ct1.Cv_ - ct2.nMoles()/st.nMoles()*ct2.Cv_,
+        ct1.nMoles()/st.nMoles()*ct1.Hf_ - ct2.nMoles()/st.nMoles()*ct2.Hf_
     );
 }
 
@@ -177,8 +190,8 @@ inline Foam::eConstThermo Foam::operator*
     return eConstThermo
     (
         s*((const specieThermo&)ct),
-        ct.CV,
-        ct.Hf
+        ct.Cv_,
+        ct.Hf_
     );
 }
 
diff --git a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.C b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.C
index 8b4736a051f97794e017a9606bd9a8d0515a6371..912ae4892fe31e5a6e1e7c0db267134b3bdabea7 100644
--- a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.C
+++ b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.C
@@ -22,28 +22,19 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    Constant properties thermodynamics package
-    templated ito the equationOfState.
-
 \*---------------------------------------------------------------------------*/
 
 #include "hConstThermo.H"
 #include "IOstreams.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class equationOfState>
-hConstThermo<equationOfState>::hConstThermo(Istream& is)
+Foam::hConstThermo<equationOfState>::hConstThermo(Istream& is)
 :
     equationOfState(is),
-    CP(readScalar(is)),
-    Hf(readScalar(is))
+    Cp_(readScalar(is)),
+    Hf_(readScalar(is))
 {
     is.check("hConstThermo::hConstThermo(Istream& is)");
 }
@@ -52,18 +43,18 @@ hConstThermo<equationOfState>::hConstThermo(Istream& is)
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class equationOfState>
-Ostream& operator<<(Ostream& os, const hConstThermo<equationOfState>& ct)
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const hConstThermo<equationOfState>& ct
+)
 {
     os  << static_cast<const equationOfState&>(ct) << tab
-        << ct.CP << tab << ct.Hf;
+        << ct.Cp_ << tab << ct.Hf_;
 
     os.check("Ostream& operator<<(Ostream& os, const hConstThermo& ct)");
     return os;
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
index 1355cb67d3df9882329b60839f15337ec8026941..f1d9018446e728f0300e6e7615605bb58c94c72d 100644
--- a/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
@@ -94,13 +94,13 @@ class hConstThermo
 {
     // Private data
 
-        scalar CP;
-        scalar Hf;
+        scalar Cp_;
+        scalar Hf_;
 
 
     // Private member functions
 
-        //- construct from components
+        //- Construct from components
         inline hConstThermo
         (
             const equationOfState& st,
@@ -136,6 +136,12 @@ public:
             //- Enthalpy [J/kmol]
             inline scalar h(const scalar T) const;
 
+            //- Sensible enthalpy [J/kmol]
+            inline scalar hs(const scalar T) const;
+
+            //- Chemical enthalpy [J/kmol]
+            inline scalar hc() const;
+
             //- Entropy [J/(kmol K)]
             inline scalar s(const scalar T) const;
 
diff --git a/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H b/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
index acc976de05716b5730998ed01b5706e32ccbd013..b628384eb9cb002d03e1b3b57b9a950d8ca2763b 100644
--- a/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
@@ -35,8 +35,8 @@ inline Foam::hConstThermo<equationOfState>::hConstThermo
 )
 :
     equationOfState(st),
-    CP(cp),
-    Hf(hf)
+    Cp_(cp),
+    Hf_(hf)
 {}
 
 
@@ -50,8 +50,8 @@ inline Foam::hConstThermo<equationOfState>::hConstThermo
 )
 :
     equationOfState(name, ct),
-    CP(ct.CP),
-    Hf(ct.Hf)
+    Cp_(ct.Cp_),
+    Hf_(ct.Hf_)
 {}
 
 
@@ -80,21 +80,47 @@ Foam::hConstThermo<equationOfState>::New(Istream& is)
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class equationOfState>
-inline Foam::scalar Foam::hConstThermo<equationOfState>::cp(const scalar) const
+inline Foam::scalar Foam::hConstThermo<equationOfState>::cp
+(
+    const scalar
+) const
+{
+    return Cp_*this->W();
+}
+
+
+template<class equationOfState>
+inline Foam::scalar Foam::hConstThermo<equationOfState>::h
+(
+    const scalar T
+) const
 {
-    return CP*this->W();
+    return (Cp_*T + Hf_)*this->W();
 }
 
 
 template<class equationOfState>
-inline Foam::scalar Foam::hConstThermo<equationOfState>::h(const scalar T) const
+inline Foam::scalar Foam::hConstThermo<equationOfState>::hs
+(
+    const scalar T
+) const
 {
-    return (CP*T + Hf)*this->W();
+    return Cp_*T*this->W();
 }
 
 
 template<class equationOfState>
-inline Foam::scalar Foam::hConstThermo<equationOfState>::s(const scalar T) const
+inline Foam::scalar Foam::hConstThermo<equationOfState>::hc() const
+{
+    return Hf_*this->W();
+}
+
+
+template<class equationOfState>
+inline Foam::scalar Foam::hConstThermo<equationOfState>::s
+(
+    const scalar T
+) const
 {
     notImplemented
     (
@@ -119,8 +145,8 @@ inline void Foam::hConstThermo<equationOfState>::operator+=
     molr1 /= this->nMoles();
     scalar molr2 = ct.nMoles()/this->nMoles();
 
-    CP = molr1*CP + molr2*ct.CP;
-    Hf = molr1*Hf + molr2*ct.Hf;
+    Cp_ = molr1*Cp_ + molr2*ct.Cp_;
+    Hf_ = molr1*Hf_ + molr2*ct.Hf_;
 }
 
 
@@ -137,8 +163,8 @@ inline void Foam::hConstThermo<equationOfState>::operator-=
     molr1 /= this->nMoles();
     scalar molr2 = ct.nMoles()/this->nMoles();
 
-    CP = molr1*CP - molr2*ct.CP;
-    Hf = molr1*Hf - molr2*ct.Hf;
+    Cp_ = molr1*Cp_ - molr2*ct.Cp_;
+    Hf_ = molr1*Hf_ - molr2*ct.Hf_;
 }
 
 
@@ -160,8 +186,10 @@ inline Foam::hConstThermo<equationOfState> Foam::operator+
     return hConstThermo<equationOfState>
     (
         eofs,
-        ct1.nMoles()/eofs.nMoles()*ct1.CP + ct2.nMoles()/eofs.nMoles()*ct2.CP,
-        ct1.nMoles()/eofs.nMoles()*ct1.Hf + ct2.nMoles()/eofs.nMoles()*ct2.Hf
+        ct1.nMoles()/eofs.nMoles()*ct1.Cp_
+      + ct2.nMoles()/eofs.nMoles()*ct2.Cp_,
+        ct1.nMoles()/eofs.nMoles()*ct1.Hf_
+      + ct2.nMoles()/eofs.nMoles()*ct2.Hf_
     );
 }
 
@@ -182,8 +210,10 @@ inline Foam::hConstThermo<equationOfState> Foam::operator-
     return hConstThermo<equationOfState>
     (
         eofs,
-        ct1.nMoles()/eofs.nMoles()*ct1.CP - ct2.nMoles()/eofs.nMoles()*ct2.CP,
-        ct1.nMoles()/eofs.nMoles()*ct1.Hf - ct2.nMoles()/eofs.nMoles()*ct2.Hf
+        ct1.nMoles()/eofs.nMoles()*ct1.Cp_
+      - ct2.nMoles()/eofs.nMoles()*ct2.Cp_,
+        ct1.nMoles()/eofs.nMoles()*ct1.Hf_
+      - ct2.nMoles()/eofs.nMoles()*ct2.Hf_
     );
 }
 
@@ -198,8 +228,8 @@ inline Foam::hConstThermo<equationOfState> Foam::operator*
     return hConstThermo<equationOfState>
     (
         s*static_cast<const equationOfState&>(ct),
-        ct.CP,
-        ct.Hf
+        ct.Cp_,
+        ct.Hf_
     );
 }
 
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.C b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.C
index 5a3d05b848c56484b4124294abd748a8f088291b..13002ecc08f465b6beb38b641ff2c652cb348752 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.C
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.C
@@ -22,23 +22,15 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    JANAF tables based thermodynamics package templated ito the equationOfState.
-
 \*---------------------------------------------------------------------------*/
 
 #include "janafThermo.H"
 #include "IOstreams.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class equationOfState>
-janafThermo<equationOfState>::janafThermo(Istream& is)
+Foam::janafThermo<equationOfState>::janafThermo(Istream& is)
 :
     equationOfState(is),
     Tlow_(readScalar(is)),
@@ -103,7 +95,11 @@ janafThermo<equationOfState>::janafThermo(Istream& is)
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class equationOfState>
-Ostream& operator<<(Ostream& os, const janafThermo<equationOfState>& jt)
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const janafThermo<equationOfState>& jt
+)
 {
     os  << static_cast<const equationOfState&>(jt) << nl
         << "    " << jt.Tlow_
@@ -145,8 +141,4 @@ Ostream& operator<<(Ostream& os, const janafThermo<equationOfState>& jt)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
index db7807a6420247cc3748bfbf0c02a1cd90e99238..e9f455b0a867096940ad002eb6a5b3de981cae2f 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
@@ -150,6 +150,12 @@ public:
         //- Enthalpy [J/kmol]
         inline scalar h(const scalar T) const;
 
+        //- Sensible enthalpy [J/kmol]
+        inline scalar hs(const scalar T) const;
+
+        //- Chemical enthalpy [J/kmol]
+        inline scalar hc() const;
+
         //- Entropy [J/(kmol K)]
         inline scalar s(const scalar T) const;
 
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
index f9dc5e7a0e0dda1e946b8b8036b534a0b6dd379b..754996f27c6c6cc8797596629309270fd2e0fbc3 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
@@ -25,6 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "janafThermo.H"
+#include "specie.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -130,8 +131,7 @@ inline Foam::scalar Foam::janafThermo<equationOfState>::h
 ) const
 {
     const coeffArray& a = coeffs(T);
-    return
-    this->RR*
+    return this->RR*
     (
         ((((a[4]/5.0*T + a[3]/4.0)*T + a[2]/3.0)*T + a[1]/2.0)*T + a[0])*T
       + a[5]
@@ -139,6 +139,36 @@ inline Foam::scalar Foam::janafThermo<equationOfState>::h
 }
 
 
+template<class equationOfState>
+inline Foam::scalar Foam::janafThermo<equationOfState>::hs
+(
+    const scalar T
+) const
+{
+    const coeffArray& a = coeffs(T);
+    scalar Tr = T - specie::Tstd;
+    return this->RR*
+    (
+        ((((a[4]/5.0*Tr + a[3]/4.0)*Tr + a[2]/3.0)*Tr + a[1]/2.0)*Tr + a[0])*Tr
+    );
+}
+
+
+template<class equationOfState>
+inline Foam::scalar Foam::janafThermo<equationOfState>::hc() const
+{
+    const coeffArray& a = lowCpCoeffs_;
+    const scalar& Tstd = specie::Tstd;
+    return this->RR*
+    (
+        (
+            (((a[4]/5.0*Tstd + a[3]/4.0)*Tstd + a[2]/3.0)*Tstd + a[1]/2.0)*Tstd
+          + a[0]
+        )*Tstd
+    );
+}
+
+
 template<class equationOfState>
 inline Foam::scalar Foam::janafThermo<equationOfState>::s
 (
diff --git a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.C b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.C
index a2aedda84a177efe7d3a67646fe9a041ad867cf5..3fcf4897719ede6a04496814b29fc4a55fcbb95b 100644
--- a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.C
+++ b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.C
@@ -22,34 +22,24 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    Basic thermodynamics type based on the use of fitting functions for
-    cp, h, s obtained from the template argument type thermo.  All other
-    properties are derived from these primitive functions.
-
 \*---------------------------------------------------------------------------*/
 
 #include "specieThermo.H"
 #include "IOstreams.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 /* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
 
 template<class thermo>
-const scalar specieThermo<thermo>::tol_ = 1.0e-4;
+const Foam::scalar Foam::specieThermo<thermo>::tol_ = 1.0e-4;
 
 template<class thermo>
-const int specieThermo<thermo>::maxIter_ = 100;
+const int Foam::specieThermo<thermo>::maxIter_ = 100;
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class thermo>
-specieThermo<thermo>::specieThermo(Istream& is)
+Foam::specieThermo<thermo>::specieThermo(Istream& is)
 :
     thermo(is)
 {
@@ -60,7 +50,7 @@ specieThermo<thermo>::specieThermo(Istream& is)
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class thermo>
-Ostream& operator<<(Ostream& os, const specieThermo<thermo>& st)
+Foam::Ostream& Foam::operator<<(Ostream& os, const specieThermo<thermo>& st)
 {
     os  << static_cast<const thermo&>(st);
 
@@ -69,8 +59,4 @@ Ostream& operator<<(Ostream& os, const specieThermo<thermo>& st)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
index c0d4802c745156e0bd0cef13613f388fc73aa748..ba92b6da15b014533fd4a9e8631be1fae03144c2 100644
--- a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
+++ b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
@@ -140,6 +140,12 @@ public:
             // Enthalpy [J/kmol]
             //scalar h(const scalar) const;
 
+            // Sensible enthalpy [J/kmol]
+            //scalar hs(const scalar) const;
+
+            // Chemical enthalpy [J/kmol]
+            //scalar hc(const scalar) const;
+
             // Entropy [J/(kmol K)]
             //scalar s(const scalar) const;
 
@@ -158,6 +164,9 @@ public:
                 //- Internal energy [J/kmol]
                 inline scalar e(const scalar T) const;
 
+                //- Sensible internal energy [J/kmol]
+                inline scalar es(const scalar T) const;
+
                 //- Gibbs free energy [J/kmol]
                 inline scalar g(const scalar T) const;
 
@@ -176,6 +185,12 @@ public:
                 //- Enthalpy [J/kg]
                 inline scalar H(const scalar T) const;
 
+                //- Sensible enthalpy [J/kg]
+                inline scalar Hs(const scalar T) const;
+
+                //- Chemical enthalpy [J/kg]
+                inline scalar Hc() const;
+
                 //- Entropy [J/(kg K)]
                 inline scalar S(const scalar T) const;
 
diff --git a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
index 30f11baf42ac884c755646c0925b85648c2f507c..13adec738bc942c74ea7441cbfa509fef8c6c62d 100644
--- a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
@@ -24,20 +24,12 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
-
 #include "specieThermo.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// construct from components
 template<class thermo>
-inline specieThermo<thermo>::specieThermo
+inline Foam::specieThermo<thermo>::specieThermo
 (
     const thermo& sp
 )
@@ -46,10 +38,8 @@ inline specieThermo<thermo>::specieThermo
 {}
 
 
-// return the temperature corresponding to the value of the
-// thermodynamic property f, given the function f = F(T) and dF(T)/dT
 template<class thermo>
-inline scalar specieThermo<thermo>::T
+inline Foam::scalar Foam::specieThermo<thermo>::T
 (
     scalar f,
     scalar T0,
@@ -87,9 +77,8 @@ inline scalar specieThermo<thermo>::T
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct as named copy
 template<class thermo>
-inline specieThermo<thermo>::specieThermo
+inline Foam::specieThermo<thermo>::specieThermo
 (
     const word& name,
     const specieThermo& st
@@ -101,105 +90,114 @@ inline specieThermo<thermo>::specieThermo
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// Calculate and return derived properties
-// (These functions need not provided in derived types)
-
-// mole specific properties
-
-// Heat capacities [J/(kmol K)]
 template<class thermo>
-inline scalar specieThermo<thermo>::cv(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::cv(const scalar T) const
 {
     return this->cp(T) - this->RR;
 }
 
-// gamma = cp/cv []
+
 template<class thermo>
-inline scalar specieThermo<thermo>::gamma(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::gamma(const scalar T) const
 {
     scalar CP = this->cp(T);
     return CP/(CP - this->RR);
 }
 
-// Internal energy [J/kmol]
+
 template<class thermo>
-inline scalar specieThermo<thermo>::e(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::e(const scalar T) const
 {
     return this->h(T) - this->RR*(T - this->Tstd);
 }
 
-// Gibbs free energy [J/kmol]
+
 template<class thermo>
-inline scalar specieThermo<thermo>::g(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::es(const scalar T) const
+{
+    return this->hs(T) - this->RR*(T - this->Tstd);
+}
+
+
+template<class thermo>
+inline Foam::scalar Foam::specieThermo<thermo>::g(const scalar T) const
 {
     return this->h(T) - T*this->s(T);
 }
 
-// Helmholtz free energy [J/kmol]
+
 template<class thermo>
-inline scalar specieThermo<thermo>::a(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::a(const scalar T) const
 {
     return this->e(T) - T*this->s(T);
 }
 
 
-// mass specific properties
-
-// Heat capacity at constant pressure [J/(kg K)]
 template<class thermo>
-inline scalar specieThermo<thermo>::Cp(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::Cp(const scalar T) const
 {
     return this->cp(T)/this->W();
 }
 
-// Heat capacity at constant pressure [J/(kg K)]
+
 template<class thermo>
-inline scalar specieThermo<thermo>::Cv(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::Cv(const scalar T) const
 {
     return this->cv(T)/this->W();
 }
 
-// Enthalpy [J/kg]
+
 template<class thermo>
-inline scalar specieThermo<thermo>::H(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::H(const scalar T) const
 {
     return this->h(T)/this->W();
 }
 
-// Entropy [J/kg]
+
+template<class thermo>
+inline Foam::scalar Foam::specieThermo<thermo>::Hs(const scalar T) const
+{
+    return this->hs(T)/this->W();
+}
+
+
+template<class thermo>
+inline Foam::scalar Foam::specieThermo<thermo>::Hc() const
+{
+    return this->hc()/this->W();
+}
+
+
 template<class thermo>
-inline scalar specieThermo<thermo>::S(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::S(const scalar T) const
 {
     return this->s(T)/this->W();
 }
 
-// Internal energy [J/kg]
+
 template<class thermo>
-inline scalar specieThermo<thermo>::E(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::E(const scalar T) const
 {
     return this->e(T)/this->W();
 }
 
-// Gibbs free energy [J/kg]
+
 template<class thermo>
-inline scalar specieThermo<thermo>::G(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::G(const scalar T) const
 {
     return this->g(T)/this->W();
 }
 
-// Helmholtz free energy [J/kg]
+
 template<class thermo>
-inline scalar specieThermo<thermo>::A(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::A(const scalar T) const
 {
     return this->a(T)/this->W();
 }
 
 
-// Equilibrium reaction thermodynamics
-
-// Equilibrium constant []
 template<class thermo>
-inline scalar specieThermo<thermo>::K(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::K(const scalar T) const
 {
     scalar arg = -this->nMoles()*this->g(T)/(this->RR*T);
 
@@ -213,21 +211,16 @@ inline scalar specieThermo<thermo>::K(const scalar T) const
     }
 }
 
-//- Equilibrium constant [] i.t.o. partial pressures
-//  = PIi(pi/Pstd)^nui
-//  For low pressures (where the gas mixture is near perfect) Kp = K
+
 template<class thermo>
-inline scalar specieThermo<thermo>::Kp(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::Kp(const scalar T) const
 {
     return K(T);
 }
 
 
-//- Equilibrium constant i.t.o. concentration
-//  For low pressures (where the gas mixture is near perfect)
-//  Kc = Kp(pstd/(RR*T))^nu
 template<class thermo>
-inline scalar specieThermo<thermo>::Kc(const scalar T) const
+inline Foam::scalar Foam::specieThermo<thermo>::Kc(const scalar T) const
 {
     if (equal(this->nMoles(), SMALL))
     {
@@ -240,11 +233,12 @@ inline scalar specieThermo<thermo>::Kc(const scalar T) const
 }
 
 
-//- Equilibrium constant [] i.t.o. mole-fractions
-//  For low pressures (where the gas mixture is near perfect) 
-//  Kx = Kp(pstd/p)^nui
 template<class thermo>
-inline scalar specieThermo<thermo>::Kx(const scalar T, const scalar p) const
+inline Foam::scalar Foam::specieThermo<thermo>::Kx
+(
+    const scalar T,
+    const scalar p
+) const
 {
     if (equal(this->nMoles(), SMALL))
     {
@@ -257,11 +251,8 @@ inline scalar specieThermo<thermo>::Kx(const scalar T, const scalar p) const
 }
 
 
-//- Equilibrium constant [] i.t.o. number of moles
-//  For low pressures (where the gas mixture is near perfect) 
-//  Kn = Kp(n*pstd/p)^nui where n = number of moles in mixture
 template<class thermo>
-inline scalar specieThermo<thermo>::Kn
+inline Foam::scalar Foam::specieThermo<thermo>::Kn
 (
     const scalar T,
     const scalar p,
@@ -279,19 +270,23 @@ inline scalar specieThermo<thermo>::Kn
 }
 
 
-// Energy->temperature  inversion functions
-
-// Temperature from Enthalpy given an initial temperature T0
 template<class thermo>
-inline scalar specieThermo<thermo>::TH(const scalar h, const scalar T0) const
+inline Foam::scalar Foam::specieThermo<thermo>::TH
+(
+    const scalar h,
+    const scalar T0
+) const
 {
     return T(h, T0, &specieThermo<thermo>::H, &specieThermo<thermo>::Cp);
 }
 
 
-// Temperature from internal energy given an initial temperature T0
 template<class thermo>
-inline scalar specieThermo<thermo>::TE(const scalar e, const scalar T0) const
+inline Foam::scalar Foam::specieThermo<thermo>::TE
+(
+    const scalar e,
+    const scalar T0
+) const
 {
     return T(e, T0, &specieThermo<thermo>::E, &specieThermo<thermo>::Cv);
 }
@@ -300,19 +295,25 @@ inline scalar specieThermo<thermo>::TE(const scalar e, const scalar T0) const
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class thermo>
-inline void specieThermo<thermo>::operator+=(const specieThermo<thermo>& st)
+inline void Foam::specieThermo<thermo>::operator+=
+(
+    const specieThermo<thermo>& st
+)
 {
     thermo::operator+=(st);
 }
 
 template<class thermo>
-inline void specieThermo<thermo>::operator-=(const specieThermo<thermo>& st)
+inline void Foam::specieThermo<thermo>::operator-=
+(
+    const specieThermo<thermo>& st
+)
 {
     thermo::operator-=(st);
 }
 
 template<class thermo>
-inline void specieThermo<thermo>::operator*=(const scalar s)
+inline void Foam::specieThermo<thermo>::operator*=(const scalar s)
 {
     thermo::operator*=(s);
 }
@@ -321,7 +322,7 @@ inline void specieThermo<thermo>::operator*=(const scalar s)
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 template<class thermo>
-inline specieThermo<thermo> operator+
+inline Foam::specieThermo<thermo> Foam::operator+
 (
     const specieThermo<thermo>& st1,
     const specieThermo<thermo>& st2
@@ -335,7 +336,7 @@ inline specieThermo<thermo> operator+
 
 
 template<class thermo>
-inline specieThermo<thermo> operator-
+inline Foam::specieThermo<thermo> Foam::operator-
 (
     const specieThermo<thermo>& st1,
     const specieThermo<thermo>& st2
@@ -349,7 +350,7 @@ inline specieThermo<thermo> operator-
 
 
 template<class thermo>
-inline specieThermo<thermo> operator*
+inline Foam::specieThermo<thermo> Foam::operator*
 (
     const scalar s,
     const specieThermo<thermo>& st
@@ -363,7 +364,7 @@ inline specieThermo<thermo> operator*
 
 
 template<class thermo>
-inline specieThermo<thermo> operator==
+inline Foam::specieThermo<thermo> Foam::operator==
 (
     const specieThermo<thermo>& st1,
     const specieThermo<thermo>& st2
@@ -373,8 +374,4 @@ inline specieThermo<thermo> operator==
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //