diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/EEqns.H
index 588bd7ddf4..ff72d8d94e 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/EEqns.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/EEqns.H
@@ -47,3 +47,9 @@
 }
 
 fluid.correctThermo();
+
+Info<< "  phase1.thermo().T(): " << min(phase1.thermo().T()).value()
+    << " - " << max(phase1.thermo().T()).value() << endl;
+
+Info<< "  phase2.thermo().T(): " << min(phase2.thermo().T()).value()
+    << " - " << max(phase2.thermo().T()).value() << endl;
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C
index 7e1693dfed..9b65bb82f7 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C
@@ -60,13 +60,13 @@ Foam::heatTransferModels::RanzMarshall::~RanzMarshall()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::heatTransferModels::RanzMarshall::K() const
+Foam::heatTransferModels::RanzMarshall::K(const scalar residualAlpha) const
 {
     volScalarField Nu(scalar(2) + 0.6*sqrt(pair_.Re())*cbrt(pair_.Pr()));
 
     return
         6.0
-       *max(pair_.dispersed(), residualAlpha_)
+       *max(pair_.dispersed(), residualAlpha)
        *pair_.continuous().kappa()
        *Nu
        /sqr(pair_.dispersed().d());
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.H
index c02126528a..a4987ab0af 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.H
@@ -79,7 +79,7 @@ public:
     // Member Functions
 
         //- The heat transfer function K used in the enthalpy equation
-        tmp<volScalarField> K() const;
+        tmp<volScalarField> K(const scalar residualAlpha) const;
 };
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/heatTransferModel/heatTransferModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/heatTransferModel/heatTransferModel.C
index 1621b6be19..3a75fbf3d1 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/heatTransferModel/heatTransferModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/heatTransferModel/heatTransferModel.C
@@ -65,4 +65,13 @@ Foam::heatTransferModel::~heatTransferModel()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::heatTransferModel::K() const
+{
+    return K(residualAlpha_.value());
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/heatTransferModel/heatTransferModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/heatTransferModel/heatTransferModel.H
index 6276e8aee1..130b82198b 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/heatTransferModel/heatTransferModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/heatTransferModel/heatTransferModel.H
@@ -117,7 +117,13 @@ public:
         //- The heat transfer function K used in the enthalpy equation
         //    ddt(alpha1*rho1*ha) + ... = ... K*(Ta - Tb)
         //    ddt(alpha2*rho2*hb) + ... = ... K*(Tb - Ta)
-        virtual tmp<volScalarField> K() const = 0;
+        tmp<volScalarField> K() const;
+
+        //- The heat transfer function K used in the enthalpy equation
+        //    ddt(alpha1*rho1*ha) + ... = ... K*(Ta - Tb)
+        //    ddt(alpha2*rho2*hb) + ... = ... K*(Tb - Ta)
+        //  with a specified residual volume fraction
+        virtual tmp<volScalarField> K(const scalar residualAlpha) const = 0;
 };
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/sphericalHeatTransfer/sphericalHeatTransfer.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/sphericalHeatTransfer/sphericalHeatTransfer.C
index 4eec6de028..0672f16ef8 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/sphericalHeatTransfer/sphericalHeatTransfer.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/sphericalHeatTransfer/sphericalHeatTransfer.C
@@ -65,11 +65,14 @@ Foam::heatTransferModels::sphericalHeatTransfer::~sphericalHeatTransfer()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::heatTransferModels::sphericalHeatTransfer::K() const
+Foam::heatTransferModels::sphericalHeatTransfer::K
+(
+    const scalar residualAlpha
+) const
 {
     return
         60.0
-       *max(pair_.dispersed(), residualAlpha_)
+       *max(pair_.dispersed(), residualAlpha)
        *pair_.continuous().kappa()
        /sqr(pair_.dispersed().d());
 }
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/sphericalHeatTransfer/sphericalHeatTransfer.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/sphericalHeatTransfer/sphericalHeatTransfer.H
index 88d098668f..f3164df9bd 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/sphericalHeatTransfer/sphericalHeatTransfer.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/sphericalHeatTransfer/sphericalHeatTransfer.H
@@ -79,7 +79,7 @@ public:
     // Member Functions
 
         //- The heat transfer function K used in the enthalpy equation
