diff --git a/src/thermophysicalModels/basic/basicThermo/basicThermo.C b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
index 3de20158ad..f03bcc8289 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.C
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
@@ -143,7 +143,6 @@ void Foam::basicThermo::eBoundaryCorrection(volScalarField& e)
     }
 }
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::basicThermo::basicThermo(const fvMesh& mesh)
@@ -300,6 +299,50 @@ Foam::tmp<Foam::scalarField> Foam::basicThermo::h
 }
 
 
+Foam::volScalarField& Foam::basicThermo::hs()
+{
+    notImplemented("basicThermo::hs()");
+    return const_cast<volScalarField&>(volScalarField::null());
+}
+
+
+const Foam::volScalarField& Foam::basicThermo::hs() const
+{
+    notImplemented("basicThermo::hs() const");
+    return volScalarField::null();
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::basicThermo::hs
+(
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    notImplemented
+    (
+        "basicThermo::hs"
+        "(const scalarField& T, const labelList& cells) const"
+    );
+    return tmp<scalarField>(NULL);
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::basicThermo::hs
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    notImplemented
+    (
+        "basicThermo::hs"
+        "(const scalarField& T, const label patchi) const"
+    );
+    return tmp<scalarField>(NULL);
+}
+
+
 Foam::volScalarField& Foam::basicThermo::e()
 {
     notImplemented("basicThermo::e()");
diff --git a/src/thermophysicalModels/basic/basicThermo/basicThermo.H b/src/thermophysicalModels/basic/basicThermo/basicThermo.H
index 1b54258525..1c907edcf9 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.H
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.H
@@ -141,7 +141,7 @@ public:
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const;
 
-            //- Enthalpy [J/kg]
+            //- Total enthalpy [J/kg]
             //  Non-const access allowed for transport equations
             virtual volScalarField& h();
 
@@ -162,6 +162,27 @@ public:
                 const label patchi
             ) const;
 
+            //- Sensible enthalpy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& hs();
+
+            //- Enthalpy [J/kg]
+            virtual const volScalarField& hs() const;
+
+            //- Enthalpy for cell-set [J/kg]
+            virtual tmp<scalarField> hs
+            (
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+            //- Enthalpy for patch [J/kg]
+            virtual tmp<scalarField> hs
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
             //- Internal energy [J/kg]
             //  Non-const access allowed for transport equations
             virtual volScalarField& e();
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C
index 19e345abe0..4679801fd3 100644
--- a/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C
@@ -30,14 +30,9 @@ License
 #include "volFields.H"
 #include "basicThermo.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
+Foam::fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF
@@ -47,7 +42,7 @@ fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
 {}
 
 
-fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
+Foam::fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
 (
     const fixedEnthalpyFvPatchScalarField& ptf,
     const fvPatch& p,
@@ -59,7 +54,7 @@ fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
 {}
 
 
-fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
+Foam::fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF,
@@ -70,7 +65,7 @@ fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
 {}
 
 
-fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
+Foam::fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
 (
     const fixedEnthalpyFvPatchScalarField& tppsf
 )
@@ -79,7 +74,7 @@ fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
 {}
 
 
-fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
+Foam::fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
 (
     const fixedEnthalpyFvPatchScalarField& tppsf,
     const DimensionedField<scalar, volMesh>& iF
@@ -91,7 +86,7 @@ fixedEnthalpyFvPatchScalarField::fixedEnthalpyFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void fixedEnthalpyFvPatchScalarField::updateCoeffs()
+void Foam::fixedEnthalpyFvPatchScalarField::updateCoeffs()
 {
     if (updated())
     {
@@ -109,7 +104,14 @@ void fixedEnthalpyFvPatchScalarField::updateCoeffs()
         const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);
     Tw.evaluate();
 
-    operator==(thermo.h(Tw, patchi));
+    if (dimensionedInternalField().name() == "h")
+    {
+        operator==(thermo.h(Tw, patchi));
+    }
+    else
+    {
+        operator==(thermo.hs(Tw, patchi));
+    }
 
     fixedValueFvPatchScalarField::updateCoeffs();
 }
@@ -117,10 +119,14 @@ void fixedEnthalpyFvPatchScalarField::updateCoeffs()
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makePatchTypeField(fvPatchScalarField, fixedEnthalpyFvPatchScalarField);
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        fixedEnthalpyFvPatchScalarField
+    );
+}
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
 
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.H b/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.H
index 725ca85df2..654b5c1af6 100644
--- a/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.H
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.H
@@ -44,7 +44,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-              Class fixedEnthalpyFvPatchScalarField Declaration
+                Class fixedEnthalpyFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
 class fixedEnthalpyFvPatchScalarField
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C
index ed3a4a6ca3..cb9b74020f 100644
--- a/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C
@@ -30,14 +30,9 @@ License
 #include "volFields.H"
 #include "basicThermo.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
+Foam::gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF
@@ -47,7 +42,7 @@ gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
 {}
 
 
-gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
+Foam::gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
 (
     const gradientEnthalpyFvPatchScalarField& ptf,
     const fvPatch& p,
@@ -59,7 +54,7 @@ gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
 {}
 
 
-gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
+Foam::gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF,
@@ -70,7 +65,7 @@ gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
 {}
 
 
-gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
+Foam::gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
 (
     const gradientEnthalpyFvPatchScalarField& tppsf
 )
@@ -79,7 +74,7 @@ gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
 {}
 
 
-gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
+Foam::gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
 (
     const gradientEnthalpyFvPatchScalarField& tppsf,
     const DimensionedField<scalar, volMesh>& iF
@@ -91,7 +86,7 @@ gradientEnthalpyFvPatchScalarField::gradientEnthalpyFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void gradientEnthalpyFvPatchScalarField::updateCoeffs()
+void Foam::gradientEnthalpyFvPatchScalarField::updateCoeffs()
 {
     if (updated())
     {
@@ -110,12 +105,24 @@ void gradientEnthalpyFvPatchScalarField::updateCoeffs()
 
     Tw.evaluate();
 
-    gradient() = thermo.Cp(Tw, patchi)*Tw.snGrad()
-      + patch().deltaCoeffs()*
+    if (dimensionedInternalField().name() == "h")
+    {
+        gradient() = thermo.Cp(Tw, patchi)*Tw.snGrad()
+        + patch().deltaCoeffs()*
         (
             thermo.h(Tw, patchi)
           - thermo.h(Tw, patch().faceCells())
         );
+    }
+    else
+    {
+        gradient() = thermo.Cp(Tw, patchi)*Tw.snGrad()
+        + patch().deltaCoeffs()*
+        (
+            thermo.hs(Tw, patchi)
+          - thermo.hs(Tw, patch().faceCells())
+        );
+    }
 
     fixedGradientFvPatchScalarField::updateCoeffs();
 }
@@ -123,10 +130,14 @@ void gradientEnthalpyFvPatchScalarField::updateCoeffs()
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makePatchTypeField(fvPatchScalarField, gradientEnthalpyFvPatchScalarField);
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        gradientEnthalpyFvPatchScalarField
+    );
+}
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
 
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/mixedEnthalpy/mixedEnthalpyFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/mixedEnthalpy/mixedEnthalpyFvPatchScalarField.C
index 3d64c9a0d2..cfe404d38e 100644
--- a/src/thermophysicalModels/basic/derivedFvPatchFields/mixedEnthalpy/mixedEnthalpyFvPatchScalarField.C
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/mixedEnthalpy/mixedEnthalpyFvPatchScalarField.C
@@ -30,14 +30,9 @@ License
 #include "volFields.H"
 #include "basicThermo.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
+Foam::mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF
@@ -51,7 +46,7 @@ mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
 }
 
 
-mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
+Foam::mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
 (
     const mixedEnthalpyFvPatchScalarField& ptf,
     const fvPatch& p,
@@ -63,7 +58,7 @@ mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
 {}
 
 
-mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
+Foam::mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF,
@@ -74,7 +69,7 @@ mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
 {}
 
 
-mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
+Foam::mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
 (
     const mixedEnthalpyFvPatchScalarField& tppsf
 )
@@ -83,7 +78,7 @@ mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
 {}
 
 
-mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
+Foam::mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
 (
     const mixedEnthalpyFvPatchScalarField& tppsf,
     const DimensionedField<scalar, volMesh>& iF
@@ -95,7 +90,7 @@ mixedEnthalpyFvPatchScalarField::mixedEnthalpyFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void mixedEnthalpyFvPatchScalarField::updateCoeffs()
+void Foam::mixedEnthalpyFvPatchScalarField::updateCoeffs()
 {
     if (updated())
     {
@@ -117,13 +112,27 @@ void mixedEnthalpyFvPatchScalarField::updateCoeffs()
     Tw.evaluate();
 
     valueFraction() = Tw.valueFraction();
-    refValue() = thermo.h(Tw.refValue(), patchi);
-    refGrad() = thermo.Cp(Tw, patchi)*Tw.refGrad()
-      + patch().deltaCoeffs()*
-        (
+
+    if (dimensionedInternalField().name() == "h")
+    {
+        refValue() = thermo.h(Tw.refValue(), patchi);
+        refGrad() = thermo.Cp(Tw, patchi)*Tw.refGrad()
+        + patch().deltaCoeffs()*
+         (
             thermo.h(Tw, patchi)
           - thermo.h(Tw, patch().faceCells())
-        );
+         );
+    }
+    else
+    {
+        refValue() = thermo.hs(Tw.refValue(), patchi);
+        refGrad() = thermo.Cp(Tw, patchi)*Tw.refGrad()
+        + patch().deltaCoeffs()*
+         (
+            thermo.hs(Tw, patchi)
+          - thermo.hs(Tw, patch().faceCells())
+         );
+    }
 
     mixedFvPatchScalarField::updateCoeffs();
 }
@@ -131,10 +140,14 @@ void mixedEnthalpyFvPatchScalarField::updateCoeffs()
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makePatchTypeField(fvPatchScalarField, mixedEnthalpyFvPatchScalarField);
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        mixedEnthalpyFvPatchScalarField
+    );
+}
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
 
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.C b/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.C
index be93dc011e..d1f844f592 100644
--- a/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.C
+++ b/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.C
@@ -156,4 +156,22 @@ Foam::tmp<Foam::fvScalarMatrix> Foam::radiation::radiationModel::Sh
 }
 
 
+Foam::tmp<Foam::fvScalarMatrix> Foam::radiation::radiationModel::Shs
+(
+    basicThermo& thermo
+) const
+{
+    volScalarField& hs = thermo.hs();
+    const volScalarField cp = thermo.Cp();
+    const volScalarField T3 = pow3(T_);
+
+    return
+    (
+        Ru()
+      - fvm::Sp(4.0*Rp()*T3/cp, hs)
+      - Rp()*T3*(T_ - 4.0*hs/cp)
+    );
+}
+
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H b/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H
index aae443461e..3e1a8019cf 100644
--- a/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H
+++ b/src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H
@@ -163,7 +163,7 @@ public:
             virtual void calculate() = 0;
 
             //- Read radiationProperties dictionary
-            virtual bool read();
+            virtual bool read() = 0;
 
 
         // Access
@@ -176,6 +176,9 @@ public:
 
             //- Enthalpy source term
             virtual tmp<fvScalarMatrix> Sh(basicThermo& thermo) const;
+
+            //- Sensible enthalpy source term
+            virtual tmp<fvScalarMatrix> Shs(basicThermo& thermo) const;
 };
 
 
@@ -189,4 +192,3 @@ public:
 #endif
 
 // ************************************************************************* //
-
diff --git a/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/hCombustionThermo.H b/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/hCombustionThermo.H
index b4b66b55bd..d6b1ff7069 100644
--- a/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/hCombustionThermo.H
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/hCombustionThermo.H
@@ -128,9 +128,6 @@ public:
             }
 
 
