diff --git a/applications/solvers/combustion/dieselEngineFoam/createFields.H b/applications/solvers/combustion/dieselEngineFoam/createFields.H
index a75c6a7d07..37aeca4c09 100644
--- a/applications/solvers/combustion/dieselEngineFoam/createFields.H
+++ b/applications/solvers/combustion/dieselEngineFoam/createFields.H
@@ -6,7 +6,7 @@ autoPtr<psiChemistryModel> pChemistry
 );
 psiChemistryModel& chemistry = pChemistry();
 
-hCombustionThermo& thermo = chemistry.thermo();
+hsCombustionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -50,7 +50,7 @@ volVectorField U
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
 const volScalarField& T = thermo.T();
-volScalarField& h = thermo.h();
+volScalarField& hs = thermo.hs();
 
 
 #include "compressibleCreatePhi.H"
@@ -92,4 +92,18 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
+
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);
diff --git a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
index 6507c117a0..2e6a9c4b05 100644
--- a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
+++ b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    dieselFoam
+    dieselEngineFoam
 
 Description
     Solver for diesel engine spray and combustion.
@@ -103,13 +103,15 @@ int main(int argc, char *argv[])
             kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
         }
 
+        chemistrySh = kappa*chemistry.Sh()();
+
         #include "rhoEqn.H"
         #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
diff --git a/applications/solvers/combustion/dieselEngineFoam/hEqn.H b/applications/solvers/combustion/dieselEngineFoam/hEqn.H
deleted file mode 100644
index 0406f7fbd4..0000000000
--- a/applications/solvers/combustion/dieselEngineFoam/hEqn.H
+++ /dev/null
@@ -1,13 +0,0 @@
-{
-    solve
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-       DpDt
-     + dieselSpray.heatTransferSource()
-    );
-
-    thermo.correct();
-}
diff --git a/applications/solvers/combustion/dieselEngineFoam/hsEqn.H b/applications/solvers/combustion/dieselEngineFoam/hsEqn.H
new file mode 100644
index 0000000000..7ae59feb81
--- /dev/null
+++ b/applications/solvers/combustion/dieselEngineFoam/hsEqn.H
@@ -0,0 +1,14 @@
+{
+    solve
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+       DpDt
+     + dieselSpray.heatTransferSource()().dimensionedInternalField()
+     + chemistrySh
+    );
+
+    thermo.correct();
+}
diff --git a/applications/solvers/combustion/dieselFoam/dieselFoam.C b/applications/solvers/combustion/dieselFoam/dieselFoam.C
index 4769a00c77..1e6ddfe0c1 100644
--- a/applications/solvers/combustion/dieselFoam/dieselFoam.C
+++ b/applications/solvers/combustion/dieselFoam/dieselFoam.C
@@ -100,7 +100,7 @@ int main(int argc, char *argv[])
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
diff --git a/applications/solvers/combustion/fireFoam/Allwclean b/applications/solvers/combustion/fireFoam/Allwclean
index 5c39709363..713a055201 100755
--- a/applications/solvers/combustion/fireFoam/Allwclean
+++ b/applications/solvers/combustion/fireFoam/Allwclean
@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wclean libso combustionModels
-wclean all sensibleEnthalpyCombustionThermophysicalModels
 wclean
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/combustion/fireFoam/Allwmake b/applications/solvers/combustion/fireFoam/Allwmake
index ada80d6666..f32c72860d 100755
--- a/applications/solvers/combustion/fireFoam/Allwmake
+++ b/applications/solvers/combustion/fireFoam/Allwmake
@@ -2,7 +2,6 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-sensibleEnthalpyCombustionThermophysicalModels/Allwmake
 wmake libso combustionModels
 wmake
 
diff --git a/applications/solvers/combustion/fireFoam/Make/options b/applications/solvers/combustion/fireFoam/Make/options
index 07d596f7f0..3b19a78173 100644
--- a/applications/solvers/combustion/fireFoam/Make/options
+++ b/applications/solvers/combustion/fireFoam/Make/options
@@ -1,7 +1,5 @@
 EXE_INC = \
     -I./combustionModels/lnInclude \
-    -I./sensibleEnthalpyCombustionThermophysicalModels/basic/lnInclude \
-    -I./sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
@@ -11,8 +9,6 @@ EXE_INC = \
 
 
 EXE_LIBS = \
-    -lbasicSensibleEnthalpyThermophysicalModels \
-    -lreactionSensibleEnthalpyThermophysicalModels \
     -lcombustionModels \
     -lspecie \
     -lreactionThermophysicalModels \
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/Allwmake b/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/Allwmake
deleted file mode 100755
index e59cf84715..0000000000
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/Allwmake
+++ /dev/null
@@ -1,8 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # run from this directory
-set -x
-
-wmake libso basic
-wmake libso reactionThermo
-
-# ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/Make/files b/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/Make/files
deleted file mode 100644
index ca1ec78b98..0000000000
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/Make/files
+++ /dev/null
@@ -1,4 +0,0 @@
-psiThermo/basicSensiblePsiThermo/basicSensiblePsiThermo.C
-psiThermo/basicSensiblePsiThermo/newBasicSensiblePsiThermo.C
-
-LIB = $(FOAM_LIBBIN)/libbasicSensibleEnthalpyThermophysicalModels
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/Make/options b/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/Make/options
deleted file mode 100644
index 8bdd75aecb..0000000000
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/Make/options
+++ /dev/null
@@ -1,7 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
-
-LIB_LIBS = \
-    -lfiniteVolume
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/makeBasicSensiblePsiThermo.H b/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/makeBasicSensiblePsiThermo.H
deleted file mode 100644
index 7d0a078719..0000000000
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/makeBasicSensiblePsiThermo.H
+++ /dev/null
@@ -1,66 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-InClass
-    Foam::basicSensiblePsiThermo
-
-Description
-    Macros for creating 'basic' compresibility-based thermo packages
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef makeBasicSensiblePsiThermo_H
-#define makeBasicSensiblePsiThermo_H
-
-#include "basicSensiblePsiThermo.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#define makeBasicSensiblePsiThermo(Cthermo,Mixture,Transport,Thermo,EqnOfState)       \
-                                                                              \
-typedef Cthermo<Mixture<Transport<specieThermo<Thermo<EqnOfState> > > > >     \
-    Cthermo##Mixture##Transport##Thermo##EqnOfState;                          \
-                                                                              \
-defineTemplateTypeNameAndDebugWithName                                        \
-(                                                                             \
-    Cthermo##Mixture##Transport##Thermo##EqnOfState,                          \
-    #Cthermo                                                                  \
-        "<"#Mixture"<"#Transport"<specieThermo<"#Thermo"<"#EqnOfState">>>>>", \
-    0                                                                         \
-);                                                                            \
-                                                                              \
-addToRunTimeSelectionTable                                                    \
-(                                                                             \
-    basicSensiblePsiThermo,                                                           \
-    Cthermo##Mixture##Transport##Thermo##EqnOfState,                          \
-    fvMesh                                                                    \
-)
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/newBasicSensiblePsiThermo.C b/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/newBasicSensiblePsiThermo.C
deleted file mode 100644
index 14fcbbff6a..0000000000
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/newBasicSensiblePsiThermo.C
+++ /dev/null
@@ -1,75 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "basicSensiblePsiThermo.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-Foam::autoPtr<Foam::basicSensiblePsiThermo> Foam::basicSensiblePsiThermo::New
-(
-    const fvMesh& mesh
-)
-{
-    word thermoTypeName;
-
-    // Enclose the creation of the thermophysicalProperties to ensure it is
-    // deleted before the turbulenceModel is created otherwise the dictionary
-    // is entered in the database twice
-    {
-        IOdictionary thermoDict
-        (
-            IOobject
-            (
-                "thermophysicalProperties",
-                mesh.time().constant(),
-                mesh,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE
-            )
-        );
-
-        thermoDict.lookup("thermoType") >> thermoTypeName;
-    }
-
-    Info<< "Selecting thermodynamics package " << thermoTypeName << endl;
-
-    fvMeshConstructorTable::iterator cstrIter =
-        fvMeshConstructorTablePtr_->find(thermoTypeName);
-
-    if (cstrIter == fvMeshConstructorTablePtr_->end())
-    {
-        FatalErrorIn("basicSensiblePsiThermo::New(const fvMesh&)")
-            << "Unknown basicSensiblePsiThermo type " << thermoTypeName << nl << nl
-            << "Valid basicSensiblePsiThermo types are:" << nl
-            << fvMeshConstructorTablePtr_->toc() << nl
-            << exit(FatalError);
-    }
-
-    return autoPtr<basicSensiblePsiThermo>(cstrIter()(mesh));
-}
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/Make/files b/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/Make/files
deleted file mode 100644
index 39b2fa53b1..0000000000
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/Make/files
+++ /dev/null
@@ -1,5 +0,0 @@
-combustionThermo/hsCombustionThermo/hsCombustionThermo.C
-combustionThermo/hsCombustionThermo/newhsCombustionThermo.C
-combustionThermo/hsCombustionThermo/hsCombustionThermos.C
-
-LIB = $(FOAM_LIBBIN)/libreactionSensibleEnthalpyThermophysicalModels
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/Make/options b/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/Make/options
deleted file mode 100644
index 9551ad0315..0000000000
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/Make/options
+++ /dev/null
@@ -1,9 +0,0 @@
-EXE_INC = \
-    -I../basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-LIB_LIBS = \
-    -lfiniteVolume
diff --git a/applications/solvers/combustion/reactingFoam/Make/options b/applications/solvers/combustion/reactingFoam/Make/options
index e982f635a1..d85350297a 100644
--- a/applications/solvers/combustion/reactingFoam/Make/options
+++ b/applications/solvers/combustion/reactingFoam/Make/options
@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I../XiFoam \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
diff --git a/applications/solvers/combustion/reactingFoam/UEqn.H b/applications/solvers/combustion/reactingFoam/UEqn.H
new file mode 100644
index 0000000000..9697c6e1ed
--- /dev/null
+++ b/applications/solvers/combustion/reactingFoam/UEqn.H
@@ -0,0 +1,15 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho*g
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }
diff --git a/applications/solvers/combustion/reactingFoam/chemistry.H b/applications/solvers/combustion/reactingFoam/chemistry.H
index a1a978210d..8bfa872d7b 100644
--- a/applications/solvers/combustion/reactingFoam/chemistry.H
+++ b/applications/solvers/combustion/reactingFoam/chemistry.H
@@ -21,4 +21,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }
diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H
index 0b2238031c..bd8f7ea1fb 100644
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/createFields.H
@@ -5,7 +5,7 @@ autoPtr<psiChemistryModel> pChemistry
 );
 psiChemistryModel& chemistry = pChemistry();
 
-hCombustionThermo& thermo = chemistry.thermo();
+hsCombustionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -40,8 +40,8 @@ volVectorField U
 
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
-volScalarField& h = thermo.h();
-
+volScalarField& hs = thermo.hs();
+const volScalarField& T = thermo.T();
 
 #include "compressibleCreatePhi.H"
 
@@ -81,5 +81,18 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
 
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);
diff --git a/applications/solvers/combustion/reactingFoam/hsEqn.H b/applications/solvers/combustion/reactingFoam/hsEqn.H
new file mode 100644
index 0000000000..e3fa4e7a13
--- /dev/null
+++ b/applications/solvers/combustion/reactingFoam/hsEqn.H
@@ -0,0 +1,20 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+//      - fvm::laplacian(turbulence->muEff(), hs)  // unit lewis no.
+     ==
+        DpDt
+      + chemistrySh
+    );
+
+    hsEqn.relax();
+    hsEqn.solve();
+
+    thermo.correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
new file mode 100644
index 0000000000..9443f909a3
--- /dev/null
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -0,0 +1,68 @@
+rho = thermo.rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+U = rUA*UEqn.H();
+
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    phi =
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index 2d4ae7589a..e9fe757e3b 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C
@@ -73,9 +73,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-
-            #define Db turbulence->alphaEff()
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
diff --git a/applications/solvers/combustion/rhoReactingFoam/Make/options b/applications/solvers/combustion/rhoReactingFoam/Make/options
index e982f635a1..d85350297a 100644
--- a/applications/solvers/combustion/rhoReactingFoam/Make/options
+++ b/applications/solvers/combustion/rhoReactingFoam/Make/options
@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I../XiFoam \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
diff --git a/applications/solvers/combustion/rhoReactingFoam/UEqn.H b/applications/solvers/combustion/rhoReactingFoam/UEqn.H
new file mode 100644
index 0000000000..9697c6e1ed
--- /dev/null
+++ b/applications/solvers/combustion/rhoReactingFoam/UEqn.H
@@ -0,0 +1,15 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho*g
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }
diff --git a/applications/solvers/combustion/rhoReactingFoam/chemistry.H b/applications/solvers/combustion/rhoReactingFoam/chemistry.H
index a1a978210d..8bfa872d7b 100644
--- a/applications/solvers/combustion/rhoReactingFoam/chemistry.H
+++ b/applications/solvers/combustion/rhoReactingFoam/chemistry.H
@@ -21,4 +21,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }
diff --git a/applications/solvers/combustion/rhoReactingFoam/createFields.H b/applications/solvers/combustion/rhoReactingFoam/createFields.H
index d2f6f5e4d4..d4119b8c92 100644
--- a/applications/solvers/combustion/rhoReactingFoam/createFields.H
+++ b/applications/solvers/combustion/rhoReactingFoam/createFields.H
@@ -5,7 +5,7 @@ autoPtr<rhoChemistryModel> pChemistry
 );
 rhoChemistryModel& chemistry = pChemistry();
 
-hReactionThermo& thermo = chemistry.thermo();
+hsReactionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -40,7 +40,8 @@ volVectorField U
 
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
-volScalarField& h = thermo.h();
+volScalarField& hs = thermo.hs();
+const volScalarField& T = thermo.T();
 
 
 #include "compressibleCreatePhi.H"
@@ -81,5 +82,18 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
 
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);
diff --git a/applications/solvers/combustion/rhoReactingFoam/hsEqn.H b/applications/solvers/combustion/rhoReactingFoam/hsEqn.H
new file mode 100644
index 0000000000..81bd6ea9d7
--- /dev/null
+++ b/applications/solvers/combustion/rhoReactingFoam/hsEqn.H
@@ -0,0 +1,19 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+      + chemistrySh
+    );
+
+    hsEqn.relax();
+    hsEqn.solve();
+
+    thermo.correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
diff --git a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
index cc37dd09c1..79b76e014d 100644
--- a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
@@ -74,7 +74,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/chemistry.H b/applications/solvers/lagrangian/coalChemistryFoam/chemistry.H
index 5489dfefba..50fc7f575b 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/chemistry.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/chemistry.H
@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
index 389a065935..872664f3a5 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@@ -84,16 +84,8 @@ int main(int argc, char *argv[])
 
         coalParcels.evolve();
 
-        coalParcels.info();
-
-        Info<< endl;
-
         limestoneParcels.evolve();
 
-        limestoneParcels.info();
-
-        Info<< endl;
-
         #include "chemistry.H"
         #include "rhoEqn.H"
 