-        tmp<volScalarField> K() const;
+        tmp<volScalarField> K(const scalar residualAlpha) const;
 };
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
index da6c62ade5..36e31630a4 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
@@ -223,6 +223,66 @@ Foam::BlendedInterfacialModel<ModelType>::K() const
 }
 
 
+template<class ModelType>
+Foam::tmp<Foam::volScalarField>
+Foam::BlendedInterfacialModel<ModelType>::K(const scalar residualAlpha) const
+{
+    tmp<volScalarField> f1, f2;
+
+    if (model_.valid() || model1In2_.valid())
+    {
+        f1 = blending_.f1(phase1_, phase2_);
+    }
+
+    if (model_.valid() || model2In1_.valid())
+    {
+        f2 = blending_.f2(phase1_, phase2_);
+    }
+
+    tmp<volScalarField> x
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                ModelType::typeName + ":K",
+                phase1_.mesh().time().timeName(),
+                phase1_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            phase1_.mesh(),
+            dimensionedScalar("zero", ModelType::dimK, 0)
+        )
+    );
+
+    if (model_.valid())
+    {
+        x() += model_->K(residualAlpha)*(scalar(1) - f1() - f2());
+    }
+    if (model1In2_.valid())
+    {
+        x() += model1In2_->K(residualAlpha)*f1;
+    }
+    if (model2In1_.valid())
+    {
+        x() += model2In1_->K(residualAlpha)*f2;
+    }
+
+    if
+    (
+        correctFixedFluxBCs_
+     && (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    )
+    {
+        correctFixedFluxBCs(x());
+    }
+
+    return x;
+}
+
+
 template<class ModelType>
 Foam::tmp<Foam::surfaceScalarField>
 Foam::BlendedInterfacialModel<ModelType>::Kf() const
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.H
index bd2ec60235..4be8f0b347 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.H
@@ -136,6 +136,9 @@ public:
         //- Return the blended force coefficient
         tmp<volScalarField> K() const;
 
+        //- Return the blended force coefficient with a specified residual alpha
+        tmp<volScalarField> K(const scalar residualAlpha) const;
+
         //- Return the face blended force coefficient
         tmp<surfaceScalarField> Kf() const;
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/Make/files b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/Make/files
index 538bd7d637..9f8f4f1fdb 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/Make/files
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/Make/files
@@ -29,4 +29,6 @@ twoPhaseSystem/twoPhaseSystem.C
 twoPhaseSystem/newTwoPhaseSystem.C
 twoPhaseSystem/twoPhaseSystems.C
 
+reactionThermo/hRefConstThermos.C
+
 LIB = $(FOAM_LIBBIN)/libreactingTwoPhaseSystem
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
index ff4fb10e9d..8b2bbdb99e 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
@@ -192,11 +192,11 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
         const volVectorField& U2(pair.phase2().U());
 
         const volScalarField dmdt(this->dmdt(pair));
-        const volScalarField dmdt12(posPart(dmdt));
-        const volScalarField dmdt21(negPart(dmdt));
+        const volScalarField dmdt21(posPart(dmdt));
+        const volScalarField dmdt12(negPart(dmdt));
 
-        *eqns[pair.phase1().name()] += fvm::Sp(dmdt21, U1) - dmdt21*U2;
-        *eqns[pair.phase2().name()] += dmdt12*U1 - fvm::Sp(dmdt12, U2);
+        *eqns[pair.phase1().name()] += dmdt21*U2 - fvm::Sp(dmdt21, U1);
+        *eqns[pair.phase2().name()] -= dmdt12*U1 - fvm::Sp(dmdt12, U2);
     }
 
     return eqnsPtr;
@@ -283,7 +283,7 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
         }
     }
 