-        //- Sensible enthalpy [J/kg]
-        virtual tmp<volScalarField> hs() const = 0;
-
         //- Chemical enthalpy [J/kg]
         virtual tmp<volScalarField> hc() const = 0;
 
diff --git a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.C b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.C
index 5a788049a6..d354b3a986 100644
--- a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.C
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.C
@@ -155,53 +155,6 @@ void Foam::hPsiMixtureThermo<MixtureType>::correct()
 }
 
 
-template<class MixtureType>
-Foam::tmp<Foam::volScalarField>
-Foam::hPsiMixtureThermo<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::hPsiMixtureThermo<MixtureType>::hc() const
diff --git a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.H b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.H
index 6450328447..a30a0c1465 100644
--- a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.H
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.H
@@ -94,9 +94,6 @@ public:
         //- Update properties
         virtual void correct();
 
-        //- Sensible enthalpy [J/kg]
-        virtual tmp<volScalarField> hs() const;
-
         //- Chemical enthalpy [J/kg]
         virtual tmp<volScalarField> hc() const;
 
diff --git a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
index 0af469ecf9..6d7245c044 100644
--- a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
@@ -174,54 +174,6 @@ 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
diff --git a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.H b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.H
index 5e9834e86c..c55e0a9216 100644
--- a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.H
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.H
@@ -36,8 +36,6 @@ SourceFiles
 #ifndef hhuMixtureThermo_H
 #define hhuMixtureThermo_H
 