@@ -102,16 +94,13 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
                 #include "pEqn.H"
             }
-
-            Info<< "T gas min/max   = " << min(T).value() << ", "
-                << max(T).value() << endl;
         }
 
         turbulence->correct();
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
index d0a284a0e1..65be09775c 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
@@ -6,7 +6,7 @@
     );
     psiChemistryModel& chemistry = pChemistry();
 
-    hCombustionThermo& thermo = chemistry.thermo();
+    hsCombustionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,7 +22,7 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -32,7 +32,7 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
 
     volScalarField rho
     (
@@ -133,5 +133,19 @@
         "energy",
         mesh,
         dimEnergy/dimTime/dimVolume,
-        "h"
+        "hs"
+    );
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
     );
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/hEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/hEqn.H
deleted file mode 100644
index 2eaf6315cf..0000000000
--- a/applications/solvers/lagrangian/coalChemistryFoam/hEqn.H
+++ /dev/null
@@ -1,22 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        DpDt
-     +  coalParcels.Sh()
-     +  limestoneParcels.Sh()
-     +  enthalpySource.Su()
-     +  radiation->Sh(thermo)
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-}
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
new file mode 100644
index 0000000000..ced6b64004
--- /dev/null
+++ b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
@@ -0,0 +1,26 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+     +  coalParcels.Sh()
+     +  limestoneParcels.Sh()
+     +  enthalpySource.Su()
+     +  radiation->Shs(thermo)
+     +  chemistrySh
+    );
+
+    hsEqn.relax();
+
+    hsEqn.solve();
+
+    thermo.correct();
+
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H
index 5489dfefba..50fc7f575b 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H
@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
index 68c0f8aee8..cf3f484bdb 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
@@ -6,7 +6,7 @@
     );
     rhoChemistryModel& chemistry = pChemistry();
 
-    hReactionThermo& thermo = chemistry.thermo();
+    hsReactionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,7 +22,7 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -88,4 +88,18 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+    );
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
similarity index 83%
rename from applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hEqn.H
rename to applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
index 0715989df7..73c276c294 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
@@ -32,14 +32,15 @@
     {
         solve
         (
-            fvm::ddt(rho, h)
-          + mvConvection->fvmDiv(phi, h)
-          - fvm::laplacian(turbulence->alphaEff(), h)
+            fvm::ddt(rho, hs)
+          + mvConvection->fvmDiv(phi, hs)
+          - fvm::laplacian(turbulence->alphaEff(), hs)
          ==
             pWork()
           + parcels.Sh()
-          + radiation->Sh(thermo)
+          + radiation->Shs(thermo)
           + energySource.Su()
+          + chemistrySh
         );
 
         thermo.correct();
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
index b345ffe28a..436fda7a8b 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
@@ -33,7 +33,7 @@ Description
     The solver includes:
     - reacting multiphase parcel cloud
     - porous media
-    - point mass sources
+    - mass, momentum and energy sources
     - polynomial based, incompressible thermodynamics (f(T))
 
     Note: ddtPhiCorr not used here when porous zones are active
@@ -89,13 +89,11 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        parcels.info();
-
         #include "chemistry.H"
         #include "rhoEqn.H"
         #include "UEqn.H"
         #include "YEqn.H"
-        #include "hEqn.H"
+        #include "hsEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/reactingParcelFoam/chemistry.H
index 5489dfefba..50fc7f575b 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/chemistry.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/chemistry.H
@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
index e825c2cffb..e071ac8b3d 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
@@ -6,7 +6,7 @@
     );
     psiChemistryModel& chemistry = pChemistry();
 
-    hCombustionThermo& thermo = chemistry.thermo();
+    hsCombustionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,7 +22,7 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -94,4 +94,18 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistry::Sh", dimEnergy/dimTime/dimVolume, 0.0)
+    );
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/hEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/hEqn.H
deleted file mode 100644
index d686e452f4..0000000000
--- a/applications/solvers/lagrangian/reactingParcelFoam/hEqn.H
+++ /dev/null
@@ -1,20 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        DpDt
-     +  parcels.Sh()
-     +  radiation->Sh(thermo)
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-}
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
new file mode 100644
index 0000000000..ce12ec3ad4
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
@@ -0,0 +1,24 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+     +  parcels.Sh()
+     +  radiation->Shs(thermo)
+     +  chemistrySh
+    );
+
+    hEqn.relax();
+
+    hEqn.solve();
+
+    thermo.correct();
+
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index 4dfd367d09..ad3b5b461c 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -74,8 +74,6 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        parcels.info();
-
         #include "chemistry.H"
         #include "rhoEqn.H"
 
@@ -88,12 +86,9 @@ int main(int argc, char *argv[])
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-                #include "hEqn.H"
+                #include "hsEqn.H"
                 #include "pEqn.H"
             }
-
-            Info<< "T gas min/max   = " << min(T).value() << ", "
-                << max(T).value() << endl;
         }
 
         turbulence->correct();
diff --git a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C
index ee0088c397..ef1ccbd150 100644
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C
@@ -61,7 +61,6 @@ int main(int argc, char *argv[])
 
         Info<< "Evolving " << kinematicCloud.name() << endl;
         kinematicCloud.evolve();
-        kinematicCloud.info();
 
         runTime.write();
 
diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
index d105a8492a..d2dcd818ca 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
@@ -270,11 +270,31 @@ int main(int argc, char *argv[])
                 procMeshes.boundaryProcAddressing()
             );
 
-            pointReconstructor.reconstructFields<scalar>(objects);
-            pointReconstructor.reconstructFields<vector>(objects);
-            pointReconstructor.reconstructFields<sphericalTensor>(objects);
-            pointReconstructor.reconstructFields<symmTensor>(objects);
-            pointReconstructor.reconstructFields<tensor>(objects);
+            pointReconstructor.reconstructFields<scalar>
+            (
+                objects,
+                selectedFields
+            );
+            pointReconstructor.reconstructFields<vector>
+            (
+                objects,
+                selectedFields
+            );
+            pointReconstructor.reconstructFields<sphericalTensor>
+            (
+                objects,
+                selectedFields
+            );
+            pointReconstructor.reconstructFields<symmTensor>
+            (
+                objects,
+                selectedFields
+            );
+            pointReconstructor.reconstructFields<tensor>
+            (
+                objects,
+                selectedFields
+            );
         }
         else
         {
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
index 48435c749b..19920b008a 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
@@ -30,13 +30,10 @@ License
 #include "fvMesh.H"
 #include "globalMeshData.H"
 #include "PstreamCombineReduceOps.H"
-#include "processorPolyPatch.H"
 #include "cellModeller.H"
 #include "IOmanip.H"
 #include "itoa.H"
 #include "ensightWriteBinary.H"
-#include "globalIndex.H"
-#include "PackedBoolList.H"
 #include "mapDistribute.H"
 
 #include <fstream>
@@ -1111,14 +1108,11 @@ void Foam::ensightMesh::writeAllFacePrimsBinary
             writeEnsDataBinary(key,ensightGeometryFile);
             writeEnsDataBinary(nPrims,ensightGeometryFile);
 
-            if (&prims != NULL)
-            {
-                writeFacePrimsBinary
-                (
-                    UIndirectList<face>(patchFaces, prims)(),
-                    ensightGeometryFile
-                );
-            }
+            writeFacePrimsBinary
+            (
+                UIndirectList<face>(patchFaces, prims)(),
+                ensightGeometryFile
+            );
 
             for (int slave=1; slave<Pstream::nProcs(); slave++)
             {
@@ -1348,17 +1342,18 @@ void Foam::ensightMesh::writeAscii
 
                 // Renumber the patch points/faces into unique points
                 labelList pointToGlobal;
-                labelList uniquePointLabels;
+                labelList uniqueMeshPointLabels;
                 autoPtr<globalIndex> globalPointsPtr =
                 mesh_.globalData().mergePoints
                 (
                     p.meshPoints(),
                     p.meshPointMap(),
                     pointToGlobal,
-                    uniquePointLabels
+                    uniqueMeshPointLabels
                 );
 
-                pointField uniquePoints(p.localPoints(), uniquePointLabels);
+
+                pointField uniquePoints(mesh_.points(), uniqueMeshPointLabels);
                 // Renumber the patch faces
                 faceList patchFaces(p.localFaces());
                 forAll(patchFaces, i)
@@ -1601,16 +1596,17 @@ void Foam::ensightMesh::writeBinary
 
                 // Renumber the patch points/faces into unique points
                 labelList pointToGlobal;
-                labelList uniquePointLabels;
+                labelList uniqueMeshPointLabels;
                 autoPtr<globalIndex> globalPointsPtr =
                 mesh_.globalData().mergePoints
                 (
                     p.meshPoints(),
                     p.meshPointMap(),
                     pointToGlobal,
-                    uniquePointLabels
+                    uniqueMeshPointLabels
                 );
-                pointField uniquePoints(p.localPoints(), uniquePointLabels);
+
+                pointField uniquePoints(mesh_.points(), uniqueMeshPointLabels);
                 // Renumber the patch faces
                 faceList patchFaces(p.localFaces());
                 forAll(patchFaces, i)
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
index c3650ccb5a..6f9f311d11 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
@@ -1705,7 +1705,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
     const labelList& meshPoints,
     const Map<label>& meshPointMap,
     labelList& pointToGlobal,
-    labelList& uniquePoints
+    labelList& uniqueMeshPoints
 ) const
 {
     const indirectPrimitivePatch& cpp = coupledPatch();
@@ -1713,125 +1713,175 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
     const mapDistribute& pointSlavesMap = globalPointSlavesMap();
 
 
-    // 1. Count number of masters on my processor.
-    label nCoupledMaster = 0;
-    label nCoupledSlave = 0;
-    PackedBoolList isMaster(meshPoints.size(), 1);
+    // The patch points come in two variants:
+    // - not on a coupled patch so guaranteed unique
+    // - on a coupled patch
+    // If the point is on a coupled patch the problem is that the
+    // master-slave structure (globalPointSlaves etc.) assigns one of the
+    // coupled points to be the master but this master point is not
+    // necessarily on the patch itself! (it might just be connected to the
+    // patch point via coupled patches).
+    // So this means that all master point loops should be over the
+    // master-slave structure, not over the patch points and that the unique
+    // point returned is a mesh point.
+    // (unless we want to do the whole master-slave analysis again for the
+    //  current patch only).
 
-    forAll(meshPoints, localPointI)
+
+    // Determine mapping from coupled point to patch point and vice versa
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList patchToCoupled(meshPoints.size(), -1);
+    label nCoupled = 0;
+    labelList coupledToPatch(pointSlavesMap.constructSize(), -1);
+
+    // Note: loop over patch since usually smaller
+    forAll(meshPoints, patchPointI)
     {
-        label meshPointI = meshPoints[localPointI];
+        label meshPointI = meshPoints[patchPointI];
 
         Map<label>::const_iterator iter = cpp.meshPointMap().find(meshPointI);
 
         if (iter != cpp.meshPointMap().end())
         {
-            // My localPointI is a coupled point.
+            patchToCoupled[patchPointI] = iter();
+            coupledToPatch[iter()] = patchPointI;
+            nCoupled++;
+        }
+    }
 
-            label coupledPointI = iter();
+    //Pout<< "Patch:" << nl
+    //    << "    points:" << meshPoints.size() << nl
+    //    << "    of which on coupled patch:" << nCoupled << endl;
 
-            if (pointSlaves[coupledPointI].size() > 0)
+
+    // Pull coupled-to-patch information to master
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    pointSlavesMap.distribute(coupledToPatch);
+
+
+    // Check on master whether point is anywhere on patch
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // List of master points that are on the patch
+    DynamicList<label> masterPoints(pointSlaves.size());
+
+    forAll(pointSlaves, coupledPointI)
+    {
+        const labelList& slaves = pointSlaves[coupledPointI];
+
+        if (slaves.size() > 0)
+        {
+            // I am master. Is this point on the patch on myself or on any
+            // any slave?
+            if (coupledToPatch[coupledPointI] != -1)
             {
-                nCoupledMaster++;
+                masterPoints.append(coupledPointI);
             }
             else
             {
-                isMaster[localPointI] = 0;
-                nCoupledSlave++;
-            }
-        }
-    }
-
-    label myUniquePoints = meshPoints.size() - nCoupledSlave;
-
-    Pout<< "Points :" << nl
-        << "    patch            : " << meshPoints.size() << nl
-        << "    of which coupled : " << nCoupledMaster+nCoupledSlave << nl
-        << "    of which master  : " << nCoupledMaster << nl
-        << "    of which slave   : " << nCoupledSlave << nl
-        << endl;
-
-
-    // 2. Create global indexing for unique points.
-    autoPtr<globalIndex> globalPointsPtr(new globalIndex(myUniquePoints));
-
-
-    // 3. Assign global point numbers. Keep slaves unset.
-    pointToGlobal.setSize(meshPoints.size());
-    pointToGlobal = -1;
-    uniquePoints.setSize(myUniquePoints);
-    label nMaster = 0;
-
-    forAll(isMaster, localPointI)
-    {
-        if (isMaster[localPointI])
-        {
-            pointToGlobal[localPointI] = globalPointsPtr().toGlobal(nMaster);
-            uniquePoints[nMaster] = localPointI;
-            nMaster++;
-        }
-    }
-
-
-    // 4. Push global index for coupled points to slaves.
-    {
-        labelList masterToGlobal(pointSlavesMap.constructSize(), -1);
-
-        forAll(meshPoints, localPointI)
-        {
-            label meshPointI = meshPoints[localPointI];
-
-            Map<label>::const_iterator iter = cpp.meshPointMap().find
-            (
-                meshPointI
-            );
-
-            if (iter != cpp.meshPointMap().end())
-            {
-                // My localPointI is a coupled point.
-                label coupledPointI = iter();
-
-                const labelList& slaves = pointSlaves[coupledPointI];
-
-                if (slaves.size() > 0)
+                forAll(slaves, i)
                 {
-                    // Duplicate master globalpoint into slave slots
-                    masterToGlobal[coupledPointI] = pointToGlobal[localPointI];
-                    forAll(slaves, i)
+                    if (coupledToPatch[slaves[i]] != -1)
                     {
-                        masterToGlobal[slaves[i]] = pointToGlobal[localPointI];
+                        masterPoints.append(coupledPointI);
+                        break;
                     }
                 }
             }
         }
+    }
 
-        // Send back
-        pointSlavesMap.reverseDistribute(cpp.nPoints(), masterToGlobal);
 
-        // On slave copy master index into overal map.
-        forAll(meshPoints, localPointI)
+    // Create global indexing
+    // ~~~~~~~~~~~~~~~~~~~~~~
+    // 1. patch points that are not on coupled patch:
+    //      meshPoints.size()-nCoupled
+    // 2. master points that are on patch:
+    //      masterPoints.size()
+    label myUniquePoints = meshPoints.size()-nCoupled+masterPoints.size();
+    autoPtr<globalIndex> globalPointsPtr(new globalIndex(myUniquePoints));
+
+    //Pout<< "CoupledPatch:" << nl
+    //    << "    points:" << cpp.nPoints() << nl
+    //    << "    of which on patch:" << masterPoints.size() << endl;
+
+
+    // Allocate unique points
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    pointToGlobal.setSize(meshPoints.size());
+    pointToGlobal = -1;
+    uniqueMeshPoints.setSize(myUniquePoints);
+
+    // Allocate globals for uncoupled patch points
+    label nMaster = 0;
+    forAll(patchToCoupled, patchPointI)
+    {
+        if (patchToCoupled[patchPointI] == -1)
         {
-            label meshPointI = meshPoints[localPointI];
+            // Allocate global point
+            label globalI = globalPointsPtr().toGlobal(nMaster);
+            pointToGlobal[patchPointI] = globalI;
+            uniqueMeshPoints[nMaster] = meshPoints[patchPointI];
+            nMaster++;
+        }
+    }
 
-            Map<label>::const_iterator iter = cpp.meshPointMap().find
-            (
-                meshPointI
-            );
+    // Allocate globals for master 
+    labelList masterToGlobal(pointSlavesMap.constructSize(), -456);
 
-            if (iter != cpp.meshPointMap().end())
+    forAll(masterPoints, i)
+    {
+        label coupledPointI = masterPoints[i];
+
+        // Allocate global point
+        label globalI = globalPointsPtr().toGlobal(nMaster);
+        if (coupledToPatch[coupledPointI] != -1)
+        {
+            pointToGlobal[coupledToPatch[coupledPointI]] = globalI;
+        }
+        uniqueMeshPoints[nMaster] = cpp.meshPoints()[coupledPointI];
+        nMaster++;
+
+        // Put global into slave slots
+        const labelList& slaves = pointSlaves[coupledPointI];
+        masterToGlobal[coupledPointI] = globalI;    // not really necessary
+        forAll(slaves, i)
+        {
+            masterToGlobal[slaves[i]] = globalI;
+        }
+    }
+
+
+    if (nMaster != myUniquePoints)
+    {
+        FatalErrorIn("globalMeshData::mergePoints(..)")
+            << "problem." << abort(FatalError);
+    }
+
+
+    // Send back (from slave slots) to originating processor
+    pointSlavesMap.reverseDistribute(cpp.nPoints(), masterToGlobal);
+
+    // On slaves take over global number
+    forAll(patchToCoupled, patchPointI)
+    {
+        label coupledPointI = patchToCoupled[patchPointI];
+
+        if (coupledPointI != -1)
+        {
+            const labelList& slaves = pointSlaves[coupledPointI];
+
+            if (slaves.size() == 0)
             {
-                // My localPointI is a coupled point.
-                label coupledPointI = iter();
-                const labelList& slaves = pointSlaves[coupledPointI];
-
-                if (slaves.size() == 0)
-                {
-                    pointToGlobal[localPointI] = masterToGlobal[coupledPointI];
-                }
+                pointToGlobal[patchPointI] = masterToGlobal[coupledPointI];
             }
         }
     }
 
+
     return globalPointsPtr;
 }
 
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
index f397d69ae1..a20caaa010 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
@@ -560,13 +560,17 @@ public:
                 ) const;
 
                 //- Helper for merging patch point data. Takes maps from
