diff --git a/src/thermophysicalModels/basic/basicThermo/basicThermo.C b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
index e2fc874b5ef6ba930a1c85a745cf2032891c6e86..16ae3ebb8a5e73ea9c926899be792996797463e5 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.C
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.C
@@ -24,6 +24,14 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "basicThermo.H"
+#include "zeroGradientFvPatchFields.H"
+#include "fixedEnergyFvPatchScalarField.H"
+#include "gradientEnergyFvPatchScalarField.H"
+#include "mixedEnergyFvPatchScalarField.H"
+#include "fixedJumpFvPatchFields.H"
+#include "fixedJumpAMIFvPatchFields.H"
+#include "energyJumpFvPatchScalarField.H"
+#include "energyJumpAMIFvPatchScalarField.H"
 
 
 /* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
@@ -36,6 +44,84 @@ namespace Foam
 
 const Foam::word Foam::basicThermo::dictName("thermophysicalProperties");
 
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::wordList Foam::basicThermo::heBoundaryBaseTypes()
+{
+    const volScalarField::GeometricBoundaryField& tbf =
+        this->T_.boundaryField();
+
+    wordList hbt(tbf.size(), word::null);
+
+    forAll(tbf, patchi)
+    {
+        if (isA<fixedJumpFvPatchScalarField>(tbf[patchi]))
+        {
+            const fixedJumpFvPatchScalarField& pf =
+                dynamic_cast<const fixedJumpFvPatchScalarField&>(tbf[patchi]);
+
+            hbt[patchi] = pf.interfaceFieldType();
+        }
+        else if (isA<fixedJumpAMIFvPatchScalarField>(tbf[patchi]))
+        {
+            const fixedJumpAMIFvPatchScalarField& pf =
+                dynamic_cast<const fixedJumpAMIFvPatchScalarField&>
+                (
+                    tbf[patchi]
+                );
+
+            hbt[patchi] = pf.interfaceFieldType();
+        }
+    }
+
+    return hbt;
+}
+
+
+Foam::wordList Foam::basicThermo::heBoundaryTypes()
+{
+    const volScalarField::GeometricBoundaryField& tbf =
+        this->T_.boundaryField();
+
+    wordList hbt = tbf.types();
+
+    forAll(tbf, patchi)
+    {
+        if (isA<fixedValueFvPatchScalarField>(tbf[patchi]))
+        {
+            hbt[patchi] = fixedEnergyFvPatchScalarField::typeName;
+        }
+        else if
+        (
+            isA<zeroGradientFvPatchScalarField>(tbf[patchi])
+         || isA<fixedGradientFvPatchScalarField>(tbf[patchi])
+        )
+        {
+            hbt[patchi] = gradientEnergyFvPatchScalarField::typeName;
+        }
+        else if (isA<mixedFvPatchScalarField>(tbf[patchi]))
+        {
+            hbt[patchi] = mixedEnergyFvPatchScalarField::typeName;
+        }
+        else if (isA<fixedJumpFvPatchScalarField>(tbf[patchi]))
+        {
+            hbt[patchi] = energyJumpFvPatchScalarField::typeName;
+        }
+        else if (isA<fixedJumpAMIFvPatchScalarField>(tbf[patchi]))
+        {
+            hbt[patchi] = energyJumpAMIFvPatchScalarField::typeName;
+        }
+        else if (tbf[patchi].type() == "energyRegionCoupledFvPatchScalarField")
+        {
+            hbt[patchi] = "energyRegionCoupledFvPatchScalarField";
+        }
+    }
+
+    return hbt;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::volScalarField& Foam::basicThermo::lookupOrConstruct
@@ -395,6 +481,12 @@ const Foam::volScalarField& Foam::basicThermo::T() const
 }
 
 
+Foam::volScalarField& Foam::basicThermo::T()
+{
+    return T_;
+}
+
+
 const Foam::volScalarField& Foam::basicThermo::alpha() const
 {
     return alpha_;
diff --git a/src/thermophysicalModels/basic/basicThermo/basicThermo.H b/src/thermophysicalModels/basic/basicThermo/basicThermo.H
index e652c30ee6cc3216a0227160a9450e99e5223f78..c9de45836480da66f2c1f5537efe7ed0e3ac9602 100644
--- a/src/thermophysicalModels/basic/basicThermo/basicThermo.H
+++ b/src/thermophysicalModels/basic/basicThermo/basicThermo.H
@@ -77,6 +77,9 @@ protected:
         //- Should the dpdt term be included in the enthalpy equation
         Switch dpdt_;
 
+
+    // Protected Member Functions
+
         //- Construct as copy (not implemented)
         basicThermo(const basicThermo&);
 
@@ -86,6 +89,14 @@ protected:
             const char* name
         ) const;
 
+        //- Return the enthalpy/internal energy field boundary types
+        //  by interrogating the temperature field boundary types
+        wordList heBoundaryTypes();
+
+        //- Return the enthalpy/internal energy field boundary base types
+        //  by interrogating the temperature field boundary types
+        wordList heBoundaryBaseTypes();
+
 
 public:
 
@@ -314,6 +325,10 @@ public:
             //- Temperature [K]
             virtual const volScalarField& T() const;
 
+            //- Temperature [K]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& T();
+
             //- Heat capacity at constant pressure [J/kg/K]
             virtual tmp<volScalarField> Cp() const = 0;
 
diff --git a/src/thermophysicalModels/basic/heThermo/heThermo.C b/src/thermophysicalModels/basic/heThermo/heThermo.C
index 5f49a4ba752e03e483865019a6f751a273925262..bec062cc5ac33be33ee157c6d09159ca45fd5d09 100644
--- a/src/thermophysicalModels/basic/heThermo/heThermo.C
+++ b/src/thermophysicalModels/basic/heThermo/heThermo.C
@@ -24,94 +24,11 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "heThermo.H"
-#include "zeroGradientFvPatchFields.H"
-#include "fixedEnergyFvPatchScalarField.H"
 #include "gradientEnergyFvPatchScalarField.H"
 #include "mixedEnergyFvPatchScalarField.H"
-#include "fixedJumpFvPatchFields.H"
-#include "fixedJumpAMIFvPatchFields.H"
-#include "energyJumpFvPatchScalarField.H"
-#include "energyJumpAMIFvPatchScalarField.H"
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-template<class BasicThermo, class MixtureType>
-Foam::wordList Foam::heThermo<BasicThermo, MixtureType>::heBoundaryBaseTypes()
-{
-    const volScalarField::GeometricBoundaryField& tbf =
-        this->T_.boundaryField();
-
-    wordList hbt(tbf.size(), word::null);
-
-    forAll(tbf, patchi)
-    {
-        if (isA<fixedJumpFvPatchScalarField>(tbf[patchi]))
-        {
-            const fixedJumpFvPatchScalarField& pf =
-                dynamic_cast<const fixedJumpFvPatchScalarField&>(tbf[patchi]);
-
-            hbt[patchi] = pf.interfaceFieldType();
-        }
-        else if (isA<fixedJumpAMIFvPatchScalarField>(tbf[patchi]))
-        {
-            const fixedJumpAMIFvPatchScalarField& pf =
-                dynamic_cast<const fixedJumpAMIFvPatchScalarField&>
-                (
-                    tbf[patchi]
-                );
-
-            hbt[patchi] = pf.interfaceFieldType();
-        }
-    }
-
-    return hbt;
-}
-
-
-template<class BasicThermo, class MixtureType>
-Foam::wordList Foam::heThermo<BasicThermo, MixtureType>::heBoundaryTypes()
-{
-    const volScalarField::GeometricBoundaryField& tbf =
-        this->T_.boundaryField();
-
-    wordList hbt = tbf.types();
-
-    forAll(tbf, patchi)
-    {
-        if (isA<fixedValueFvPatchScalarField>(tbf[patchi]))
-        {
-            hbt[patchi] = fixedEnergyFvPatchScalarField::typeName;
-        }
-        else if
-        (
-            isA<zeroGradientFvPatchScalarField>(tbf[patchi])
-         || isA<fixedGradientFvPatchScalarField>(tbf[patchi])
-        )
-        {
-            hbt[patchi] = gradientEnergyFvPatchScalarField::typeName;
-        }
-        else if (isA<mixedFvPatchScalarField>(tbf[patchi]))
-        {
-            hbt[patchi] = mixedEnergyFvPatchScalarField::typeName;
-        }
-        else if (isA<fixedJumpFvPatchScalarField>(tbf[patchi]))
-        {
-            hbt[patchi] = energyJumpFvPatchScalarField::typeName;
-        }
-        else if (isA<fixedJumpAMIFvPatchScalarField>(tbf[patchi]))
-        {
-            hbt[patchi] = energyJumpAMIFvPatchScalarField::typeName;
-        }
-        else if (tbf[patchi].type() == "energyRegionCoupledFvPatchScalarField")
-        {
-            hbt[patchi] = "energyRegionCoupledFvPatchScalarField";
-        }
-    }
-
-    return hbt;
-}
-
-
 template<class BasicThermo, class MixtureType>
 void Foam::heThermo<BasicThermo, MixtureType>::
 heBoundaryCorrection(volScalarField& h)
diff --git a/src/thermophysicalModels/basic/heThermo/heThermo.H b/src/thermophysicalModels/basic/heThermo/heThermo.H
index f77d017e1396afd3840611449829eb89fbc6653a..c600258b65f4218a8b6a796fa0490aa84944c2ed 100644
--- a/src/thermophysicalModels/basic/heThermo/heThermo.H
+++ b/src/thermophysicalModels/basic/heThermo/heThermo.H
@@ -64,14 +64,6 @@ protected:
 
         // Enthalpy/Internal energy
 
-            //- Return the enthalpy/internal energy field boundary types
-            //  by interrogating the temperature field boundary types
-            wordList heBoundaryTypes();
-
-            //- Return the enthalpy/internal energy field boundary base types
-            //  by interrogating the temperature field boundary types
-            wordList heBoundaryBaseTypes();
-
             //- Correct the enthalpy/internal energy field boundaries
             void heBoundaryCorrection(volScalarField& he);