-//#include "hPsiMixtureThermo.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -94,9 +92,6 @@ public:
         //- Update properties
         virtual void correct();
 
-        //- Sensible enthalpy [J/kg]
-        virtual tmp<volScalarField> hs() const;
-
         //- Chemical enthalpy [J/kg]
         virtual tmp<volScalarField> hc() const;
 
diff --git a/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/hReactionThermo.H b/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/hReactionThermo.H
index 3d5c37614c..f2319efa4f 100644
--- a/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/hReactionThermo.H
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/hReactionThermo.H
@@ -128,9 +128,6 @@ public:
             }
 
 
-        //- Sensible enthalpy [J/kg]
-        virtual tmp<volScalarField> hs() const = 0;
-
         //- Chemical enthalpy [J/kg]
         virtual tmp<volScalarField> hc() const = 0;
 
diff --git a/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hRhoMixtureThermo/hRhoMixtureThermo.C b/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hRhoMixtureThermo/hRhoMixtureThermo.C
index 9cbf9202de..e12fadf8a1 100644
--- a/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hRhoMixtureThermo/hRhoMixtureThermo.C
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hRhoMixtureThermo/hRhoMixtureThermo.C
@@ -154,53 +154,6 @@ void Foam::hRhoMixtureThermo<MixtureType>::correct()
 }
 
 
-template<class MixtureType>
-Foam::tmp<Foam::volScalarField>
-Foam::hRhoMixtureThermo<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::hRhoMixtureThermo<MixtureType>::hc() const
diff --git a/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hRhoMixtureThermo/hRhoMixtureThermo.H b/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hRhoMixtureThermo/hRhoMixtureThermo.H
index 570f696bac..28e95270d9 100644
--- a/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hRhoMixtureThermo/hRhoMixtureThermo.H
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hRhoMixtureThermo/hRhoMixtureThermo.H
@@ -94,9 +94,6 @@ public:
         //- Update properties
         virtual void correct();
 
-        //- Sensible enthalpy [J/kg]
-        virtual tmp<volScalarField> hs() const;
-
         //- Chemical enthalpy [J/kg]
         virtual tmp<volScalarField> hc() const;
 
diff --git a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
index 560f8f4145..7ec9d97959 100644
--- a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
+++ b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
@@ -242,6 +242,9 @@ public:
             //- Temperature from Enthalpy given an initial temperature T0
             inline scalar TH(const scalar H, const scalar T0) const;
 
+            //- Temperature from sensible Enthalpy given an initial T0
+            inline scalar THs(const scalar Hs, const scalar T0) const;
+
             //- Temperature from internal energy given an initial temperature T0
             inline scalar TE(const scalar E, const scalar T0) const;
 
diff --git a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
index 215ea7b203..788f6b9055 100644
--- a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
@@ -281,6 +281,17 @@ inline Foam::scalar Foam::specieThermo<thermo>::TH
 }
 
 
+template<class thermo>
+inline Foam::scalar Foam::specieThermo<thermo>::THs
+(
+    const scalar hs,
+    const scalar T0
+) const
+{
+    return T(hs, T0, &specieThermo<thermo>::Hs, &specieThermo<thermo>::Cp);
+}
+
+
 template<class thermo>
 inline Foam::scalar Foam::specieThermo<thermo>::TE
 (