-                //  local points to/from mesh
+                //  local points to/from mesh. Determines
+                //  - my unique points. These are mesh points, not patch points
+                //    since the master might not be on the patch.
+                //  - global numbering over all unique indices.
+                //  - the global number for all local points.
                 autoPtr<globalIndex> mergePoints
                 (
                     const labelList& meshPoints,
                     const Map<label>& meshPointMap,
                     labelList& pointToGlobal,
-                    labelList& uniquePoints
+                    labelList& uniqueMeshPoints
                 ) const;
 
 
diff --git a/src/lagrangian/dieselSpray/parcel/parcel.C b/src/lagrangian/dieselSpray/parcel/parcel.C
index c091325845..5c555e1940 100644
--- a/src/lagrangian/dieselSpray/parcel/parcel.C
+++ b/src/lagrangian/dieselSpray/parcel/parcel.C
@@ -247,7 +247,7 @@ bool Foam::parcel::move(spray& sDB)
         {
             oMass[i] = m()*oYf[i];
             label j = sDB.liquidToGasIndex()[i];
-            oHg += oYf[i]*sDB.gasProperties()[j].H(T());
+            oHg += oYf[i]*sDB.gasProperties()[j].Hs(T());
         }
 
         vector oMom = m()*U();
@@ -273,7 +273,7 @@ bool Foam::parcel::move(spray& sDB)
         {
             nMass[i] = m()*nYf[i];
             label j = sDB.liquidToGasIndex()[i];
-            nHg += nYf[i]*sDB.gasProperties()[j].H(T());
+            nHg += nYf[i]*sDB.gasProperties()[j].Hs(T());
         }
 
         vector nMom = m()*U();
@@ -286,11 +286,9 @@ bool Foam::parcel::move(spray& sDB)
         {
             sDB.srhos()[i][celli] += oMass[i] - nMass[i];
         }
-        sDB.sms()[celli]   += oMom - nMom;
+        sDB.sms()[celli] += oMom - nMom;
 
-        sDB.shs()[celli]   +=
-            oTotMass*(oH + oPE)
-          - m()*(nH + nPE);
+        sDB.shs()[celli] += oTotMass*(oH + oPE) - m()*(nH + nPE);
 
         // Remove evaporated mass from stripped mass
         ms() -= ms()*(oTotMass-m())/oTotMass;
@@ -446,7 +444,7 @@ void Foam::parcel::updateParcelProperties
     for (label i=0; i<Nf; i++)
     {
         label j = sDB.liquidToGasIndex()[i];
-        oldhg += Yf0[i]*sDB.gasProperties()[j].H(T());
+        oldhg += Yf0[i]*sDB.gasProperties()[j].Hs(T());
     }
 
     scalar oldhv = fuels.hl(pg, T(), X());
@@ -478,7 +476,7 @@ void Foam::parcel::updateParcelProperties
             for (label i=0; i<Nf; i++)
             {
                 label j = sDB.liquidToGasIndex()[i];
-                newhg += Ynew[i]*sDB.gasProperties()[j].H(Tnew);
+                newhg += Ynew[i]*sDB.gasProperties()[j].Hs(Tnew);
             }
 
             newhv = fuels.hl(pg, Tnew, X());
@@ -616,7 +614,8 @@ void Foam::parcel::updateParcelProperties
                     {
                         if (n>100)
                         {
-                            Info<< "n = " << n << ", T = " << Td << ", pv = " << pAtSurface << endl;
+                            Info<< "n = " << n << ", T = " << Td << ", pv = "
+                                << pAtSurface << endl;
                         }
                     }
                 }
diff --git a/src/lagrangian/dieselSpray/parcel/parcel.H b/src/lagrangian/dieselSpray/parcel/parcel.H
index c6084d08ee..b917802927 100644
--- a/src/lagrangian/dieselSpray/parcel/parcel.H
+++ b/src/lagrangian/dieselSpray/parcel/parcel.H
@@ -185,7 +185,8 @@ public:
             //- Return the names of the liquid components
             inline const List<word>& liquidNames() const;
 
-            //- Return the names of the liquid fuel components - identical with liquidNames
+            //- Return the names of the liquid fuel components
+            //  - identical with liquidNames
             inline const List<word>& fuelNames() const;
 
             //- Return diameter of droplets in parcel
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index f0c3c5f033..29130c69dc 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -260,6 +260,9 @@ void Foam::KinematicCloud<ParcelType>::evolve()
         evolveCloud();
 
         postEvolve();
+
+        info();
+        Info<< endl;
     }
 }
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index 210a589902..f0ccd58f19 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -272,6 +272,9 @@ void Foam::ReactingCloud<ParcelType>::evolve()
         evolveCloud();
 
         postEvolve();
+
+        info();
+        Info<< endl;
     }
 }
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
index b783cb2b08..3df450724b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
@@ -228,6 +228,9 @@ void Foam::ReactingMultiphaseCloud<ParcelType>::evolve()
         evolveCloud();
 
         postEvolve();
+
+        info();
+        Info<< endl;
     }
 }
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index 0d12d286b9..b3135ce81c 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -160,20 +160,6 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
         ),
         this->mesh(),
         dimensionedScalar("zero", dimEnergy, 0.0)
-    ),
-    hcTrans_
-    (
-        IOobject
-        (
-            this->name() + "hcTrans",
-            this->db().time().timeName(),
-            this->db(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        this->mesh(),
-        dimensionedScalar("zero", dimEnergy, 0.0)
     )
 {
     if (readFields)
@@ -220,7 +206,6 @@ void Foam::ThermoCloud<ParcelType>::resetSourceTerms()
 {
     KinematicCloud<ParcelType>::resetSourceTerms();
     hsTrans_.field() = 0.0;
-    hcTrans_.field() = 0.0;
 }
 
 
@@ -234,6 +219,9 @@ void Foam::ThermoCloud<ParcelType>::evolve()
         evolveCloud();
 
         postEvolve();
+
+        info();
+        Info<< endl;
     }
 }
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
index c0e2d5864a..b5d9f27b8b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
@@ -108,15 +108,9 @@ protected:
 
         // Sources
 
-            //- Sensible enthalpy transfer
+            //- Sensible enthalpy transfer [J/kg]
             DimensionedField<scalar, volMesh> hsTrans_;
 
-            //- Chemical enthalpy transfer
-            // - If solving for total enthalpy, the carrier phase enthalpy will
-            //   receive the full enthalpy of reaction via creation of reaction
-            //   products
-            DimensionedField<scalar, volMesh> hcTrans_;
-
 
     // Protected member functions
 
@@ -194,19 +188,10 @@ public:
 
                 // Enthalpy
 
-                    //- Return reference to sensible enthalpy source
+                    //- Sensible enthalpy transfer [J/kg]
                     inline DimensionedField<scalar, volMesh>& hsTrans();
 
-                    //- Return tmp total sensible enthalpy source term
-                    inline tmp<DimensionedField<scalar, volMesh> > Shs() const;
-
-                    //- Return reference to chemical enthalpy source
-                    inline DimensionedField<scalar, volMesh>& hcTrans();
-
-                    //- Return tmp chemical enthalpy source term
-                    inline tmp<DimensionedField<scalar, volMesh> > Shc() const;
-
-                    //- Return tmp total enthalpy source term
+                    //- Return enthalpy source [J/kg/m3/s]
                     inline tmp<DimensionedField<scalar, volMesh> > Sh() const;
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
index f537f6283a..249906311d 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
@@ -85,82 +85,6 @@ Foam::ThermoCloud<ParcelType>::hsTrans()
 }
 
 
-template<class ParcelType>
-inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ThermoCloud<ParcelType>::Shs() const
-{
-    tmp<DimensionedField<scalar, volMesh> > tShs
-    (
-        new DimensionedField<scalar, volMesh>
-        (
-            IOobject
-            (
-                this->name() + "Shs",
-                this->db().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            this->mesh(),
-            dimensionedScalar
-            (
-                "zero",
-                dimMass/dimLength/pow3(dimTime),
-                0.0
-            )
-        )
-    );
-
-    scalarField& Shs = tShs().field();
-    Shs = hsTrans_/(this->mesh().V()*this->db().time().deltaT());
-
-    return tShs;
-}
-
-
-template<class ParcelType>
-inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::ThermoCloud<ParcelType>::hcTrans()
-{
-    return hcTrans_;
-}
-
-
-template<class ParcelType>
-inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ThermoCloud<ParcelType>::Shc() const
-{
-    tmp<DimensionedField<scalar, volMesh> > tShc
-    (
-        new DimensionedField<scalar, volMesh>
-        (
-            IOobject
-            (
-                this->name() + "Shc",
-                this->db().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            this->mesh(),
-            dimensionedScalar
-            (
-                "zero",
-                dimMass/dimLength/pow3(dimTime),
-                0.0
-            )
-        )
-    );
-
-    scalarField& Shc = tShc().field();
-    Shc = hcTrans_/(this->mesh().V()*this->db().time().deltaT());
-
-    return tShc;
-}
-
-
 template<class ParcelType>
 inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
 Foam::ThermoCloud<ParcelType>::Sh() const
@@ -178,19 +102,10 @@ Foam::ThermoCloud<ParcelType>::Sh() const
                 IOobject::AUTO_WRITE,
                 false
             ),
-            this->mesh(),
-            dimensionedScalar
-            (
-                "zero",
-                dimMass/dimLength/pow3(dimTime),
-                0.0
-            )
+            hsTrans_/(this->mesh().V()*this->db().time().deltaT())
         )
     );
 
-    scalarField& Sh = tSh().field();
-    Sh = (hsTrans_ + hcTrans_)/(this->mesh().V()*this->db().time().deltaT());
-
     return tSh;
 }
 