-    // Source term due to mass trasfer
+    // Source term due to mass transfer
     forAllConstIter
     (
         phaseSystem::phasePairTable,
@@ -308,17 +308,19 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
         const volScalarField& K2(phase2.K());
 
         const volScalarField dmdt(this->dmdt(pair));
-        const volScalarField dmdt12(posPart(dmdt));
-        const volScalarField dmdt21(negPart(dmdt));
+        const volScalarField dmdt21(posPart(dmdt));
+        const volScalarField dmdt12(negPart(dmdt));
         const volScalarField& Tf(*Tf_[pair]);
 
         *eqns[phase1.name()] +=
-            fvm::Sp(dmdt21, he1) + dmdt21*K1
-          - dmdt21*(phase2.thermo().he(phase2.thermo().p(), Tf) + K2);
+            dmdt21*(phase1.thermo().he(phase1.thermo().p(), Tf))
+          - fvm::Sp(dmdt21, he1)
+          + dmdt21*(K2 - K1);
 
-        *eqns[phase2.name()] +=
-            dmdt12*(phase1.thermo().he(phase1.thermo().p(), Tf) + K1)
-          - fvm::Sp(dmdt12, he2) - dmdt12*K2;
+        *eqns[phase2.name()] -=
+            dmdt12*(phase2.thermo().he(phase2.thermo().p(), Tf))
+          - fvm::Sp(dmdt12, he2)
+          + dmdt12*(K1 - K2);
     }
 
     return eqnsPtr;
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index bbce1962a1..ffdb41326e 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ThermalPhaseChangePhaseSystem.H"
+#include "fvCFD.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -36,7 +37,8 @@ ThermalPhaseChangePhaseSystem
 :
     HeatAndMassTransferPhaseSystem<BasePhaseSystem>(mesh),
     volatile_(this->lookup("volatile")),
-    saturationModel_(saturationModel::New(this->subDict("saturationModel")))
+    saturationModel_(saturationModel::New(this->subDict("saturationModel"))),
+    massTransfer_(this->lookup("massTransfer"))
 {}
 
 
@@ -143,37 +145,79 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
         const phaseModel& phase1 = pair.phase1();
         const phaseModel& phase2 = pair.phase2();
 
-        volScalarField& Tf = *this->Tf_[pair];
-        Tf = saturationModel_->Tsat(phase1.thermo().p());
-
-        Info<< "Tf." << pair.name()
-            << ": min = " << min(Tf.internalField())
-            << ", mean = " << average(Tf.internalField())
-            << ", max = " << max(Tf.internalField())
-            << endl;
-
-        volScalarField& dmdt(*this->dmdt_[pair]);
-
-        volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
-        volScalarField H2(this->heatTransferModels_[pair][pair.second()]->K());
-
         const volScalarField& T1(phase1.thermo().T());
         const volScalarField& T2(phase2.thermo().T());
 
         const volScalarField& he1(phase1.thermo().he());
         const volScalarField& he2(phase2.thermo().he());
 
-        volScalarField hef2(phase2.thermo().he(phase2.thermo().p(), Tf));
-        volScalarField hef1(phase1.thermo().he(phase1.thermo().p(), Tf));
+        volScalarField& dmdt(*this->dmdt_[pair]);
 