@@ -205,7 +120,7 @@ Foam::ThermoCloud<ParcelType>::Ep() const
         (
             IOobject
             (
-                this->name() + "radiationEp",
+                this->name() + "radiation::Ep",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -248,7 +163,7 @@ Foam::ThermoCloud<ParcelType>::ap() const
         (
             IOobject
             (
-                this->name() + "radiationAp",
+                this->name() + "radiation::ap",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -291,7 +206,7 @@ Foam::ThermoCloud<ParcelType>::sigmap() const
         (
             IOobject
             (
-                this->name() + "radiationSigmap",
+                this->name() + "radiation::sigmap",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index c1914e5ed5..8d32c6a709 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -267,7 +267,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         YLiquid_,
         dMassPC,
         Sh,
-        dhsTrans,
         Ne,
         NCpW,
         Cs
@@ -296,7 +295,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         canCombust_,
         dMassDV,
         Sh,
-        dhsTrans,
         Ne,
         NCpW,
         Cs
@@ -398,19 +396,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(GAS, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassGas[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
         forAll(YLiquid_, i)
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(LIQ, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassLiquid[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
 /*
         // No mapping between solid components and carrier phase
@@ -418,19 +408,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(SLD, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassSolid[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
 */
         forAll(dMassSRCarrier, i)
         {
             td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassSRCarrier[i]
-               *td.cloud().mcCarrierThermo().speciesData()[i].Hc();
         }
 
         // Update momentum transfer
@@ -476,7 +458,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 */
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
             td.cloud().hsTrans()[cellI] +=
-                np0*mass1*HEff(td, pc, T1, idG, idL, idS);
+                np0*mass1*HEff(td, pc, T1, idG, idL, idS); // using total h
         }
     }
 
@@ -520,7 +502,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     bool& canCombust,
     scalarField& dMassDV,
     scalar& Sh,
-    scalar& dhsTrans,
     scalar& N,
     scalar& NCpW,
     scalarField& Cs
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index 1bc1ee99d7..3951688cc3 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -238,7 +238,6 @@ protected:
             bool& canCombust,          // 'can combust' flag
             scalarField& dMassDV,      // mass transfer - local to particle
             scalar& Sh,                // explicit particle enthalpy source
-            scalar& dhsTrans,          // sensible enthalpy transfer to carrier
             scalar& N,                 // flux of species emitted from particle
             scalar& NCpW,              // sum of N*Cp*W of emission species
             scalarField& Cs            // carrier conc. of emission species
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 92525ea2b1..edfb377176 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -287,7 +287,6 @@ void Foam::ReactingParcel<ParcelType>::calc
         Y_,
         dMassPC,
         Sh,
-        dhsTrans,
         Ne,
         NCpW,
         Cs
@@ -341,10 +340,6 @@ void Foam::ReactingParcel<ParcelType>::calc
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(0, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassPC[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
 
         // Update momentum transfer
@@ -371,7 +366,7 @@ void Foam::ReactingParcel<ParcelType>::calc
                 td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
             }
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
-            td.cloud().hcTrans()[cellI] +=
+            td.cloud().hsTrans()[cellI] +=
                 np0*mass1*td.cloud().composition().H(0, Y_, pc_, T1);
         }
     }
@@ -417,7 +412,6 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     const scalarField& YComponents,
     scalarField& dMassPC,
     scalar& Sh,
-    scalar& dhsTrans,               // TODO: not used
     scalar& N,
     scalar& NCpW,
     scalarField& Cs
@@ -469,6 +463,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
             td.cloud().composition().localToGlobalCarrierId(idPhase, i);
         const label idl = td.cloud().composition().globalIds(idPhase)[i];
 
+        // Calculate enthalpy transfer
         if
         (
             td.cloud().phaseChange().enthalpyTransfer()
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index 41d7c6d30c..1c7ed466d2 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -209,7 +209,6 @@ protected:
             const scalarField& YComponents, // component mass fractions
             scalarField& dMassPC,      // mass transfer - local to particle
             scalar& Sh,                // explicit particle enthalpy source
-            scalar& dhsTrans,          // sensible enthalpy transfer to carrier
             scalar& N,                 // flux of species emitted from particle
             scalar& NCpW,              // sum of N*Cp*W of emission species
             scalarField& Cs            // carrier conc. of emission species
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
index 6ae3dabb6f..e06b4925be 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
@@ -410,6 +410,130 @@ Foam::scalar Foam::CompositionModel<CloudType>::H
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::CompositionModel<CloudType>::Hs
+(
+    const label phaseI,
+    const scalarField& Y,
+    const scalar p,
+    const scalar T
+) const
+{
+    const phaseProperties& props = phaseProps_[phaseI];
+    scalar HsMixture = 0.0;
+    switch (props.phase())
+    {
+        case phaseProperties::GAS:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HsMixture += Y[i]*mcCarrierThermo_.speciesData()[gid].Hs(T);
+            }
+            break;
+        }
+        case phaseProperties::LIQUID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HsMixture +=
+                    Y[i]
+                   *(
+                       this->liquids().properties()[gid].h(p, T)
+                     - this->liquids().properties()[gid].h(p, 298.25)
+                    );
+            }
+            break;
+        }
+        case phaseProperties::SOLID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HsMixture += Y[i]*this->solids().properties()[gid].cp()*T;
+            }
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "Foam::scalar Foam::CompositionModel<CloudType>::Hs"
+                "("
+                "    const label, "
+                "    const scalarField&, "
+                "    const scalar, "
+                "    const scalar"
+                ") const"
+            )   << "Unknown phase enumeration" << nl << abort(FatalError);
+        }
+    }
+
+    return HsMixture;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::CompositionModel<CloudType>::Hc
+(
+    const label phaseI,
+    const scalarField& Y,
+    const scalar p,
+    const scalar T
+) const
+{
+    const phaseProperties& props = phaseProps_[phaseI];
+    scalar HcMixture = 0.0;
+    switch (props.phase())
+    {
+        case phaseProperties::GAS:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HcMixture += Y[i]*mcCarrierThermo_.speciesData()[gid].Hc();
+            }
+            break;
+        }
+        case phaseProperties::LIQUID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HcMixture +=
+                    Y[i]*this->liquids().properties()[gid].h(p, 298.15);
+            }
+            break;
+        }
+        case phaseProperties::SOLID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HcMixture += Y[i]*this->solids().properties()[gid].Hf();
+            }
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "Foam::scalar Foam::CompositionModel<CloudType>::Hc"
+                "("
+                "    const label, "
+                "    const scalarField&, "
+                "    const scalar, "
+                "    const scalar"
+                ") const"
+            )   << "Unknown phase enumeration" << nl << abort(FatalError);
+        }
+    }
+
+    return HcMixture;
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::CompositionModel<CloudType>::cp
 (
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
index b3f0949def..a3a8e9c724 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
@@ -226,7 +226,7 @@ public:
 
         // Evaluation
 
-            //- Return enthalpy for the phase phaseI
+            //- Return total enthalpy for the phase phaseI
             virtual scalar H
             (
                 const label phaseI,
@@ -235,6 +235,24 @@ public:
                 const scalar T
             ) const;
 
+            //- Return sensible enthalpy for the phase phaseI
+            virtual scalar Hs
+            (
+                const label phaseI,
+                const scalarField& Y,
+                const scalar p,
+                const scalar T
+            ) const;
+
+            //- Return chemical enthalpy for the phase phaseI
+            virtual scalar Hc
+            (
+                const label phaseI,
+                const scalarField& Y,
+                const scalar p,
+                const scalar T
+            ) const;
+
             //- Return specific heat caoacity for the phase phaseI
             virtual scalar cp
             (
diff --git a/src/parallel/reconstruct/pointFieldReconstructor.H b/src/parallel/reconstruct/pointFieldReconstructor.H
index 7ab04c3d41..3be8ba3a7a 100644
--- a/src/parallel/reconstruct/pointFieldReconstructor.H
+++ b/src/parallel/reconstruct/pointFieldReconstructor.H
@@ -138,7 +138,11 @@ public:
 
         //- Reconstruct and write all fields
         template<class Type>
-        void reconstructFields(const IOobjectList& objects);
+        void reconstructFields
+        (
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
+        );
 };
 
 
diff --git a/src/parallel/reconstruct/pointFieldReconstructorReconstructFields.C b/src/parallel/reconstruct/pointFieldReconstructorReconstructFields.C
index c96fcfa7ed..9b8e161483 100644
--- a/src/parallel/reconstruct/pointFieldReconstructorReconstructFields.C
+++ b/src/parallel/reconstruct/pointFieldReconstructorReconstructFields.C
@@ -143,7 +143,8 @@ Foam::pointFieldReconstructor::reconstructField(const IOobject& fieldIoObject)
 template<class Type>
 void Foam::pointFieldReconstructor::reconstructFields
 (
-    const IOobjectList& objects
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
 )
 {
     word fieldClassName
@@ -157,16 +158,18 @@ void Foam::pointFieldReconstructor::reconstructFields
     {
         Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
 
-        for
-        (
-            IOobjectList::iterator fieldIter = fields.begin();
-            fieldIter != fields.end();
-            ++fieldIter
-        )
+        forAllConstIter(IOobjectList, fields, fieldIter)
         {
-            Info<< "        " << fieldIter()->name() << endl;
+            if
+            (
+                !selectedFields.size()
+             || selectedFields.found(fieldIter()->name())
+            )
+            {
+                Info<< "        " << fieldIter()->name() << endl;
 
-            reconstructField<Type>(*fieldIter())().write();
+                reconstructField<Type>(*fieldIter())().write();
+            }
         }
 
         Info<< endl;
diff --git a/src/thermophysicalModels/basic/Make/files b/src/thermophysicalModels/basic/Make/files
index dcf9c8f4c0..59f8f822ea 100644
--- a/src/thermophysicalModels/basic/Make/files
+++ b/src/thermophysicalModels/basic/Make/files
@@ -6,11 +6,13 @@ basicThermo/basicThermo.C
 psiThermo/basicPsiThermo/basicPsiThermo.C
 psiThermo/basicPsiThermo/newBasicPsiThermo.C
 psiThermo/hPsiThermo/hPsiThermos.C
+psiThermo/hsPsiThermo/hsPsiThermos.C
 psiThermo/ePsiThermo/ePsiThermos.C
 
 rhoThermo/basicRhoThermo/basicRhoThermo.C
 rhoThermo/basicRhoThermo/newBasicRhoThermo.C
 rhoThermo/hRhoThermo/hRhoThermos.C
+rhoThermo/hsRhoThermo/hsRhoThermos.C
 
 derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C
 derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C
diff --git a/src/thermophysicalModels/basic/basicThermo/basicThermo.C b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
index f03bcc8289..67f0de7092 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.C
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
@@ -343,6 +343,13 @@ Foam::tmp<Foam::scalarField> Foam::basicThermo::hs
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::basicThermo::hc() const
+{
+    notImplemented("basicThermo::hc()");
+    return volScalarField::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 1c907edcf9..4f28779a30 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.H
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.H
@@ -145,17 +145,17 @@ public:
             //  Non-const access allowed for transport equations
             virtual volScalarField& h();
 
-            //- Enthalpy [J/kg]
+            //- Total enthalpy [J/kg]
             virtual const volScalarField& h() const;
 
-            //- Enthalpy for cell-set [J/kg]
+            //- Total enthalpy for cell-set [J/kg]
             virtual tmp<scalarField> h
             (
                 const scalarField& T,
                 const labelList& cells
             ) const;
 
-            //- Enthalpy for patch [J/kg]
+            //- Total enthalpy for patch [J/kg]
             virtual tmp<scalarField> h
             (
                 const scalarField& T,
@@ -166,23 +166,26 @@ public:
             //  Non-const access allowed for transport equations
             virtual volScalarField& hs();
 
-            //- Enthalpy [J/kg]
+            //- Sensible enthalpy [J/kg]
             virtual const volScalarField& hs() const;
 
-            //- Enthalpy for cell-set [J/kg]
+            //- Sensible enthalpy for cell-set [J/kg]
             virtual tmp<scalarField> hs
             (
                 const scalarField& T,
                 const labelList& cells
             ) const;
 
-            //- Enthalpy for patch [J/kg]
+            //- Sensible enthalpy for patch [J/kg]
             virtual tmp<scalarField> hs
             (
                 const scalarField& T,
                 const label patchi
             ) const;
 
+            //- Chemical enthalpy [J/kg]
+            virtual tmp<volScalarField> hc() const;
+
             //- Internal energy [J/kg]
             //  Non-const access allowed for transport equations
             virtual volScalarField& e();
@@ -197,7 +200,7 @@ public:
                 const labelList& cells
             ) const;
 
-            //-Internal energy for patch [J/kg]
+            //- Internal energy for patch [J/kg]
             virtual tmp<scalarField> e
             (
                 const scalarField& T,
diff --git a/src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermo.C b/src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermo.C
new file mode 100644
index 0000000000..2683e3f851
--- /dev/null
+++ b/src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermo.C
@@ -0,0 +1,347 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 "hsPsiThermo.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::hsPsiThermo<MixtureType>::calculate()
+{
+    const scalarField& hsCells = hs_.internalField();
+    const scalarField& pCells = this->p_.internalField();
+
+    scalarField& TCells = this->T_.internalField();
+    scalarField& psiCells = this->psi_.internalField();
+    scalarField& muCells = this->mu_.internalField();
+    scalarField& alphaCells = this->alpha_.internalField();
+
+    forAll(TCells, celli)
+    {
+        const typename MixtureType::thermoType& mixture_ =
+            this->cellMixture(celli);
+
+        TCells[celli] = mixture_.THs(hsCells[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 = this->p_.boundaryField()[patchi];
+        fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
+        fvPatchScalarField& ppsi = this->psi_.boundaryField()[patchi];
+
+        fvPatchScalarField& phs = hs_.boundaryField()[patchi];
+
+        fvPatchScalarField& pmu = this->mu_.boundaryField()[patchi];
+        fvPatchScalarField& palpha = this->alpha_.boundaryField()[patchi];
+
+        if (pT.fixesValue())
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+
+                phs[facei] = mixture_.Hs(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_.THs(phs[facei], pT[facei]);
+
+                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                palpha[facei] = mixture_.alpha(pT[facei]);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::hsPsiThermo<MixtureType>::hsPsiThermo(const fvMesh& mesh)
+:
+    basicPsiThermo(mesh),
+    MixtureType(*this, mesh),
+
+    hs_
+    (
+        IOobject
+        (
+            "hs",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimEnergy/dimMass,
+        this->hBoundaryTypes()
+    )
+{
+    scalarField& hsCells = hs_.internalField();
+    const scalarField& TCells = this->T_.internalField();
+
+    forAll(hsCells, celli)
+    {
+        hsCells[celli] = this->cellMixture(celli).Hs(TCells[celli]);
+    }
+
+    forAll(hs_.boundaryField(), patchi)
+    {
+        hs_.boundaryField()[patchi] ==
+            hs(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    hBoundaryCorrection(hs_);
+
+    calculate();
+
+    // Switch on saving old time
+    this->psi_.oldTime();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::hsPsiThermo<MixtureType>::~hsPsiThermo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::hsPsiThermo<MixtureType>::correct()
+{
+    if (debug)
+    {
+        Info<< "entering hsPsiThermo<MixtureType>::correct()" << endl;
+    }
+
+    // force the saving of the old-time values
+    this->psi_.oldTime();
+
+    calculate();
+
+    if (debug)
+    {
+        Info<< "exiting hsPsiThermo<MixtureType>::correct()" << endl;
+    }
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsPsiThermo<MixtureType>::hs
+(
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
+
+    forAll(T, celli)
+    {
+        hs[celli] = this->cellMixture(cells[celli]).Hs(T[celli]);
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsPsiThermo<MixtureType>::hs
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
+
+    forAll(T, facei)
+    {
+        hs[facei] = this->patchFaceMixture(patchi, facei).Hs(T[facei]);
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsPsiThermo<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::hsPsiThermo<MixtureType>::Cp() const
+{
+    const fvMesh& mesh = this->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),
+            this->T_.boundaryField().types()
+        )
+    );
+
+    volScalarField& cp = tCp();
+
+    forAll(this->T_, celli)
+    {
+        cp[celli] = this->cellMixture(celli).Cp(this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        const fvPatchScalarField& pT = this->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::hsPsiThermo<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::hsPsiThermo<MixtureType>::Cv() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCv
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cv",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimEnergy/dimMass/dimTemperature
+        )
+    );
+
+    volScalarField& cv = tCv();
+
+    forAll(this->T_, celli)
+    {
+        cv[celli] = this->cellMixture(celli).Cv(this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        cv.boundaryField()[patchi] =
+            Cv(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    return tCv;
+}
+
+
+template<class MixtureType>
+bool Foam::hsPsiThermo<MixtureType>::read()
+{
+    if (basicPsiThermo::read())
+    {
+        MixtureType::read(*this);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermo.H b/src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermo.H
new file mode 100644
index 0000000000..0c24ad84ce
--- /dev/null
+++ b/src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermo.H
@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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::hsPsiThermo
+
+Description
+    Sensible enthalpy for a mixture based on compressibility
+
+SourceFiles
+    hsPsiThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef hsPsiThermo_H
+#define hsPsiThermo_H
+
+#include "basicPsiThermo.H"
+#include "basicMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class hsPsiThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MixtureType>
+class hsPsiThermo
+:
+    public basicPsiThermo,
+    public MixtureType
+{
+    // Private data
+
+        //- Sensible enthalpy field [J/kg]
+        volScalarField hs_;
+
+
+    // Private member functions
+
+        //- Calculate the thermo variables
+        void calculate();
+
+        //- Construct as copy (not implemented)
+        hsPsiThermo(const hsPsiThermo<MixtureType>&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("hsPsiThermo");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        hsPsiThermo(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~hsPsiThermo();
+
+
+    // Member functions
+
+        //- Return the compostion of the mixture
+        virtual basicMixture& composition()
+        {
+            return *this;
+        }
+
+        //- Return the compostion of the mixture
+        virtual const basicMixture& composition() const
+        {
+            return *this;
+        }
+
+        //- Update properties
+        virtual void correct();
+
+
+        // Access to thermodynamic state variables
+
+            //- Sensible enthalpy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& hs()
+            {
+                return hs_;
+            }
+
+            //- Sensible enthalpy [J/kg]
+            virtual const volScalarField& hs() const
+            {
+                return hs_;
+            }
+
+
+        // Fields derived from thermodynamic state variables
+
+            //- 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;
+
+            //- 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 "hsPsiThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermos.C b/src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermos.C
similarity index 75%
rename from applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermos.C
rename to src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermos.C
index bf3bc4ebfe..3fad4dd806 100644
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermos.C
+++ b/src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermos.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,10 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "makeHsCombustionThermo.H"
-
-#include "hsCombustionThermo.H"
-#include "hsPsiMixtureThermo.H"
+#include "makeBasicPsiThermo.H"
 
 #include "perfectGas.H"
 
@@ -35,32 +32,47 @@ License
 #include "janafThermo.H"
 #include "specieThermo.H"
 
+#include "constTransport.H"
 #include "sutherlandTransport.H"
 
-#include "veryInhomogeneousMixture.H"
-
-#include "multiComponentMixture.H"
-
-#include "thermoPhysicsTypes.H"
+#include "hsPsiThermo.H"
+#include "pureMixture.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
 
-
-makeHsCombustionThermo
+makeBasicPsiThermo
 (
-    hsCombustionThermo,
-    hsPsiMixtureThermo,
-    veryInhomogeneousMixture,
+    hsPsiThermo,
+    pureMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicPsiThermo
+(
+    hsPsiThermo,
+    pureMixture,
+    sutherlandTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicPsiThermo
+(
+    hsPsiThermo,
+    pureMixture,
     sutherlandTransport,
     janafThermo,
     perfectGas
 );
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermo.C b/src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermo.C
new file mode 100644
index 0000000000..1e5e118077
--- /dev/null
+++ b/src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermo.C
@@ -0,0 +1,346 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 "hsRhoThermo.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::hsRhoThermo<MixtureType>::calculate()
+{
+    const scalarField& hsCells = this->hs_.internalField();
+    const scalarField& pCells = this->p_.internalField();
+
+    scalarField& TCells = this->T_.internalField();
+    scalarField& psiCells = this->psi_.internalField();
+    scalarField& rhoCells = this->rho_.internalField();
+    scalarField& muCells = this->mu_.internalField();
+    scalarField& alphaCells = this->alpha_.internalField();
+
+    forAll(TCells, celli)
+    {
+        const typename MixtureType::thermoType& mixture_ =
+            this->cellMixture(celli);
+
+        TCells[celli] = mixture_.THs(hsCells[celli], TCells[celli]);
+        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
+        rhoCells[celli] = mixture_.rho(pCells[celli], TCells[celli]);
+
+        muCells[celli] = mixture_.mu(TCells[celli]);
+        alphaCells[celli] = mixture_.alpha(TCells[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+        fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
+        fvPatchScalarField& ppsi = this->psi_.boundaryField()[patchi];
+        fvPatchScalarField& prho = this->rho_.boundaryField()[patchi];
+
+        fvPatchScalarField& phs = this->hs_.boundaryField()[patchi];
+
+        fvPatchScalarField& pmu = this->mu_.boundaryField()[patchi];
+        fvPatchScalarField& palpha = this->alpha_.boundaryField()[patchi];
+
+        if (pT.fixesValue())
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+
+                phs[facei] = mixture_.Hs(pT[facei]);
+
+                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
+                prho[facei] = mixture_.rho(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_.THs(phs[facei], pT[facei]);
+
+                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
+                prho[facei] = mixture_.rho(pp[facei], pT[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                palpha[facei] = mixture_.alpha(pT[facei]);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::hsRhoThermo<MixtureType>::hsRhoThermo(const fvMesh& mesh)
+:
+    basicRhoThermo(mesh),
+    MixtureType(*this, mesh),
+
+    hs_
+    (
+        IOobject
+        (
+            "hs",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimEnergy/dimMass,
+        this->hBoundaryTypes()
+    )
+{
+    scalarField& hsCells = hs_.internalField();
+    const scalarField& TCells = this->T_.internalField();
+
+    forAll(hsCells, celli)
+    {
+        hsCells[celli] = this->cellMixture(celli).Hs(TCells[celli]);
+    }
+
+    forAll(hs_.boundaryField(), patchi)
+    {
+        hs_.boundaryField()[patchi] ==
+            hs(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    hBoundaryCorrection(hs_);
+
+    calculate();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::hsRhoThermo<MixtureType>::~hsRhoThermo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::hsRhoThermo<MixtureType>::correct()
+{
+    if (debug)
+    {
+        Info<< "entering hsRhoThermo<MixtureType>::correct()" << endl;
+    }
+
+    calculate();
+
+    if (debug)
+    {
+        Info<< "exiting hsRhoThermo<MixtureType>::correct()" << endl;
+    }
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsRhoThermo<MixtureType>::hs
+(
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
+
+    forAll(T, celli)
+    {
+        hs[celli] = this->cellMixture(cells[celli]).Hs(T[celli]);
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsRhoThermo<MixtureType>::hs
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
+
+    forAll(T, facei)
+    {
+        hs[facei] = this->patchFaceMixture(patchi, facei).Hs(T[facei]);
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsRhoThermo<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::hsRhoThermo<MixtureType>::Cp() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCp
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cp",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimEnergy/dimMass/dimTemperature,
+            this->T_.boundaryField().types()
+        )
+    );
+
+    volScalarField& cp = tCp();
+
+    forAll(this->T_, celli)
+    {
+        cp[celli] = this->cellMixture(celli).Cp(this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        const fvPatchScalarField& pT = this->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::hsRhoThermo<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::hsRhoThermo<MixtureType>::Cv() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCv
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cv",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimEnergy/dimMass/dimTemperature
+        )
+    );
+
+    volScalarField& cv = tCv();
+
+    forAll(this->T_, celli)
+    {
+        cv[celli] = this->cellMixture(celli).Cv(this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        cv.boundaryField()[patchi] =
+            Cv(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    return tCv;
+}
+
+
+template<class MixtureType>
+bool Foam::hsRhoThermo<MixtureType>::read()
+{
+    if (basicRhoThermo::read())
+    {
+        MixtureType::read(*this);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermo.H b/src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermo.H
new file mode 100644
index 0000000000..8ce55f07a6
--- /dev/null
+++ b/src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermo.H
@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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::hsRhoThermo
+
+Description
+    Sensible enthalpy for a mixture based on density
+
+SourceFiles
+    hsRhoThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef hsRhoThermo_H
+#define hsRhoThermo_H
+
+#include "basicRhoThermo.H"
+#include "basicMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class hsRhoThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MixtureType>
+class hsRhoThermo
+:
+    public basicRhoThermo,
+    public MixtureType
+{
+    // Private data
+
+        //- Sensible enthalpy field [J/kg]
+        volScalarField hs_;
+
+
+    // Private member functions
+
+        //- Calculate the thermo variables
+        void calculate();
+
+        //- Construct as copy (not implemented)
+        hsRhoThermo(const hsRhoThermo<MixtureType>&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("hsRhoThermo");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        hsRhoThermo(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~hsRhoThermo();
+
+
+    // 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
+
+            //- Sensible enthalpy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& hs()
+            {
+                return hs_;
+            }
+
+            //- Sensible enthalpy [J/kg]
+            virtual const volScalarField& hs() const
+            {
+                return hs_;
+            }
+
+
+        // Fields derived from thermodynamic state variables
+
+            //- Sensible enthalpy for cell-set [J/kg]
+            virtual tmp<scalarField> hs
+            (
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+            //- Sensible enthalpy for patch [J/kg]
+            virtual tmp<scalarField> hs
+            (
+                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 "hsRhoThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermos.C b/src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermos.C
new file mode 100644
index 0000000000..ba1741c537
--- /dev/null
+++ b/src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermos.C
@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 "makeBasicRhoThermo.H"
+
+#include "perfectGas.H"
+
+#include "hConstThermo.H"
+#include "janafThermo.H"
+#include "specieThermo.H"
+
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+
+#include "hsRhoThermo.H"
+#include "pureMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+makeBasicRhoThermo
+(
+    hsRhoThermo,
+    pureMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicRhoThermo
+(
+    hsRhoThermo,
+    pureMixture,
+    sutherlandTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicRhoThermo
+(
+    hsRhoThermo,
+    pureMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C
index 1b9b15fc8e..fd9c206053 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C
@@ -466,7 +466,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
         this->thermo().rho()
     );
 
-    tmp<volScalarField> tsource
+    tmp<volScalarField> ttc
     (
         new volScalarField
         (
@@ -484,7 +484,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
         )
     );
 
-    scalarField& t = tsource();
+    scalarField& tc = ttc();
 
     label nReaction = reactions_.size();
 
@@ -517,17 +517,58 @@ Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
                 forAll(R.rhs(), s)
                 {
                     scalar sr = R.rhs()[s].stoichCoeff;
-                    t[celli] += sr*pf*cf;
+                    tc[celli] += sr*pf*cf;
                 }
             }
-            t[celli] = nReaction*cSum/t[celli];
+            tc[celli] = nReaction*cSum/tc[celli];
         }
     }
 
 
-    tsource().correctBoundaryConditions();
+    ttc().correctBoundaryConditions();
 
-    return tsource;
+    return ttc;
+}
+
+
+template<class CompType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::ODEChemistryModel<CompType, ThermoType>::Sh() const
+{
+    tmp<volScalarField> tSh
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Sh",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    if (this->chemistry_)
+    {
+        scalarField& Sh = tSh();
+
+        forAll(Y_, i)
+        {
+            forAll(Sh, cellI)
+            {
+                scalar hi = specieThermo_[i].Hc();
+                Sh[cellI] -= hi*RR_[i][cellI];
+            }
+        }
+    }
+
+    return tSh;
 }
 
 
@@ -545,37 +586,19 @@ Foam::ODEChemistryModel<CompType, ThermoType>::dQ() const
                 this->mesh_.time().timeName(),
                 this->mesh_,
                 IOobject::NO_READ,
-                IOobject::NO_WRITE
+                IOobject::NO_WRITE,
+                false
             ),
             this->mesh_,
-            dimensionedScalar
-            (
-                "zero",
-                dimensionSet(0, 2, -3 , 0, 0, 0, 0),
-                0.0
-            )
+            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0),
+            zeroGradientFvPatchScalarField::typeName
         )
     );
 
     if (this->chemistry_)
     {
-        scalarField& dQ = tdQ();
-
-        scalarField rhoEff(dQ.size(), 0.0);
-
-        forAll(Y_, i)
-        {
-            forAll(dQ, cellI)
-            {
-                scalar Ti = this->thermo().T()[cellI];
-                scalar pi = this->thermo().p()[cellI];
-                rhoEff[cellI] += Y_[i][cellI]*specieThermo_[i].rho(pi, Ti);
-                scalar hi = specieThermo_[i].H(Ti);
-                dQ[cellI] -= hi*RR_[i][cellI];
-            }
-        }
-
-        dQ /= rhoEff;
+        volScalarField& dQ = tdQ();
+        dQ.dimensionedInternalField() = this->mesh_.V()*Sh()();
     }
 
     return tdQ;
@@ -678,6 +701,9 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
 
     scalar deltaTMin = GREAT;
 
+    tmp<volScalarField> thc = this->thermo().hc();
+    const scalarField& hc = thc();
+
     forAll(rho, celli)
     {
         for (label i=0; i<nSpecie_; i++)
@@ -687,7 +713,7 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
 
         scalar rhoi = rho[celli];
         scalar Ti = this->thermo().T()[celli];
-        scalar hi = this->thermo().h()[celli];
+        scalar hi = this->thermo().hs()[celli] + hc[celli];
         scalar pi = this->thermo().p()[celli];
 
         scalarField c(nSpecie_);
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.H
index 6a5820ff19..59ab6baed5 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.H
@@ -39,9 +39,9 @@ SourceFiles
 #ifndef ODEChemistryModel_H
 #define ODEChemistryModel_H
 
-#include "hCombustionThermo.H"
 #include "Reaction.H"
 #include "ODE.H"
+#include "volFieldsFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,6 +49,8 @@ namespace Foam
 {
 
 // Forward declaration of classes
+class fvMesh;
+
 template<class CompType, class ThermoType>
 class chemistrySolver;
 
@@ -180,6 +182,9 @@ public:
             //- Return the chemical time scale
             virtual tmp<volScalarField> tc() const;
 
+            //- Return source for enthalpy equation [kg/m/s3]
+            virtual tmp<volScalarField> Sh() const;
+
             //- Return the heat release, i.e. enthalpy/sec [m2/s3]
             virtual tmp<volScalarField> dQ() const;
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModelI.H
index b5264a28d4..712b067e29 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModelI.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModelI.H
@@ -24,6 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "volFields.H"
 #include "zeroGradientFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -96,7 +97,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::RR
                 IOobject::NO_WRITE
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimensionSet(1, -3, -1, 0, 0), 0.0),
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0),
             zeroGradientFvPatchScalarField::typeName
         )
     );
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
index c8fcdd9863..290c3b62e1 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
@@ -120,7 +120,7 @@ public:
 
             // Fields
 
-                //- Return const access to chemical source terms
+                //- Return const access to chemical source terms [kg/m3/s]
                 virtual tmp<volScalarField> RR(const label i) const = 0;
 
 
@@ -133,7 +133,10 @@ public:
                 //- Return the chemical time scale
                 virtual tmp<volScalarField> tc() const = 0;
 
-                //- Return the heat release
+                //- Return source for enthalpy equation [kg/m/s3]
+                virtual tmp<volScalarField> Sh() const = 0;
+
+                //- Return the heat release, i.e. enthalpy/sec [m2/s3]
                 virtual tmp<volScalarField> dQ() const = 0;
 };
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.C
index 7381d12a1d..11978f2ddb 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.C
@@ -45,7 +45,7 @@ Foam::psiChemistryModel::psiChemistryModel
 )
 :
     basicChemistryModel(mesh),
-    thermo_(hCombustionThermo::NewType(mesh, thermoTypeName))
+    thermo_(hsCombustionThermo::NewType(mesh, thermoTypeName))
 {}
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.H
index 4cf7af58ec..617d627aab 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.H
@@ -41,7 +41,7 @@ SourceFiles
 #include "basicChemistryModel.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
-#include "hCombustionThermo.H"
+#include "hsCombustionThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,7 +73,7 @@ protected:
     // Protected data
 
         //- Thermo package
-        autoPtr<hCombustionThermo> thermo_;
+        autoPtr<hsCombustionThermo> thermo_;
 
 
 public:
@@ -114,10 +114,10 @@ public:
     // Member Functions
 
         //- Return access to the thermo package
-        inline hCombustionThermo& thermo();
+        inline hsCombustionThermo& thermo();
 
         //- Return const access to the thermo package
-        inline const hCombustionThermo& thermo() const;
+        inline const hsCombustionThermo& thermo() const;
 };
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModelI.H
index a6122b8b68..bea7731228 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModelI.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModelI.H
@@ -26,13 +26,13 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::hCombustionThermo& Foam::psiChemistryModel::thermo()
+inline Foam::hsCombustionThermo& Foam::psiChemistryModel::thermo()
 {
     return thermo_();
 }
 
 
-inline const Foam::hCombustionThermo& Foam::psiChemistryModel::thermo() const
+inline const Foam::hsCombustionThermo& Foam::psiChemistryModel::thermo() const
 {
     return thermo_();
 }
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.C
index 73a3f00e9c..9684adaf77 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.C
@@ -45,7 +45,7 @@ Foam::rhoChemistryModel::rhoChemistryModel
 )
 :
     basicChemistryModel(mesh),
-    thermo_(hReactionThermo::NewType(mesh, thermoTypeName))
+    thermo_(hsReactionThermo::NewType(mesh, thermoTypeName))
 {}
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.H
index 8991f250f9..4f22904b8f 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.H
@@ -41,7 +41,7 @@ SourceFiles
 #include "basicChemistryModel.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
-#include "hReactionThermo.H"
+#include "hsReactionThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,7 +73,7 @@ protected:
     // Protected data
 
         //- Thermo package
-        autoPtr<hReactionThermo> thermo_;
+        autoPtr<hsReactionThermo> thermo_;
 
 
 public:
@@ -114,10 +114,10 @@ public:
     // Member Functions
 
         //- Return access to the thermo package
-        inline hReactionThermo& thermo();
+        inline hsReactionThermo& thermo();
 
         //- Return const access to the thermo package
-        inline const hReactionThermo& thermo() const;
+        inline const hsReactionThermo& thermo() const;
 };
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModelI.H
index d96171001b..cff2455ebc 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModelI.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModelI.H
@@ -26,13 +26,13 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::hReactionThermo& Foam::rhoChemistryModel::thermo()
+inline Foam::hsReactionThermo& Foam::rhoChemistryModel::thermo()
 {
     return thermo_();
 }
 
 
-inline const Foam::hReactionThermo& Foam::rhoChemistryModel::thermo() const
+inline const Foam::hsReactionThermo& Foam::rhoChemistryModel::thermo() const
 {
     return thermo_();
 }
diff --git a/src/thermophysicalModels/reactionThermo/Make/files b/src/thermophysicalModels/reactionThermo/Make/files
index c2f9612ec5..471dbb114d 100644
--- a/src/thermophysicalModels/reactionThermo/Make/files
+++ b/src/thermophysicalModels/reactionThermo/Make/files
@@ -8,6 +8,10 @@ combustionThermo/hCombustionThermo/hCombustionThermo.C
 combustionThermo/hCombustionThermo/newhCombustionThermo.C
 combustionThermo/hCombustionThermo/hCombustionThermos.C
 
+combustionThermo/hsCombustionThermo/hsCombustionThermo.C
+combustionThermo/hsCombustionThermo/newhsCombustionThermo.C
+combustionThermo/hsCombustionThermo/hsCombustionThermos.C
+
 combustionThermo/hhuCombustionThermo/hhuCombustionThermo.C
 combustionThermo/hhuCombustionThermo/newhhuCombustionThermo.C
 combustionThermo/hhuCombustionThermo/hhuCombustionThermos.C
@@ -16,6 +20,9 @@ reactionThermo/hReactionThermo/hReactionThermo.C
 reactionThermo/hReactionThermo/newhReactionThermo.C
 reactionThermo/hReactionThermo/hReactionThermos.C
 
+reactionThermo/hsReactionThermo/hsReactionThermo.C
+reactionThermo/hsReactionThermo/newhsReactionThermo.C
+reactionThermo/hsReactionThermo/hsReactionThermos.C
 
 derivedFvPatchFields/fixedUnburntEnthalpy/fixedUnburntEnthalpyFvPatchScalarField.C
 derivedFvPatchFields/gradientUnburntEnthalpy/gradientUnburntEnthalpyFvPatchScalarField.C
diff --git a/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/hCombustionThermo.H b/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/hCombustionThermo.H
index d6b1ff7069..cfcf5c3834 100644
--- a/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/hCombustionThermo.H
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/hCombustionThermo.H
@@ -128,9 +128,6 @@ public:
             }
 
 
-        //- Chemical enthalpy [J/kg]
-        virtual tmp<volScalarField> hc() const = 0;
-
         //- Update properties
         virtual void correct() = 0;
 };
diff --git a/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/makeCombustionThermo.H b/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/makeCombustionThermo.H
index 4ee0c2a022..39079cfb31 100644
--- a/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/makeCombustionThermo.H
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/makeCombustionThermo.H
@@ -25,8 +25,6 @@ License
 InClass
     Foam::hCombustionThermo
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef makeCombustionThermo_H
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.C b/src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.C
similarity index 96%
rename from applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.C
rename to src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.C
index 2e50a2a64c..55d8984301 100644
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.C
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.C
@@ -39,7 +39,7 @@ namespace Foam
 
 Foam::hsCombustionThermo::hsCombustionThermo(const fvMesh& mesh)
 :
-    basicSensiblePsiThermo(mesh),
+    basicPsiThermo(mesh),
 
     hs_
     (
@@ -52,7 +52,7 @@ Foam::hsCombustionThermo::hsCombustionThermo(const fvMesh& mesh)
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionSet(0, 2, -2, 0, 0),
+        dimEnergy/dimMass,
         this->hBoundaryTypes()
     )
 {}
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.H b/src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.H
similarity index 94%
rename from applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.H
rename to src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.H
index 44e4dfe433..2d7be88a74 100644
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.H
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.H
@@ -26,7 +26,7 @@ Class
     Foam::hsCombustionThermo
 
 Description
-    Foam::hsCombustionThermo
+    Sensible enthalpy variant of combustionThermo
 
 SourceFiles
     hsCombustionThermo.C
@@ -36,7 +36,7 @@ SourceFiles
 #ifndef hsCombustionThermo_H
 #define hsCombustionThermo_H
 
-#include "basicSensiblePsiThermo.H"
+#include "basicPsiThermo.H"
 #include "basicMultiComponentMixture.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
@@ -52,14 +52,14 @@ namespace Foam
 
 class hsCombustionThermo
 :
-    public basicSensiblePsiThermo
+    public basicPsiThermo
 {
 
 protected:
 
     // Protected data
 
-        //- sensible enthalpy field
+        //- Sensible enthalpy field
         volScalarField hs_;
 
 
@@ -127,8 +127,6 @@ public:
                 return hs_;
             }
 
-        //- Chemical enthalpy [J/kg]
-        virtual tmp<volScalarField> hc() const = 0;
 
         //- Update properties
         virtual void correct() = 0;
diff --git a/src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermos.C b/src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermos.C
new file mode 100644
index 0000000000..e3332c7bb8
--- /dev/null
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermos.C
@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeHsCombustionThermo.H"
+
+#include "hsCombustionThermo.H"
+#include "hsPsiMixtureThermo.H"
+
+#include "perfectGas.H"
+
+#include "hConstThermo.H"
+#include "janafThermo.H"
+#include "specieThermo.H"
+
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+
+#include "dieselMixture.H"
+#include "homogeneousMixture.H"
+#include "inhomogeneousMixture.H"
+#include "veryInhomogeneousMixture.H"
+
+#include "reactingMixture.H"
+#include "multiComponentMixture.H"
+
+#include "thermoPhysicsTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    homogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    inhomogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    veryInhomogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    homogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    inhomogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    veryInhomogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    dieselMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+// Multi-component thermo
+
+makeHsCombustionMixtureThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    multiComponentMixture,
+    gasThermoPhysics
+);
+
+
+// Multi-component reaction thermo
+
+makeHsCombustionMixtureThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    reactingMixture,
+    gasThermoPhysics
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/makeHsCombustionThermo.H b/src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/makeHsCombustionThermo.H
similarity index 96%
rename from applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/makeHsCombustionThermo.H
rename to src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/makeHsCombustionThermo.H
index a1507ba35d..41962dfc54 100644
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/makeHsCombustionThermo.H
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/makeHsCombustionThermo.H
@@ -25,15 +25,13 @@ License
 InClass
     Foam::hsCombustionThermo
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef makeHsCombustionThermo_H
 #define makeHsCombustionThermo_H
 
 #include "addToRunTimeSelectionTable.H"
-#include "basicSensiblePsiThermo.H"
+#include "basicPsiThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,7 +51,7 @@ defineTemplateTypeNameAndDebugWithName                                        \
                                                                               \
 addToRunTimeSelectionTable                                                    \
 (                                                                             \
-    basicSensiblePsiThermo,                                                   \
+    basicPsiThermo,                                                           \
     MixtureThermo##Mixture##Transport##Thermo##EqnOfState,                    \
     fvMesh                                                                    \
 );                                                                            \
@@ -80,7 +78,7 @@ defineTemplateTypeNameAndDebugWithName                                        \
                                                                               \
 addToRunTimeSelectionTable                                                    \
 (                                                                             \
-    basicSensiblePsiThermo,                                                   \
+    basicPsiThermo,                                                           \
     MixtureThermo##Mixture##ThermoPhys,                                       \
     fvMesh                                                                    \
 );                                                                            \
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/newhsCombustionThermo.C b/src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/newhsCombustionThermo.C
similarity index 100%
rename from applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/newhsCombustionThermo.C
rename to src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/newhsCombustionThermo.C
diff --git a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.C b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.C
index d354b3a986..b3623a7b17 100644
--- a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.C
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.C
@@ -280,7 +280,7 @@ Foam::hPsiMixtureThermo<MixtureType>::Cp() const
                 IOobject::NO_WRITE
             ),
             mesh,
-            dimensionSet(0, 2, -2, -1, 0)
+            dimEnergy/dimMass/dimTemperature
         )
     );
 
diff --git a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.H b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.H
index a30a0c1465..bf932f8a9b 100644
--- a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.H
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.H
@@ -94,12 +94,12 @@ public:
         //- Update properties
         virtual void correct();
 
-        //- Chemical enthalpy [J/kg]
-        virtual tmp<volScalarField> hc() const;
-
 
         // Fields derived from thermodynamic state variables
 
+            //- Chemical enthalpy [J/kg]
+            virtual tmp<volScalarField> hc() const;
+
             //- Enthalpy for cell-set [J/kg]
             virtual tmp<scalarField> h
             (
diff --git a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
index 6d7245c044..2f6de08965 100644
--- a/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
@@ -296,7 +296,7 @@ Foam::hhuMixtureThermo<MixtureType>::Cp() const
                 IOobject::NO_WRITE
             ),
             mesh,
-            dimensionSet(0, 2, -2, -1, 0)
+            dimEnergy/dimMass/dimTemperature
         )
     );
 
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.C b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.C
similarity index 97%
rename from applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.C
rename to src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.C
index 0e6b14e1b0..83a40c8349 100644
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.C
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.C
@@ -205,15 +205,15 @@ Foam::hsPsiMixtureThermo<MixtureType>::hs
     const labelList& cells
 ) const
 {
-    tmp<scalarField> th(new scalarField(T.size()));
-    scalarField& hs = th();
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
 
     forAll(T, celli)
     {
         hs[celli] = this->cellMixture(cells[celli]).Hs(T[celli]);
     }
 
-    return th;
+    return ths;
 }
 
 
@@ -225,15 +225,15 @@ Foam::hsPsiMixtureThermo<MixtureType>::hs
     const label patchi
 ) const
 {
-    tmp<scalarField> th(new scalarField(T.size()));
-    scalarField& hs = th();
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
 
     forAll(T, facei)
     {
         hs[facei] = this->patchFaceMixture(patchi, facei).Hs(T[facei]);
     }
 
-    return th;
+    return ths;
 }
 
 
@@ -277,7 +277,7 @@ Foam::hsPsiMixtureThermo<MixtureType>::Cp() const
                 IOobject::NO_WRITE
             ),
             mesh,
-            dimensionSet(0, 2, -2, -1, 0)
+            dimEnergy/dimMass/dimTemperature
         )
     );
 
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.H b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.H
similarity index 96%
rename from applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.H
rename to src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.H
index e114b84734..b022b0eeb2 100644
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.H
+++ b/src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.H
@@ -44,7 +44,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class hsPsiMixtureThermo Declaration
+                   Class hsPsiMixtureThermo Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class MixtureType>
@@ -95,11 +95,11 @@ public:
         virtual void correct();
 
 
-        //- Chemical enthalpy [J/kg]
-        virtual tmp<volScalarField> hc() const;
-
         // Fields derived from thermodynamic state variables
 
+            //- Chemical enthalpy [J/kg]
+            virtual tmp<volScalarField> hc() const;
+
             //- Sensible enthalpy for cell-set [J/kg]
             virtual tmp<scalarField> hs
             (
diff --git a/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/hReactionThermo.H b/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/hReactionThermo.H
index f2319efa4f..ef801b2e96 100644
--- a/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/hReactionThermo.H
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/hReactionThermo.H
@@ -128,9 +128,6 @@ public:
             }
 
 
-        //- Chemical enthalpy [J/kg]
-        virtual tmp<volScalarField> hc() const = 0;
-
         //- Update properties
         virtual void correct() = 0;
 };
diff --git a/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/makeReactionThermo.H b/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/makeReactionThermo.H
index a74bded0bf..214602c844 100644
--- a/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/makeReactionThermo.H
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/makeReactionThermo.H
@@ -25,8 +25,6 @@ License
 InClass
     Foam::hReactionThermo
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef makeReactionThermo_H
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/basicSensiblePsiThermo.C b/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermo.C
similarity index 71%
rename from applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/basicSensiblePsiThermo.C
rename to src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermo.C
index 5d41d00c93..7f255de38a 100644
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/basicSensiblePsiThermo.C
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermo.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,28 +24,43 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "basicSensiblePsiThermo.H"
+#include "hsReactionThermo.H"
+#include "fvMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(basicSensiblePsiThermo, 0);
-    defineRunTimeSelectionTable(basicSensiblePsiThermo, fvMesh);
+    defineTypeNameAndDebug(hsReactionThermo, 0);
+    defineRunTimeSelectionTable(hsReactionThermo, fvMesh);
 }
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::basicSensiblePsiThermo::basicSensiblePsiThermo(const fvMesh& mesh)
+Foam::hsReactionThermo::hsReactionThermo(const fvMesh& mesh)
 :
-    basicThermo(mesh)
+    basicRhoThermo(mesh),
+
+    hs_
+    (
+        IOobject
+        (
+            "hs",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimEnergy/dimMass,
+        this->hBoundaryTypes()
+    )
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::basicSensiblePsiThermo::~basicSensiblePsiThermo()
+Foam::hsReactionThermo::~hsReactionThermo()
 {}
 
 
diff --git a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/basicSensiblePsiThermo.H b/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermo.H
similarity index 53%
rename from applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/basicSensiblePsiThermo.H
rename to src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermo.H
index f2eadec321..ae038b0b31 100644
--- a/applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/basicSensiblePsiThermo.H
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,21 +23,22 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::basicSensiblePsiThermo
+    Foam::hsReactionThermo
 
 Description
-    Basic thermodynamic properties based on compressibility
+    Sensible enthalpy variant for reactionThermo
 
 SourceFiles
-    basicSensiblePsiThermo.C
-    newBasicSensiblePsiThermo.C
+    hReactionThermo.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef basicSensiblePsiThermo_H
-#define basicSensiblePsiThermo_H
+#ifndef hsReactionThermo_H
+#define hsReactionThermo_H
 
-#include "basicThermo.H"
+#include "basicRhoThermo.H"
+#include "basicMultiComponentMixture.H"
+#include "autoPtr.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -46,60 +47,89 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class basicSensiblePsiThermo Declaration
+                     Class hsReactionThermo Declaration
 \*---------------------------------------------------------------------------*/
 
-class basicSensiblePsiThermo
+class hsReactionThermo
 :
-    public basicThermo
+    public basicRhoThermo
 {
 
 protected:
 
-    // Protected member functions
+    // Protected data
 
-        //- Construct as copy (not implemented)
-        basicSensiblePsiThermo(const basicSensiblePsiThermo&);
+        //- Sensible enthalpy [J/kg]
+        volScalarField hs_;
 
 
 public:
 
     //- Runtime type information
-    TypeName("basicSensiblePsiThermo");
+    TypeName("hsReactionThermo");
 
 
-    //- Declare run-time constructor selection table
+    //- Declare run-time constructor selection tables
     declareRunTimeSelectionTable
     (
         autoPtr,
-        basicSensiblePsiThermo,
+        hsReactionThermo,
         fvMesh,
         (const fvMesh& mesh),
         (mesh)
     );
 
+
     // Constructors
 
-        //- Construct from mesh
-        basicSensiblePsiThermo(const fvMesh&);
+        //- Construct from dictionary and mesh
+        hsReactionThermo(const fvMesh&);
 
 
-    //- Selector
-    static autoPtr<basicSensiblePsiThermo> New(const fvMesh&);
+    // Selectors
+
+        //- Standard selection based on fvMesh
+        static autoPtr<hsReactionThermo> New(const fvMesh&);
+
+        //- Select and check that package contains 'thermoType'
+        static autoPtr<hsReactionThermo> NewType
+        (
+            const fvMesh&,
+            const word& thermoType
+        );
+
 
     //- Destructor
-    virtual ~basicSensiblePsiThermo();
+    virtual ~hsReactionThermo();
 
 
     // Member functions
 
-        // Fields derived from thermodynamic state variables
+        //- Return the composition of the multi-component mixture
+        virtual basicMultiComponentMixture& composition() = 0;
 
-            //- Density [kg/m^3] - uses current value of pressure
-            virtual tmp<volScalarField> rho() const
+        //- Return the composition of the multi-component mixture
+        virtual const basicMultiComponentMixture& composition() const = 0;
+
+
+        // Access to thermodynamic state variables
+
+            //- Sensible enthalpy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& hs()
             {
-                return p_*psi();
+                return hs_;
             }
+
+            //- Sensible enthalpy [J/kg]
+            virtual const volScalarField& hs() const
+            {
+                return hs_;
+            }
+
+
+        //- Update properties
+        virtual void correct() = 0;
 };
 
 
diff --git a/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermos.C b/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermos.C
new file mode 100644
index 0000000000..be0e7967de
--- /dev/null
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermos.C
@@ -0,0 +1,171 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 "makeHsReactionThermo.H"
+
+#include "hsReactionThermo.H"
+#include "hsRhoMixtureThermo.H"
+
+#include "perfectGas.H"
+
+#include "hConstThermo.H"
+#include "janafThermo.H"
+#include "specieThermo.H"
+
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+
+#include "homogeneousMixture.H"
+#include "inhomogeneousMixture.H"
+#include "veryInhomogeneousMixture.H"
+#include "dieselMixture.H"
+#include "multiComponentMixture.H"
+#include "reactingMixture.H"
+
+#include "thermoPhysicsTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    homogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    inhomogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    veryInhomogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    homogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    inhomogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    veryInhomogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    dieselMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+
+// Multi-component thermo
+
+makeHsReactionMixtureThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    multiComponentMixture,
+    icoPoly8ThermoPhysics
+);
+
+makeHsReactionMixtureThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    multiComponentMixture,
+    gasThermoPhysics
+);
+
+
+// Multi-component reaction thermo
+
+makeHsReactionMixtureThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    reactingMixture,
+    icoPoly8ThermoPhysics
+);
+
+makeHsReactionMixtureThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    reactingMixture,
+    gasThermoPhysics
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/makeHsReactionThermo.H b/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/makeHsReactionThermo.H
new file mode 100644
index 0000000000..321d2cc862
--- /dev/null
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/makeHsReactionThermo.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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::hReactionThermo
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeReactionThermo_H
+#define makeReactionThermo_H
+
+#include "addToRunTimeSelectionTable.H"
+#include "basicRhoThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeHsReactionThermo(CThermo,MixtureThermo,Mixture,Transport,Thermo,EqnOfState)\
+                                                                              \
+typedef MixtureThermo                                                         \
+    <Mixture<Transport<specieThermo<Thermo<EqnOfState> > > > >                \
+    MixtureThermo##Mixture##Transport##Thermo##EqnOfState;                    \
+                                                                              \
+defineTemplateTypeNameAndDebugWithName                                        \
+(                                                                             \
+    MixtureThermo##Mixture##Transport##Thermo##EqnOfState,                    \
+    #MixtureThermo                                                            \
+        "<"#Mixture"<"#Transport"<specieThermo<"#Thermo"<"#EqnOfState">>>>>", \
+    0                                                                         \
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    basicRhoThermo,                                                           \
+    MixtureThermo##Mixture##Transport##Thermo##EqnOfState,                    \
+    fvMesh                                                                    \
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    CThermo,                                                                  \
+    MixtureThermo##Mixture##Transport##Thermo##EqnOfState,                    \
+    fvMesh                                                                    \
+)
+
+
+#define makeHsReactionMixtureThermo(CThermo,MixtureThermo,Mixture,ThermoPhys) \
+                                                                              \
+typedef MixtureThermo<Mixture<ThermoPhys> >                                   \
+    MixtureThermo##Mixture##ThermoPhys;                                       \
+                                                                              \
+defineTemplateTypeNameAndDebugWithName                                        \
+(                                                                             \
+    MixtureThermo##Mixture##ThermoPhys,                                       \
+    #MixtureThermo"<"#Mixture"<"#ThermoPhys">>",                              \
+    0                                                                         \
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    basicRhoThermo,                                                           \
+    MixtureThermo##Mixture##ThermoPhys,                                       \
+    fvMesh                                                                    \
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    CThermo,                                                                  \
+    MixtureThermo##Mixture##ThermoPhys,                                       \
+    fvMesh                                                                    \
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/newhsReactionThermo.C b/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/newhsReactionThermo.C
new file mode 100644
index 0000000000..00457f5674
--- /dev/null
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/newhsReactionThermo.C
@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 "hsReactionThermo.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::hsReactionThermo> Foam::hsReactionThermo::New
+(
+    const fvMesh& mesh
+)
+{
+    word hsReactionThermoTypeName;
+
+    // Enclose the creation of the thermophysicalProperties to ensure it is
+    // deleted before the turbulenceModel is created otherwise the dictionary
+    // is entered in the database twice
+    {
+        IOdictionary thermoDict
+        (
+            IOobject
+            (
+                "thermophysicalProperties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+        );
+
+        thermoDict.lookup("thermoType") >> hsReactionThermoTypeName;
+    }
+
+    Info<< "Selecting thermodynamics package " << hsReactionThermoTypeName
+        << endl;
+
+    fvMeshConstructorTable::iterator cstrIter =
+        fvMeshConstructorTablePtr_->find(hsReactionThermoTypeName);
+
+    if (cstrIter == fvMeshConstructorTablePtr_->end())
+    {
+        FatalErrorIn("hsReactionThermo::New(const fvMesh&)")
+            << "Unknown hsReactionThermo type "
+            << hsReactionThermoTypeName << nl << nl
+            << "Valid hsReactionThermo types are:" << nl
+            << fvMeshConstructorTablePtr_->sortedToc() << nl
+            << exit(FatalError);
+    }
+
+    return autoPtr<hsReactionThermo>(cstrIter()(mesh));
+}
+
+
+Foam::autoPtr<Foam::hsReactionThermo> Foam::hsReactionThermo::NewType
+(
+    const fvMesh& mesh,
+    const word& thermoType
+)
+{
+    word hsReactionThermoTypeName;
+
+    // Enclose the creation of the thermophysicalProperties to ensure it is
+    // deleted before the turbulenceModel is created otherwise the dictionary
+    // is entered in the database twice
+    {
+        IOdictionary thermoDict
+        (
+            IOobject
+            (
+                "thermophysicalProperties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+        );
+
+        thermoDict.lookup("thermoType") >> hsReactionThermoTypeName;
+
+        if (hsReactionThermoTypeName.find(thermoType) == string::npos)
+        {
+            wordList allModels = fvMeshConstructorTablePtr_->sortedToc();
+            DynamicList<word> validModels;
+            forAll(allModels, i)
+            {
+                if (allModels[i].find(thermoType) != string::npos)
+                {
+                    validModels.append(allModels[i]);
+                }
+            }
+
+            FatalErrorIn
+            (
+                "autoPtr<hsReactionThermo> hsReactionThermo::NewType"
+                "("
+                    "const fvMesh&, "
+                    "const word&"
+                ")"
+            )   << "Inconsistent thermo package selected:" << nl << nl
+                << hsReactionThermoTypeName << nl << nl << "Please select a "
+                << "thermo package based on " << thermoType
+                << ". Valid options include:" << nl << validModels << nl
+                << exit(FatalError);
+        }
+    }
+
+    Info<< "Selecting thermodynamics package " << hsReactionThermoTypeName
+        << endl;
+
+    fvMeshConstructorTable::iterator cstrIter =
+        fvMeshConstructorTablePtr_->find(hsReactionThermoTypeName);
+
+    if (cstrIter == fvMeshConstructorTablePtr_->end())
+    {
+        FatalErrorIn("hsReactionThermo::New(const fvMesh&)")
+            << "Unknown hsReactionThermo type "
+            << hsReactionThermoTypeName << nl << nl
+            << "Valid hsReactionThermo types are:" << nl
+            << fvMeshConstructorTablePtr_->sortedToc() << nl
+            << exit(FatalError);
+    }
+
+    return autoPtr<hsReactionThermo>(cstrIter()(mesh));
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hsRhoMixtureThermo/hsRhoMixtureThermo.C b/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hsRhoMixtureThermo/hsRhoMixtureThermo.C
new file mode 100644
index 0000000000..e17bd8d120
--- /dev/null
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hsRhoMixtureThermo/hsRhoMixtureThermo.C
@@ -0,0 +1,320 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 "hsRhoMixtureThermo.H"
+#include "fvMesh.H"
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::hsRhoMixtureThermo<MixtureType>::calculate()
+{
+    const scalarField& hsCells = hs_.internalField();
+    const scalarField& pCells = p_.internalField();
+
+    scalarField& TCells = T_.internalField();
+    scalarField& psiCells = psi_.internalField();
+    scalarField& rhoCells = rho_.internalField();
+    scalarField& muCells = mu_.internalField();
+    scalarField& alphaCells = alpha_.internalField();
+
+    forAll(TCells, celli)
+    {
+        const typename MixtureType::thermoType& mixture =
+            this->cellMixture(celli);
+
+        TCells[celli] = mixture.THs(hsCells[celli], TCells[celli]);
+        psiCells[celli] = mixture.psi(pCells[celli], TCells[celli]);
+        rhoCells[celli] = mixture.rho(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& prho = rho_.boundaryField()[patchi];
+
+        fvPatchScalarField& phs = hs_.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);
+
+                phs[facei] = mixture.Hs(pT[facei]);
+
+                ppsi[facei] = mixture.psi(pp[facei], pT[facei]);
+                prho[facei] = mixture.rho(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.THs(phs[facei], pT[facei]);
+
+                ppsi[facei] = mixture.psi(pp[facei], pT[facei]);
+                prho[facei] = mixture.rho(pp[facei], pT[facei]);
+                pmu_[facei] = mixture.mu(pT[facei]);
+                palpha_[facei] = mixture.alpha(pT[facei]);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::hsRhoMixtureThermo<MixtureType>::hsRhoMixtureThermo(const fvMesh& mesh)
+:
+    hsReactionThermo(mesh),
+    MixtureType(*this, mesh)
+{
+    scalarField& hsCells = hs_.internalField();
+    const scalarField& TCells = T_.internalField();
+
+    forAll(hsCells, celli)
+    {
+        hsCells[celli] = this->cellMixture(celli).Hs(TCells[celli]);
+    }
+
+    forAll(hs_.boundaryField(), patchi)
+    {
+        hs_.boundaryField()[patchi] == hs(T_.boundaryField()[patchi], patchi);
+    }
+
+    hBoundaryCorrection(hs_);
+
+    calculate();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::hsRhoMixtureThermo<MixtureType>::~hsRhoMixtureThermo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::hsRhoMixtureThermo<MixtureType>::correct()
+{
+    if (debug)
+    {
+        Info<< "entering hsRhoMixtureThermo<MixtureType>::correct()" << endl;
+    }
+
+    calculate();
+
+    if (debug)
+    {
+        Info<< "exiting hsRhoMixtureThermo<MixtureType>::correct()" << endl;
+    }
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField>
+Foam::hsRhoMixtureThermo<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,
+            hs_.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::hsRhoMixtureThermo<MixtureType>::hs
+(
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
+
+    forAll(T, celli)
+    {
+        hs[celli] = this->cellMixture(cells[celli]).Hs(T[celli]);
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField>
+Foam::hsRhoMixtureThermo<MixtureType>::hs
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
+
+    forAll(T, facei)
+    {
+        hs[facei] = this->patchFaceMixture(patchi, facei).Hs(T[facei]);
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField>
+Foam::hsRhoMixtureThermo<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::hsRhoMixtureThermo<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,
+            dimEnergy/dimMass/dimTemperature
+        )
+    );
+
+    volScalarField& cp = tCp();
+
+    scalarField& cpCells = cp.internalField();
+    const scalarField& TCells = T_.internalField();
+
+    forAll(TCells, celli)
+    {
+        cpCells[celli] = this->cellMixture(celli).Cp(TCells[celli]);
+    }
+
+    forAll(T_.boundaryField(), patchi)
+    {
+        cp.boundaryField()[patchi] = Cp(T_.boundaryField()[patchi], patchi);
+    }
+
+    return tCp;
+}
+
+
+template<class MixtureType>
+bool Foam::hsRhoMixtureThermo<MixtureType>::read()
+{
+    if (hsReactionThermo::read())
+    {
+        MixtureType::read(*this);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hsRhoMixtureThermo/hsRhoMixtureThermo.H b/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hsRhoMixtureThermo/hsRhoMixtureThermo.H
new file mode 100644
index 0000000000..f11b02e599
--- /dev/null
+++ b/src/thermophysicalModels/reactionThermo/reactionThermo/mixtureThermos/hsRhoMixtureThermo/hsRhoMixtureThermo.H
@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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::hsRhoMixtureThermo
+
+Description
+    Foam::hsRhoMixtureThermo
+
+SourceFiles
+    hsRhoMixtureThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef hsRhoMixtureThermo_H
+#define hsRhoMixtureThermo_H
+
+#include "hsReactionThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class hsRhoMixtureThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MixtureType>
+class hsRhoMixtureThermo
+:
+    public hsReactionThermo,
+    public MixtureType
+{
+    // Private member functions
+
+        void calculate();
+
+        //- Construct as copy (not implemented)
+        hsRhoMixtureThermo(const hsRhoMixtureThermo<MixtureType>&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("hsRhoMixtureThermo");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        hsRhoMixtureThermo(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~hsRhoMixtureThermo();
+
+
+    // Member functions
+
+        //- Return the compostion of the multi-component mixture
+        virtual basicMultiComponentMixture& composition()
+        {
+            return *this;
+        }
+
+        //- Return the compostion of the multi-component mixture
+        virtual const basicMultiComponentMixture& composition() const
+        {
+            return *this;
+        }
+
+        //- Update properties
+        virtual void correct();
+
+
+        // Fields derived from thermodynamic state variables
+
+            //- Chemical enthalpy [J/kg]
+            virtual tmp<volScalarField> hc() const;
+
+
+            //- Sensible nthalpy for cell-set [J/kg]
+            virtual tmp<scalarField> hs
+            (
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+            //- Sensible enthalpy for patch [J/kg]
+            virtual tmp<scalarField> hs
+            (
+                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;
+
+
+        //- Read thermophysicalProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "hsRhoMixtureThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/dieselFoam/aachenBomb/constant/thermophysicalProperties b/tutorials/combustion/dieselFoam/aachenBomb/constant/thermophysicalProperties
index 0cda0abcd2..e5fbb0bd78 100644
--- a/tutorials/combustion/dieselFoam/aachenBomb/constant/thermophysicalProperties
+++ b/tutorials/combustion/dieselFoam/aachenBomb/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hPsiMixtureThermo<reactingMixture<gasThermoPhysics>>;
+thermoType      hsPsiMixtureThermo<reactingMixture<gasThermoPhysics>>;
 
 CHEMKINFile     "$FOAM_CASE/chemkin/chem.inp";
 
diff --git a/tutorials/combustion/dieselFoam/aachenBomb/system/fvSchemes b/tutorials/combustion/dieselFoam/aachenBomb/system/fvSchemes
index fb78bfbd25..1ac217a857 100644
--- a/tutorials/combustion/dieselFoam/aachenBomb/system/fvSchemes
+++ b/tutorials/combustion/dieselFoam/aachenBomb/system/fvSchemes
@@ -39,7 +39,7 @@ divSchemes
     {
         fu limitedLinear 1;
         ft limitedLinear 1;
-        h limitedLinear 1;
+        hs limitedLinear 1;
     };
     div((muEff*dev2(grad(U).T()))) Gauss linear;
 }
@@ -50,7 +50,7 @@ laplacianSchemes
     laplacian(muEff,U) Gauss linear corrected;
     laplacian(muEff,ft) Gauss linear corrected;
     laplacian(muEff,fu) Gauss linear corrected;
-    laplacian(((alphah*mut)+alpha),h) Gauss linear corrected;
+    laplacian(((alphah*mut)+alpha),hs) Gauss linear corrected;
     laplacian((rho|A(U)),p) Gauss linear corrected;
     laplacian(rhoD,k) Gauss linear corrected;
     laplacian(rhoD,epsilon) Gauss linear corrected;
diff --git a/tutorials/combustion/dieselFoam/aachenBomb/system/fvSolution b/tutorials/combustion/dieselFoam/aachenBomb/system/fvSolution
index a1ebddd2ca..90bf7e8dac 100644
--- a/tutorials/combustion/dieselFoam/aachenBomb/system/fvSolution
+++ b/tutorials/combustion/dieselFoam/aachenBomb/system/fvSolution
@@ -31,7 +31,7 @@ solvers
         tolerance       1e-06;
     }
 
-    "(U|Yi|h|k|epsilon)"
+    "(U|Yi|hs|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/thermophysicalProperties b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/thermophysicalProperties
index 1586cb545c..6faf5b59ba 100644
--- a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/thermophysicalProperties
+++ b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType hPsiMixtureThermo<reactingMixture<gasThermoPhysics>>;
+thermoType hsPsiMixtureThermo<reactingMixture<gasThermoPhysics>>;
 
 inertSpecie N2;
 
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/thermophysicalProperties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/thermophysicalProperties
index 50b814b4dc..d1bada5e1d 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/thermophysicalProperties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hPsiMixtureThermo<reactingMixture<gasThermoPhysics>>;
+thermoType      hsPsiMixtureThermo<reactingMixture<gasThermoPhysics>>;
 
 chemistryReader foamChemistryReader;
 
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSchemes b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSchemes
index 1bd64bed7b..07351539c1 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSchemes
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSchemes
@@ -32,7 +32,7 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
     div(phiU,p)     Gauss linear;
-    div(phi,h)      Gauss upwind;
+    div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(U)          Gauss linear;
@@ -49,7 +49,7 @@ laplacianSchemes
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(DREff,R) Gauss linear corrected;
     laplacian((rho*(1|A(U))),p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
+    laplacian(alphaEff,hs) Gauss linear corrected;
 }
 
 interpolationSchemes
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSolution b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSolution
index d875b9f976..d98907822a 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSolution
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSolution
@@ -33,7 +33,7 @@ solvers
         relTol          0;
     }
 
-    "(U|h|R|k|epsilon)"
+    "(U|hs|R|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/Allclean b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/Allclean
index cb015f101f..1dac791f4b 100755
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/Allclean
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/Allclean
@@ -8,5 +8,8 @@ cleanCase
 # remove 0 directory
 rm -rf 0
 
+# remove post-processing folder
+rm -rf postProcessing
+
 # copy 0.org to 0
 cp -r 0.org 0
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/thermophysicalProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/thermophysicalProperties
index 8af627dc12..faa859ed18 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/thermophysicalProperties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hRhoMixtureThermo<reactingMixture<icoPoly8ThermoPhysics>>;
+thermoType      hsRhoMixtureThermo<reactingMixture<icoPoly8ThermoPhysics>>;
 
 chemistryReader foamChemistryReader;
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSchemes
index 1bd64bed7b..07351539c1 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSchemes
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSchemes
@@ -32,7 +32,7 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
     div(phiU,p)     Gauss linear;
-    div(phi,h)      Gauss upwind;
+    div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(U)          Gauss linear;
@@ -49,7 +49,7 @@ laplacianSchemes
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(DREff,R) Gauss linear corrected;
     laplacian((rho*(1|A(U))),p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
+    laplacian(alphaEff,hs) Gauss linear corrected;
 }
 
 interpolationSchemes
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
index 8083f38051..d1da25fa5a 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution
@@ -25,7 +25,7 @@ solvers
         relTol          0;
     }
 
-    "(U|h|k|epsilon)"
+    "(U|hs|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/energySourceProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/energySourceProperties
new file mode 100644
index 0000000000..e9e651eb49
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/energySourceProperties
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      energySourcesProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+(
+    // none
+);
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/pointMassSourcesProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/massSourceProperties
similarity index 90%
rename from tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/pointMassSourcesProperties
rename to tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/massSourceProperties
index 4a838a9e79..41a0251b15 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/pointMassSourcesProperties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/massSourceProperties
@@ -11,15 +11,13 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      pointMassSourcesProperties;
+    object      massSourcesProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-active          false;
-
-pointSources
 (
-//none
+    // none
 );
 
+
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/pointMassSourcesProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/momentumSourceProperties
similarity index 90%
rename from tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/pointMassSourcesProperties
rename to tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/momentumSourceProperties
index 2698204f7c..6decc46694 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/pointMassSourcesProperties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/momentumSourceProperties
@@ -11,15 +11,13 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      pointMassSourcesProperties;
+    object      momentumSourcesProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-active          false;
-
-pointSources
 (
-    //none
+    // none
 );
 
+
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/thermophysicalProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/thermophysicalProperties
index fbb4d94eb2..e83b1d1183 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/thermophysicalProperties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hRhoMixtureThermo<reactingMixture<icoPoly8ThermoPhysics>>;
+thermoType      hsRhoMixtureThermo<reactingMixture<icoPoly8ThermoPhysics>>;
 
 chemistryReader foamChemistryReader;
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes
index 99e07cad3d..e2d5c43f93 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes
@@ -24,8 +24,6 @@ gradSchemes
 {
     default         Gauss linear;
     grad(p)         Gauss linear;
-//    grad(U)         cellLimited Gauss linear 1;
-//    grad(k)         cellLimited Gauss linear 1;
 }
 
 divSchemes
@@ -34,7 +32,7 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
     div(phiU,p)     Gauss linear;
-    div(phi,h)      Gauss upwind;
+    div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,omega) Gauss upwind;
@@ -46,16 +44,6 @@ divSchemes
 laplacianSchemes
 {
     default         Gauss linear uncorrected;
-/*
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(mut,U) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DepsilonEff,omega) Gauss linear corrected;
-    laplacian(DREff,R) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
-*/
 }
 
 interpolationSchemes
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
index 6d53c94485..496ce336c9 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
@@ -49,7 +49,7 @@ solvers
         tolerance       1e-06;
         relTol          0;
     };
-    "(Yi|h|)"
+    "(Yi|hs)"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allclean b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allclean
index 4843b685f3..fbf6469b87 100755
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allclean
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allclean
@@ -2,8 +2,9 @@
 # Source tutorial clean functions
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
-# remove old time folders
-rm -rf 0 *[1-9]* processor*
+# remove old time and post-processing folders
+rm -rf 0 *[1-9]* processor* postProcessing
+
 
 # copy 0.org to 0
 cp -r 0.org 0
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/energySourceProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/energySourceProperties
new file mode 100644
index 0000000000..e9e651eb49
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/energySourceProperties
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      energySourcesProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+(
+    // none
+);
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/massSourceProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/massSourceProperties
new file mode 100644
index 0000000000..41a0251b15
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/massSourceProperties
@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      massSourcesProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+(
+    // none
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/momentumSourceProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/momentumSourceProperties
new file mode 100644
index 0000000000..6decc46694
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/momentumSourceProperties
@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      momentumSourcesProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+(
+    // none
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/thermophysicalProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/thermophysicalProperties
index fbb4d94eb2..e83b1d1183 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/thermophysicalProperties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hRhoMixtureThermo<reactingMixture<icoPoly8ThermoPhysics>>;
+thermoType      hsRhoMixtureThermo<reactingMixture<icoPoly8ThermoPhysics>>;
 
 chemistryReader foamChemistryReader;
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes
index 12543c5fd7..fb7c45a4c8 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
     div(phiU,p)     Gauss linear;
-    div(phi,h)      Gauss upwind;
+    div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,omega) Gauss upwind;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
index 1ccb5e9656..473d8c68f8 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
@@ -49,7 +49,7 @@ solvers
         tolerance       1e-06;
         relTol          0;
     };
-    "(Yi|h|)"
+    "(Yi|hs)"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes
index 1bd64bed7b..07351539c1 100644
--- a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes
+++ b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSchemes
@@ -32,7 +32,7 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
     div(phiU,p)     Gauss linear;
-    div(phi,h)      Gauss upwind;
+    div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(U)          Gauss linear;
@@ -49,7 +49,7 @@ laplacianSchemes
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(DREff,R) Gauss linear corrected;
     laplacian((rho*(1|A(U))),p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
+    laplacian(alphaEff,hs) Gauss linear corrected;
 }
 
 interpolationSchemes
diff --git a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSolution b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSolution
index 736a2e3e70..613b764228 100644
--- a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSolution
+++ b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/system/fvSolution
@@ -25,7 +25,7 @@ solvers
         relTol          0;
     }
 
-    "(U|h|k|epsilon)"
+    "(U|hs|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createFields.H b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createFields.H
index cebd44cfd4..46bbe2d917 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createFields.H
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/createFields.H
@@ -7,7 +7,7 @@
     basicPsiThermo& thermo = pThermo();
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& psi = thermo.psi();
 
     volScalarField rho
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/hEqn.H b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/hEqn.H
deleted file mode 100644
index 2b3b60ce98..0000000000
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/hEqn.H
+++ /dev/null
@@ -1,17 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + fvm::div(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        DpDt
-      + thermoCloud1.Sh()
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-}
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/hsEqn.H b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/hsEqn.H
new file mode 100644
index 0000000000..3a1c6e5755
--- /dev/null
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/hsEqn.H
@@ -0,0 +1,17 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + fvm::div(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+      + thermoCloud1.Sh()
+    );
+
+    hsEqn.relax();
+
+    hsEqn.solve();
+
+    thermo.correct();
+}
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam.C b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam.C
index b20e73b487..f93437068f 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam.C
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam/rhoPisoTwinParcelFoam.C
@@ -86,7 +86,7 @@ int main(int argc, char *argv[])
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-                #include "hEqn.H"
+                #include "hsEqn.H"
                 #include "pEqn.H"
             }
         }
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
index 1756757d0d..9460e2e028 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
@@ -99,6 +99,7 @@ StandardWallInteractionCoeffs
 
 RanzMarshallCoeffs
 {
+    // thermal shielding
     BirdCorrection  false;
 }
 
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermophysicalProperties b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermophysicalProperties
index 36a5ad071a..a56b7fee60 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermophysicalProperties
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermophysicalProperties
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hPsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>;
+thermoType      hsPsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>;
 
 mixture         air 1 28.9 1007 0 1.84e-05 0.7;
 
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/fvSchemes b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/fvSchemes
index 1bd64bed7b..840a87f8b5 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/fvSchemes
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/fvSchemes
@@ -32,7 +32,7 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
     div(phiU,p)     Gauss linear;
-    div(phi,h)      Gauss upwind;
+    div(phi,hs)     Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(U)          Gauss linear;
@@ -49,7 +49,7 @@ laplacianSchemes
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(DREff,R) Gauss linear corrected;
     laplacian((rho*(1|A(U))),p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
+    laplacian(alphaEff,hs) Gauss linear corrected;
 }
 
 interpolationSchemes
@@ -65,7 +65,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p;
 }
 
 
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/fvSolution b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/fvSolution
index 0cdd3e30cf..481d14fbf6 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/fvSolution
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/fvSolution
@@ -33,7 +33,7 @@ solvers
         relTol          0;
     }
 
-    "(U|h|R|k|epsilon)"
+    "(U|hs|R|k|epsilon)"
     {
         $p;
         tolerance       1e-05;