-        dmdt =
-            (H2*(Tf - T1) + H1*(Tf - T2))
-           /min
+        volScalarField& Tf = *this->Tf_[pair];
+
+        volScalarField hef1(phase1.thermo().he(phase1.thermo().p(), Tf));
+        volScalarField hef2(phase2.thermo().he(phase2.thermo().p(), Tf));
+
+        if (massTransfer_ )
+        {
+            volScalarField H1
+            (
+                this->heatTransferModels_[pair][pair.first()]->K(0)
+            );
+
+            volScalarField H2
+            (
+                this->heatTransferModels_[pair][pair.second()]->K(0)
+            );
+
+            Tf = saturationModel_->Tsat(phase1.thermo().p());
+
+            dmdt =
+                (H1*(Tf - T1) + H2*(Tf - T2))
+               /min
+                (
+                    (pos(dmdt)*he2 + neg(dmdt)*hef2)
+                  - (neg(dmdt)*he1 + pos(dmdt)*hef1),
+                    0.3*mag(hef2 - hef1)
+                );
+
+            Info<< "dmdt." << pair.name()
+                << ": min = " << min(dmdt.internalField())
+                << ", mean = " << average(dmdt.internalField())
+                << ", max = " << max(dmdt.internalField())
+                << ", integral = " << fvc::domainIntegrate(dmdt).value()
+                << endl;
+        }
+        else
+        {
+            dmdt == dimensionedScalar("0", dmdt.dimensions(), 0);
+        }
+
+        volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
+        volScalarField H2(this->heatTransferModels_[pair][pair.second()]->K());
+
+        // Limit the H[12] boundary field to avoid /0
+        const scalar HLimit = 1e-4;
+        H1.boundaryField() =
+            max(H1.boundaryField(), phase1.boundaryField()*HLimit);
+        H2.boundaryField() =
+            max(H2.boundaryField(), phase2.boundaryField()*HLimit);
+
+        volScalarField mDotL
+        (
+            dmdt*
             (
                 (pos(dmdt)*he2 + neg(dmdt)*hef2)
-              - (neg(dmdt)*he1 + pos(dmdt)*hef1),
-                0.3*(hef1 - hef2)
-            );
+              - (neg(dmdt)*he1 + pos(dmdt)*hef1)
+            )
+        );
+        Tf = (H1*T1 + H2*T2 + mDotL)/(H1 + H2);
+
+        Info<< "Tf." << pair.name()
+            << ": min = " << min(Tf.internalField())
+            << ", mean = " << average(Tf.internalField())
+            << ", max = " << max(Tf.internalField())
+            << endl;
     }
 }
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
index 9a0a0d109a..edae5a6ff6 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
@@ -43,6 +43,7 @@ SourceFiles
 
 #include "HeatAndMassTransferPhaseSystem.H"
 #include "saturationModel.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,6 +70,9 @@ protected:
         //- The saturation model used to evaluate Tsat = Tf
         autoPtr<saturationModel> saturationModel_;
 
+        // Mass transfer enabled
+        Switch massTransfer_;
+
 
 public:
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
new file mode 100644
index 0000000000..25bc49fd48
--- /dev/null
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeReactionThermo.H"
+#include "makeThermo.H"
+
+#include "rhoReactionThermo.H"
+#include "heRhoThermo.H"
+
+#include "specie.H"
+#include "perfectGas.H"
+#include "perfectFluid.H"
+#include "rhoConst.H"
+#include "incompressiblePerfectGas.H"
+#include "hConstThermo.H"
+#include "janafThermo.H"
+#include "sensibleEnthalpy.H"
+#include "absoluteEnthalpy.H"
+#include "absoluteInternalEnergy.H"
+#include "thermo.H"
+
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+
+#include "pureMixture.H"
+#include "homogeneousMixture.H"
+#include "inhomogeneousMixture.H"
+#include "veryInhomogeneousMixture.H"
+#include "multiComponentMixture.H"
+#include "reactingMixture.H"
+#include "singleStepReactingMixture.H"
+
+#include "thermoPhysicsTypes.H"
+
+#include "hRefConstThermo.H"
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Thermo type typedefs:
+
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            perfectGas<specie>
+        >,
+        sensibleEnthalpy
+    >
+> constRefGasHThermoPhysics;
+
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            perfectFluid<specie>
+        >,
+        sensibleEnthalpy
+    >
+> constRefFluidHThermoPhysics;
+
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            perfectGas<specie>
+        >,
+        sensibleInternalEnergy
+    >
+> constRefGasEThermoPhysics;
+
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            perfectFluid<specie>
+        >,
+        sensibleInternalEnergy
+    >
+> constRefFluidEThermoPhysics;
+
+
+// pureMixture, sensibleEnthalpy:
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleEnthalpy,
+    hRefConstThermo,
+    perfectGas,
+    specie
+);
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleEnthalpy,
+    hRefConstThermo,
+    perfectFluid,
+    specie
+);
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleEnthalpy,
+    hRefConstThermo,
+    rhoConst,
+    specie
+);
+
+
+// pureMixture, sensibleInternalEnergy:
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    hRefConstThermo,
+    perfectGas,
+    specie
+);
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    hRefConstThermo,
+    perfectFluid,
+    specie
+);
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    hRefConstThermo,
+    rhoConst,
+    specie
+);
+
+
+// multiComponentMixture, sensibleInternalEnergy:
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    constRefGasEThermoPhysics
+);
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    constRefFluidEThermoPhysics
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //