diff --git a/applications/utilities/thermophysical/adiabaticFlameT/adiabaticFlameT.C b/applications/utilities/thermophysical/adiabaticFlameT/adiabaticFlameT.C
index e6c26f52cb7c9f4d99e550e8cf54050802948323..cfdb984b6d0720fbde66437f6a89b22b9ec03f12 100644
--- a/applications/utilities/thermophysical/adiabaticFlameT/adiabaticFlameT.C
+++ b/applications/utilities/thermophysical/adiabaticFlameT/adiabaticFlameT.C
@@ -36,15 +36,16 @@ Description
 #include "IFstream.H"
 #include "OSspecific.H"
 
+#include "specie.H"
+#include "perfectGas.H"
 #include "thermo.H"
-#include "absoluteEnthalpy.H"
 #include "janafThermo.H"
-#include "perfectGas.H"
+#include "absoluteEnthalpy.H"
 
 using namespace Foam;
 
-typedef species::thermo<janafThermo<perfectGas>, absoluteEnthalpy> thermo;
-
+typedef species::thermo<janafThermo<perfectGas<specie> >, absoluteEnthalpy>
+    thermo;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/utilities/thermophysical/equilibriumCO/equilibriumCO.C b/applications/utilities/thermophysical/equilibriumCO/equilibriumCO.C
index 57f7aad30842a8fba7b22a97f5445c113dd2e02f..8bf51d72c2496af0040785d911134816c4a4b119 100644
--- a/applications/utilities/thermophysical/equilibriumCO/equilibriumCO.C
+++ b/applications/utilities/thermophysical/equilibriumCO/equilibriumCO.C
@@ -35,15 +35,18 @@ Description
 #include "OSspecific.H"
 #include "IOmanip.H"
 
+#include "specie.H"
+#include "perfectGas.H"
 #include "thermo.H"
-#include "absoluteEnthalpy.H"
 #include "janafThermo.H"
-#include "perfectGas.H"
+#include "absoluteEnthalpy.H"
+
 #include "SLPtrList.H"
 
 using namespace Foam;
 
-typedef species::thermo<janafThermo<perfectGas>, absoluteEnthalpy> thermo;
+typedef species::thermo<janafThermo<perfectGas<specie> >, absoluteEnthalpy>
+    thermo;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:
diff --git a/applications/utilities/thermophysical/equilibriumFlameT/equilibriumFlameT.C b/applications/utilities/thermophysical/equilibriumFlameT/equilibriumFlameT.C
index 3a11369da7c1d3af2722baaaab0085d0bd3089c9..40452749a462801f2c0a213e16ddd543e0c9d8b8 100644
--- a/applications/utilities/thermophysical/equilibriumFlameT/equilibriumFlameT.C
+++ b/applications/utilities/thermophysical/equilibriumFlameT/equilibriumFlameT.C
@@ -38,14 +38,16 @@ Description
 #include "OSspecific.H"
 #include "IOmanip.H"
 
+#include "specie.H"
+#include "perfectGas.H"
 #include "thermo.H"
-#include "absoluteEnthalpy.H"
 #include "janafThermo.H"
-#include "perfectGas.H"
+#include "absoluteEnthalpy.H"
 
 using namespace Foam;
 
-typedef species::thermo<janafThermo<perfectGas>, absoluteEnthalpy> thermo;
+typedef species::thermo<janafThermo<perfectGas<specie> >, absoluteEnthalpy>
+    thermo;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/utilities/thermophysical/mixtureAdiabaticFlameT/mixtureAdiabaticFlameT.C b/applications/utilities/thermophysical/mixtureAdiabaticFlameT/mixtureAdiabaticFlameT.C
index 0b00e865cac760089b60662b90830d22d8d8d327..7f1def4476e893fb562e8c5ef26875181237fe44 100644
--- a/applications/utilities/thermophysical/mixtureAdiabaticFlameT/mixtureAdiabaticFlameT.C
+++ b/applications/utilities/thermophysical/mixtureAdiabaticFlameT/mixtureAdiabaticFlameT.C
@@ -35,16 +35,17 @@ Description
 #include "IFstream.H"
 #include "OSspecific.H"
 
+#include "specie.H"
+#include "perfectGas.H"
 #include "thermo.H"
-#include "absoluteEnthalpy.H"
 #include "janafThermo.H"
-#include "perfectGas.H"
+#include "absoluteEnthalpy.H"
 #include "mixture.H"
 
 using namespace Foam;
 
-typedef species::thermo<janafThermo<perfectGas>, absoluteEnthalpy> thermo;
-
+typedef species::thermo<janafThermo<perfectGas<specie> >, absoluteEnthalpy>
+    thermo;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/thermophysicalModels/basic/fluidThermo/makeThermo.H b/src/thermophysicalModels/basic/fluidThermo/makeThermo.H
index abdd0a81bc996bf6e5db328ec63961b8162350ce..fde9d2e1a18d310e22a5ef954f9e838e268ee752 100644
--- a/src/thermophysicalModels/basic/fluidThermo/makeThermo.H
+++ b/src/thermophysicalModels/basic/fluidThermo/makeThermo.H
@@ -48,7 +48,7 @@ typedef Cthermo                                                               \
             <                                                                 \
                 Thermo                                                        \
                 <                                                             \
-                    EqnOfState                                                \
+                    EqnOfState<specie>                                        \
                 >,                                                            \
                 Type                                                          \
             >                                                                 \
@@ -97,7 +97,7 @@ typedef polynomialTransport                                                   \
     <                                                                         \
         hPolynomialThermo                                                     \
         <                                                                     \
-            icoPolynomial<Order>,                                             \
+            icoPolynomial<specie, Order>,                                     \
             Order                                                             \
         >,                                                                    \
         Type                                                                  \
diff --git a/src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C b/src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C
index a35994147f170d675c202f5453e9fc1e033007ff..2ce94472d0499fa4fa9f0176a40dc33b126fba6a 100644
--- a/src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C
+++ b/src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C
@@ -31,6 +31,7 @@ Description
 #include "basicMixture.H"
 #include "makeBasicMixture.H"
 
+#include "specie.H"
 #include "perfectGas.H"
 #include "rhoConst.H"
 #include "incompressiblePerfectGas.H"
diff --git a/src/thermophysicalModels/basic/mixtures/basicMixture/makeBasicMixture.H b/src/thermophysicalModels/basic/mixtures/basicMixture/makeBasicMixture.H
index 91fa87ea8df8157a2d8dffccf15e256499ada6bb..1997397646a04dc9409317a426a9d50b37b30858 100644
--- a/src/thermophysicalModels/basic/mixtures/basicMixture/makeBasicMixture.H
+++ b/src/thermophysicalModels/basic/mixtures/basicMixture/makeBasicMixture.H
@@ -39,7 +39,7 @@ Description
 #define makeBasicMixture(Mixture,Transport,Type,Thermo,EqnOfState)            \
                                                                               \
 typedef                                                                       \
-    Mixture<Transport<species::thermo<Thermo<EqnOfState>, Type> > >           \
+    Mixture<Transport<species::thermo<Thermo<EqnOfState<specie> >, Type> > >  \
     Mixture##Transport##Type##Thermo##EqnOfState;                             \
                                                                               \
 defineTemplateTypeNameAndDebugWithName                                        \
@@ -55,7 +55,7 @@ typedef polynomialTransport                                                   \
     <                                                                         \
         hPolynomialThermo                                                     \
         <                                                                     \
-            icoPolynomial<Order>,                                             \
+            icoPolynomial<specie, Order>,                                     \
             Order                                                             \
         >,                                                                    \
         Type                                                                  \
diff --git a/src/thermophysicalModels/basic/psiThermo/hePsiThermo/hePsiThermos.C b/src/thermophysicalModels/basic/psiThermo/hePsiThermo/hePsiThermos.C
index c997caab4b9b5adf801a3050db8a909aca1a77c0..5046cc85ec5853fd467192d2a91f399e736212c7 100644
--- a/src/thermophysicalModels/basic/psiThermo/hePsiThermo/hePsiThermos.C
+++ b/src/thermophysicalModels/basic/psiThermo/hePsiThermo/hePsiThermos.C
@@ -26,8 +26,8 @@ License
 #include "psiThermo.H"
 #include "makeThermo.H"
 
+#include "specie.H"
 #include "perfectGas.H"
-
 #include "hConstThermo.H"
 #include "eConstThermo.H"
 #include "janafThermo.H"
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
index 008aedd64d514918e11ab305dfd06b1355fb9149..a237983ad8e128617704324fa98e953c0116ee7a 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
@@ -619,7 +619,7 @@ bool finishReaction = false;
             currentSpecieName,
             new gasThermoPhysics
             (
-                janafThermo<perfectGas>
+                janafThermo<perfectGas<specie> >
                 (
                     specie
                     (
diff --git a/src/thermophysicalModels/reactionThermo/makeReactionThermo.H b/src/thermophysicalModels/reactionThermo/makeReactionThermo.H
index c0f38eb8f0c242285780c89bb18aad3cb2e3f07e..bbc8336dc6a9b2089afd431919c44901bc9a3267 100644
--- a/src/thermophysicalModels/reactionThermo/makeReactionThermo.H
+++ b/src/thermophysicalModels/reactionThermo/makeReactionThermo.H
@@ -45,7 +45,7 @@ typedef MixtureThermo                                                         \
                 <                                                             \
                     Thermo                                                    \
                     <                                                         \
-                        EqnOfState                                            \
+                        EqnOfState<specie>                                    \
                     >,                                                        \
                     Type                                                      \
                 >                                                             \
diff --git a/src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H b/src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H
index 7c5090a99a5adb3c0198b8d16542ccc825386a06..15b0a1c6eaf7e9993dd43d1522cebf70ead5b974 100644
--- a/src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H
+++ b/src/thermophysicalModels/solidSpecie/include/solidThermoPhysicsTypes.H
@@ -32,17 +32,14 @@ Description
 #ifndef solidThermoPhysicsTypes_H
 #define solidThermoPhysicsTypes_H
 
+#include "specie.H"
 #include "rhoConst.H"
-
 #include "hConstThermo.H"
 #include "hExponentialThermo.H"
-
 #include "constIsoSolidTransport.H"
 #include "constAnIsoSolidTransport.H"
 #include "exponentialSolidTransport.H"
-
 #include "constSolidRad.H"
-
 #include "sensibleEnthalpy.H"
 #include "thermo.H"
 
@@ -59,7 +56,7 @@ namespace Foam
                 <
                     hConstThermo
                     <
-                        rhoConst
+                        rhoConst<specie>
                     >,
                     sensibleEnthalpy
                 >
@@ -76,7 +73,7 @@ namespace Foam
                 <
                     hExponentialThermo
                     <
-                        rhoConst
+                        rhoConst<specie>
                     >,
                     sensibleEnthalpy
                 >
diff --git a/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H b/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
index ca08fb88e835e868613460330ab929255e0b9cfa..25aaf9d2eeeba9b6d2a57d4d7687cf1908a0bd1f 100644
--- a/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
+++ b/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
@@ -39,7 +39,7 @@ Description
 namespace Foam
 {
 
-#define makeSolidThermo(BaseThermo,Cthermo,Mixture,Transport,Radiation,Type,Thermo,Rho)\
+#define makeSolidThermo(BaseThermo,Cthermo,Mixture,Transport,Radiation,Type,Thermo,EqnOfState)\
                                                                               \
 typedef Cthermo                                                               \
 <                                                                             \
@@ -53,7 +53,7 @@ typedef Cthermo                                                               \
                 <                                                             \
                     Thermo                                                    \
                     <                                                         \
-                        Rho                                                   \
+                        EqnOfState<specie>                                    \
                     >,                                                        \
                     Type                                                      \
                 >                                                             \
@@ -61,11 +61,11 @@ typedef Cthermo                                                               \
         >                                                                     \
     >,                                                                        \
     BaseThermo                                                                \
->   Cthermo##Mixture##Transport##Radiation##Type##Thermo##Rho##BaseThermo;    \
+>Cthermo##Mixture##Transport##Radiation##Type##Thermo##EqnOfState##BaseThermo;\
                                                                               \
 defineTemplateTypeNameAndDebugWithName                                        \
 (                                                                             \
-    Cthermo##Mixture##Transport##Radiation##Type##Thermo##Rho##BaseThermo,    \
+Cthermo##Mixture##Transport##Radiation##Type##Thermo##EqnOfState##BaseThermo, \
     #Cthermo                                                                  \
     "<"                                                                       \
         #Mixture                                                              \
@@ -76,7 +76,7 @@ defineTemplateTypeNameAndDebugWithName                                        \
                 "<"                                                           \
                     #Thermo                                                   \
                     "<"                                                       \
-                        #Rho                                                  \
+                        #EqnOfState                                           \
                     ">,"                                                      \
                     #Type                                                     \
         ">>>>",                                                               \
@@ -94,19 +94,19 @@ typedef Mixture                                                               \
             <                                                                 \
                 Thermo                                                        \
                 <                                                             \
-                    Rho                                                       \
+                    EqnOfState<specie>                                        \
                 >,                                                            \
                 Type                                                          \
             >                                                                 \
         >                                                                     \
     >                                                                         \
->   Mixture##Transport##Radiation##Type##Thermo##Rho;                         \
+>   Mixture##Transport##Radiation##Type##Thermo##EqnOfState;                  \
                                                                               \
                                                                               \
                                                                               \
 defineTemplateTypeNameAndDebugWithName                                        \
 (                                                                             \
-    Mixture##Transport##Radiation##Type##Thermo##Rho,                         \
+    Mixture##Transport##Radiation##Type##Thermo##EqnOfState,                  \
     #Mixture                                                                  \
     "<"                                                                       \
         #Transport                                                            \
@@ -115,7 +115,7 @@ defineTemplateTypeNameAndDebugWithName                                        \
             "<"                                                               \
                 #Thermo                                                       \
                 "<"                                                           \
-                    #Rho                                                      \
+                    #EqnOfState                                               \
                 ">,"                                                          \
                 #Type                                                         \
     ">>>",                                                                    \
@@ -125,20 +125,20 @@ defineTemplateTypeNameAndDebugWithName                                        \
 addToRunTimeSelectionTable                                                    \
 (                                                                             \
     BaseThermo,                                                               \
-    Cthermo##Mixture##Transport##Radiation##Type##Thermo##Rho##BaseThermo,    \
+Cthermo##Mixture##Transport##Radiation##Type##Thermo##EqnOfState##BaseThermo, \
     mesh                                                                      \
 );                                                                            \
 addToRunTimeSelectionTable                                                    \
 (                                                                             \
     basicThermo,                                                              \
-    Cthermo##Mixture##Transport##Radiation##Type##Thermo##Rho##BaseThermo,    \
+Cthermo##Mixture##Transport##Radiation##Type##Thermo##EqnOfState##BaseThermo, \
     fvMesh                                                                    \
 );                                                                            \
                                                                               \
 addToRunTimeSelectionTable                                                    \
 (                                                                             \
     BaseThermo,                                                               \
-    Cthermo##Mixture##Transport##Radiation##Type##Thermo##Rho##BaseThermo,    \
+Cthermo##Mixture##Transport##Radiation##Type##Thermo##EqnOfState##BaseThermo, \
     dictionary                                                                \
 );
 
diff --git a/src/thermophysicalModels/specie/Make/files b/src/thermophysicalModels/specie/Make/files
index 23a7bdf3a9e9482cb61403f2274b656da87dd660..739c9f91fd0b2b717bbc39e4a19af73887bf052a 100644
--- a/src/thermophysicalModels/specie/Make/files
+++ b/src/thermophysicalModels/specie/Make/files
@@ -5,9 +5,6 @@ reactions = reaction/reactions
 
 $(atomicWeights)/atomicWeights.C
 $(specie)/specie.C
-$(equationOfState)/perfectGas/perfectGas.C
-$(equationOfState)/rhoConst/rhoConst.C
-$(equationOfState)/incompressiblePerfectGas/incompressiblePerfectGas.C
 $(reactions)/makeReactionThermoReactions.C
 $(reactions)/makeLangmuirHinshelwoodReactions.C
 
diff --git a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.C b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.C
index 80b3cdca64d231ee948582510105845219943be7..01b417e0dbed2eb556e1e1dc8fcc6816e2caebcc 100644
--- a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.C
+++ b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.C
@@ -33,20 +33,20 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<int PolySize>
-icoPolynomial<PolySize>::icoPolynomial(Istream& is)
+template<class Specie, int PolySize>
+icoPolynomial<Specie, PolySize>::icoPolynomial(Istream& is)
 :
-    specie(is),
+    Specie(is),
     rhoCoeffs_("rhoCoeffs<" + Foam::name(PolySize) + '>', is)
 {
     rhoCoeffs_ *= this->W();
 }
 
 
-template<int PolySize>
-icoPolynomial<PolySize>::icoPolynomial(const dictionary& dict)
+template<class Specie, int PolySize>
+icoPolynomial<Specie, PolySize>::icoPolynomial(const dictionary& dict)
 :
-    specie(dict),
+    Specie(dict),
     rhoCoeffs_
 (
     dict.subDict("equationOfState").lookup
@@ -61,10 +61,10 @@ icoPolynomial<PolySize>::icoPolynomial(const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<int PolySize>
-void icoPolynomial<PolySize>::write(Ostream& os) const
+template<class Specie, int PolySize>
+void icoPolynomial<Specie, PolySize>::write(Ostream& os) const
 {
-    specie::write(os);
+    Specie::write(os);
 
     dictionary dict("equationOfState");
     dict.add
@@ -79,16 +79,17 @@ void icoPolynomial<PolySize>::write(Ostream& os) const
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
-template<int PolySize>
-Ostream& operator<<(Ostream& os, const icoPolynomial<PolySize>& ip)
+template<class Specie, int PolySize>
+Ostream& operator<<(Ostream& os, const icoPolynomial<Specie, PolySize>& ip)
 {
-    os  << static_cast<const specie&>(ip) << tab
+    os  << static_cast<const Specie&>(ip) << tab
         << "rhoCoeffs<" << Foam::name(PolySize) << '>' << tab
         << ip.rhoCoeffs_/ip.W();
 
     os.check
     (
-        "Ostream& operator<<(Ostream& os, const icoPolynomial<PolySize>& ip)"
+        "Ostream& operator<<"
+        "(Ostream& os, const icoPolynomial<Specie, PolySize>& ip)"
     );
 
     return os;
diff --git a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H
index 3afc213e227f76894b6977051e9359eef5ff7e41..ba92455f8276f9cf6c44ca8a4b2b3933372654de 100644
--- a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H
+++ b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H
@@ -37,7 +37,6 @@ SourceFiles
 #ifndef icoPolynomial_H
 #define icoPolynomial_H
 
-#include "specie.H"
 #include "autoPtr.H"
 #include "Polynomial.H"
 
@@ -48,42 +47,42 @@ namespace Foam
 
 // Forward declaration of friend functions and operators
 
-template<int PolySize>
+template<class Specie, int PolySize>
 class icoPolynomial;
 
-template<int PolySize>
-icoPolynomial<PolySize> operator+
+template<class Specie, int PolySize>
+icoPolynomial<Specie, PolySize> operator+
 (
-    const icoPolynomial<PolySize>&,
-    const icoPolynomial<PolySize>&
+    const icoPolynomial<Specie, PolySize>&,
+    const icoPolynomial<Specie, PolySize>&
 );
 
-template<int PolySize>
-icoPolynomial<PolySize> operator-
+template<class Specie, int PolySize>
+icoPolynomial<Specie, PolySize> operator-
 (
-    const icoPolynomial<PolySize>&,
-    const icoPolynomial<PolySize>&
+    const icoPolynomial<Specie, PolySize>&,
+    const icoPolynomial<Specie, PolySize>&
 );
 
-template<int PolySize>
-icoPolynomial<PolySize> operator*
+template<class Specie, int PolySize>
+icoPolynomial<Specie, PolySize> operator*
 (
     const scalar,
-    const icoPolynomial<PolySize>&
+    const icoPolynomial<Specie, PolySize>&
 );
 
-template<int PolySize>
-icoPolynomial<PolySize> operator==
+template<class Specie, int PolySize>
+icoPolynomial<Specie, PolySize> operator==
 (
-    const icoPolynomial<PolySize>&,
-    const icoPolynomial<PolySize>&
+    const icoPolynomial<Specie, PolySize>&,
+    const icoPolynomial<Specie, PolySize>&
 );
 
-template<int PolySize>
+template<class Specie, int PolySize>
 Ostream& operator<<
 (
     Ostream&,
-    const icoPolynomial<PolySize>&
+    const icoPolynomial<Specie, PolySize>&
 );
 
 
@@ -91,10 +90,10 @@ Ostream& operator<<
                         Class icoPolynomial Declaration
 \*---------------------------------------------------------------------------*/
 
-template<int PolySize>
+template<class Specie, int PolySize>
 class icoPolynomial
 :
-    public specie
+    public Specie
 {
     // Private data
 
@@ -110,7 +109,7 @@ public:
         //- Construct from components
         inline icoPolynomial
         (
-            const specie& sp,
+            const Specie& sp,
             const Polynomial<PolySize>& rhoPoly
         );
 
@@ -176,25 +175,25 @@ public:
 
     // Friend operators
 
-        friend icoPolynomial operator+ <PolySize>
+        friend icoPolynomial operator+ <Specie, PolySize>
         (
             const icoPolynomial&,
             const icoPolynomial&
         );
 
-        friend icoPolynomial operator- <PolySize>
+        friend icoPolynomial operator- <Specie, PolySize>
         (
             const icoPolynomial&,
             const icoPolynomial&
         );
 
-        friend icoPolynomial operator* <PolySize>
+        friend icoPolynomial operator* <Specie, PolySize>
         (
             const scalar s,
             const icoPolynomial&
         );
 
-        friend icoPolynomial operator== <PolySize>
+        friend icoPolynomial operator== <Specie, PolySize>
         (
             const icoPolynomial&,
             const icoPolynomial&
@@ -203,7 +202,11 @@ public:
 
     // Ostream Operator
 
-        friend Ostream& operator<< <PolySize>(Ostream&, const icoPolynomial&);
+        friend Ostream& operator<< <Specie, PolySize>
+        (
+            Ostream&,
+            const icoPolynomial&
+        );
 };
 
 
@@ -217,7 +220,7 @@ public:
                                                                              \
 defineTemplateTypeNameAndDebugWithName                                       \
 (                                                                            \
-    icoPolynomial<PolySize>,                                                 \
+    icoPolynomial<Specie, PolySize>,                                         \
     "icoPolynomial<"#PolySize">",                                            \
     0                                                                        \
 );
diff --git a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H
index cfa527bd0ea491638b394eb0c7a78af0c72af553..dbe3772ec0deffc7a073c56303a9e2864280bf7c 100644
--- a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H
+++ b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H
@@ -27,98 +27,113 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<int PolySize>
-inline Foam::icoPolynomial<PolySize>::icoPolynomial
+template<class Specie, int PolySize>
+inline Foam::icoPolynomial<Specie, PolySize>::icoPolynomial
 (
-    const specie& sp,
+    const Specie& sp,
     const Polynomial<PolySize>& rhoCoeffs
 )
 :
-    specie(sp),
+    Specie(sp),
     rhoCoeffs_(rhoCoeffs)
 {}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<int PolySize>
-inline Foam::icoPolynomial<PolySize>::icoPolynomial
+template<class Specie, int PolySize>
+inline Foam::icoPolynomial<Specie, PolySize>::icoPolynomial
 (
-    const icoPolynomial<PolySize>& ip
+    const icoPolynomial<Specie, PolySize>& ip
 )
 :
-    specie(ip),
+    Specie(ip),
     rhoCoeffs_(ip.rhoCoeffs_)
 {}
 
 
-template<int PolySize>
-inline Foam::icoPolynomial<PolySize>::icoPolynomial
+template<class Specie, int PolySize>
+inline Foam::icoPolynomial<Specie, PolySize>::icoPolynomial
 (
     const word& name,
-    const icoPolynomial<PolySize>& ip
+    const icoPolynomial<Specie, PolySize>& ip
 )
 :
-    specie(name, ip),
+    Specie(name, ip),
     rhoCoeffs_(ip.rhoCoeffs_)
 {}
 
 
-template<int PolySize>
-inline Foam::autoPtr<Foam::icoPolynomial<PolySize> >
-Foam::icoPolynomial<PolySize>::clone() const
+template<class Specie, int PolySize>
+inline Foam::autoPtr<Foam::icoPolynomial<Specie, PolySize> >
+Foam::icoPolynomial<Specie, PolySize>::clone() const
 {
-    return autoPtr<icoPolynomial<PolySize> >
+    return autoPtr<icoPolynomial<Specie, PolySize> >
     (
-        new icoPolynomial<PolySize>(*this)
+        new icoPolynomial<Specie, PolySize>(*this)
     );
 }
 
 
-template<int PolySize>
-inline Foam::autoPtr<Foam::icoPolynomial<PolySize> >
-Foam::icoPolynomial<PolySize>::New(Istream& is)
+template<class Specie, int PolySize>
+inline Foam::autoPtr<Foam::icoPolynomial<Specie, PolySize> >
+Foam::icoPolynomial<Specie, PolySize>::New(Istream& is)
 {
-    return autoPtr<icoPolynomial<PolySize> >(new icoPolynomial<PolySize>(is));
+    return autoPtr<icoPolynomial<Specie, PolySize> >
+    (
+        new icoPolynomial<Specie, PolySize>(is)
+    );
 }
 
 
-template<int PolySize>
-inline Foam::autoPtr<Foam::icoPolynomial<PolySize> >
-Foam::icoPolynomial<PolySize>::New(const dictionary& dict)
+template<class Specie, int PolySize>
+inline Foam::autoPtr<Foam::icoPolynomial<Specie, PolySize> >
+Foam::icoPolynomial<Specie, PolySize>::New(const dictionary& dict)
 {
-    return autoPtr<icoPolynomial<PolySize> >
+    return autoPtr<icoPolynomial<Specie, PolySize> >
     (
-        new icoPolynomial<PolySize>(dict)
+        new icoPolynomial<Specie, PolySize>(dict)
     );
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<int PolySize>
-inline Foam::scalar Foam::icoPolynomial<PolySize>::rho(scalar, scalar T) const
+template<class Specie, int PolySize>
+inline Foam::scalar Foam::icoPolynomial<Specie, PolySize>::rho
+(
+    scalar,
+    scalar T
+) const
 {
     return rhoCoeffs_.value(T)/this->W();
 }
 
 
-template<int PolySize>
-inline Foam::scalar Foam::icoPolynomial<PolySize>::psi(scalar, scalar) const
+template<class Specie, int PolySize>
+inline Foam::scalar Foam::icoPolynomial<Specie, PolySize>::psi
+(
+    scalar,
+    scalar
+) const
 {
     return 0.0;
 }
 
 
-template<int PolySize>
-inline Foam::scalar Foam::icoPolynomial<PolySize>::Z(scalar, scalar) const
+template<class Specie, int PolySize>
+inline Foam::scalar Foam::icoPolynomial<Specie, PolySize>::Z
+(
+    scalar,
+    scalar
+) const
 {
     return 0.0;
 }
 
 
-template<int PolySize>
-inline Foam::scalar Foam::icoPolynomial<PolySize>::cpMcv
+template<class Specie, int PolySize>
+inline Foam::scalar Foam::icoPolynomial<Specie, PolySize>::cpMcv
 (
     scalar p,
     scalar T
@@ -130,13 +145,14 @@ inline Foam::scalar Foam::icoPolynomial<PolySize>::cpMcv
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-template<int PolySize>
-inline Foam::icoPolynomial<PolySize>& Foam::icoPolynomial<PolySize>::operator=
+template<class Specie, int PolySize>
+inline Foam::icoPolynomial<Specie, PolySize>&
+Foam::icoPolynomial<Specie, PolySize>::operator=
 (
-    const icoPolynomial<PolySize>& ip
+    const icoPolynomial<Specie, PolySize>& ip
 )
 {
-    specie::operator=(ip);
+    Specie::operator=(ip);
 
     rhoCoeffs_ = ip.rhoCoeffs_;
 
@@ -144,15 +160,15 @@ inline Foam::icoPolynomial<PolySize>& Foam::icoPolynomial<PolySize>::operator=
 }
 
 
-template<int PolySize>
-inline void Foam::icoPolynomial<PolySize>::operator+=
+template<class Specie, int PolySize>
+inline void Foam::icoPolynomial<Specie, PolySize>::operator+=
 (
-    const icoPolynomial<PolySize>& ip
+    const icoPolynomial<Specie, PolySize>& ip
 )
 {
     scalar molr1 = this->nMoles();
 
-    specie::operator+=(ip);
+    Specie::operator+=(ip);
 
     molr1 /= this->nMoles();
     scalar molr2 = ip.nMoles()/this->nMoles();
@@ -161,15 +177,15 @@ inline void Foam::icoPolynomial<PolySize>::operator+=
 }
 
 
-template<int PolySize>
-inline void Foam::icoPolynomial<PolySize>::operator-=
+template<class Specie, int PolySize>
+inline void Foam::icoPolynomial<Specie, PolySize>::operator-=
 (
-    const icoPolynomial<PolySize>& ip
+    const icoPolynomial<Specie, PolySize>& ip
 )
 {
     scalar molr1 = this->nMoles();
 
-    specie::operator-=(ip);
+    Specie::operator-=(ip);
 
     molr1 /= this->nMoles();
     scalar molr2 = ip.nMoles()/this->nMoles();
@@ -178,75 +194,75 @@ inline void Foam::icoPolynomial<PolySize>::operator-=
 }
 
 
-template<int PolySize>
-inline void Foam::icoPolynomial<PolySize>::operator*=(const scalar s)
+template<class Specie, int PolySize>
+inline void Foam::icoPolynomial<Specie, PolySize>::operator*=(const scalar s)
 {
-    specie::operator*=(s);
+    Specie::operator*=(s);
 }
 
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-template<int PolySize>
-Foam::icoPolynomial<PolySize> Foam::operator+
+template<class Specie, int PolySize>
+Foam::icoPolynomial<Specie, PolySize> Foam::operator+
 (
-    const icoPolynomial<PolySize>& ip1,
-    const icoPolynomial<PolySize>& ip2
+    const icoPolynomial<Specie, PolySize>& ip1,
+    const icoPolynomial<Specie, PolySize>& ip2
 )
 {
     scalar nMoles = ip1.nMoles() + ip2.nMoles();
     scalar molr1 = ip1.nMoles()/nMoles;
     scalar molr2 = ip2.nMoles()/nMoles;
 
-    return icoPolynomial<PolySize>
+    return icoPolynomial<Specie, PolySize>
     (
-        static_cast<const specie&>(ip1)
-      + static_cast<const specie&>(ip2),
+        static_cast<const Specie&>(ip1)
+      + static_cast<const Specie&>(ip2),
         molr1*ip1.rhoCoeffs_ + molr2*ip2.rhoCoeffs_
     );
 }
 
 
-template<int PolySize>
-Foam::icoPolynomial<PolySize> Foam::operator-
+template<class Specie, int PolySize>
+Foam::icoPolynomial<Specie, PolySize> Foam::operator-
 (
-    const icoPolynomial<PolySize>& ip1,
-    const icoPolynomial<PolySize>& ip2
+    const icoPolynomial<Specie, PolySize>& ip1,
+    const icoPolynomial<Specie, PolySize>& ip2
 )
 {
     scalar nMoles = ip1.nMoles() + ip2.nMoles();
     scalar molr1 = ip1.nMoles()/nMoles;
     scalar molr2 = ip2.nMoles()/nMoles;
 
-    return icoPolynomial<PolySize>
+    return icoPolynomial<Specie, PolySize>
     (
-        static_cast<const specie&>(ip1)
-      - static_cast<const specie&>(ip2),
+        static_cast<const Specie&>(ip1)
+      - static_cast<const Specie&>(ip2),
         molr1*ip1.rhoCoeffs_ - molr2*ip2.rhoCoeffs_
     );
 }
 
 
-template<int PolySize>
-Foam::icoPolynomial<PolySize> Foam::operator*
+template<class Specie, int PolySize>
+Foam::icoPolynomial<Specie, PolySize> Foam::operator*
 (
     const scalar s,
-    const icoPolynomial<PolySize>& ip
+    const icoPolynomial<Specie, PolySize>& ip
 )
 {
-    return icoPolynomial<PolySize>
+    return icoPolynomial<Specie, PolySize>
     (
-        s*static_cast<const specie&>(ip),
+        s*static_cast<const Specie&>(ip),
         ip.rhoCoeffs_
     );
 }
 
 
-template<int PolySize>
-Foam::icoPolynomial<PolySize> Foam::operator==
+template<class Specie, int PolySize>
+Foam::icoPolynomial<Specie, PolySize> Foam::operator==
 (
-    const icoPolynomial<PolySize>& ip1,
-    const icoPolynomial<PolySize>& ip2
+    const icoPolynomial<Specie, PolySize>& ip1,
+    const icoPolynomial<Specie, PolySize>& ip2
 )
 {
     return ip2 - ip1;
diff --git a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.C b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.C
index 700980e44fd3a843132481262f0f46e720c7cbab..4f957aafedfaa2995d107e4e7a939e8a86ffb3fe 100644
--- a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.C
+++ b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.C
@@ -28,27 +28,37 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::incompressiblePerfectGas::incompressiblePerfectGas(Istream& is)
+template<class Specie>
+Foam::incompressiblePerfectGas<Specie>::incompressiblePerfectGas(Istream& is)
 :
-    specie(is),
+    Specie(is),
     pRef_(readScalar(is))
 {
-    is.check("incompressiblePerfectGas::incompressiblePerfectGas(Istream& is)");
+    is.check
+    (
+        "incompressiblePerfectGas<Specie>::"
+        "incompressiblePerfectGas(Istream& is)"
+    );
 }
 
 
-Foam::incompressiblePerfectGas::incompressiblePerfectGas(const dictionary& dict)
+template<class Specie>
+Foam::incompressiblePerfectGas<Specie>::incompressiblePerfectGas
+(
+    const dictionary& dict
+)
 :
-    specie(dict),
+    Specie(dict),
     pRef_(readScalar(dict.subDict("equationOfState").lookup("pRef")))
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::incompressiblePerfectGas::write(Ostream& os) const
+template<class Specie>
+void Foam::incompressiblePerfectGas<Specie>::write(Ostream& os) const
 {
-    specie::write(os);
+    Specie::write(os);
     dictionary dict("equationOfState");
     dict.add("pRef", pRef_);
 
@@ -58,14 +68,20 @@ void Foam::incompressiblePerfectGas::write(Ostream& os) const
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
-Foam::Ostream& Foam::operator<<(Ostream& os, const incompressiblePerfectGas& pg)
+template<class Specie>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const incompressiblePerfectGas<Specie>& pg
+)
 {
-    os  << static_cast<const specie&>(pg)
+    os  << static_cast<const Specie&>(pg)
         << token::SPACE << pg.pRef_;
 
     os.check
     (
-        "Ostream& operator<<(Ostream& os, const incompressiblePerfectGas& st)"
+        "Ostream& operator<<"
+        "(Ostream& os, const incompressiblePerfectGas<Specie>& st)"
     );
     return os;
 }
diff --git a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H
index c755c2d8f212ce23621d261b9feb78b633ec46e5..a752cc924a042eb78d44ba99b3fec16e122e3474 100644
--- a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H
+++ b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGas.H
@@ -38,7 +38,6 @@ SourceFiles
 #ifndef incompressiblePerfectGas_H
 #define incompressiblePerfectGas_H
 
-#include "specie.H"
 #include "autoPtr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -46,13 +45,54 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+
+template<class Specie> class incompressiblePerfectGas;
+
+template<class Specie>
+inline incompressiblePerfectGas<Specie> operator+
+(
+    const incompressiblePerfectGas<Specie>&,
+    const incompressiblePerfectGas<Specie>&
+);
+
+template<class Specie>
+inline incompressiblePerfectGas<Specie> operator-
+(
+    const incompressiblePerfectGas<Specie>&,
+    const incompressiblePerfectGas<Specie>&
+);
+
+template<class Specie>
+inline incompressiblePerfectGas<Specie> operator*
+(
+    const scalar,
+    const incompressiblePerfectGas<Specie>&
+);
+
+template<class Specie>
+inline incompressiblePerfectGas<Specie> operator==
+(
+    const incompressiblePerfectGas<Specie>&,
+    const incompressiblePerfectGas<Specie>&
+);
+
+template<class Specie>
+Ostream& operator<<
+(
+    Ostream&,
+    const incompressiblePerfectGas<Specie>&
+);
+
+
 /*---------------------------------------------------------------------------*\
                     Class incompressiblePerfectGas Declaration
 \*---------------------------------------------------------------------------*/
 
+template<class Specie>
 class incompressiblePerfectGas
 :
-    public specie
+    public Specie
 {
     // Private data
 
@@ -65,7 +105,7 @@ public:
     // Constructors
 
         //- Construct from components
-        inline incompressiblePerfectGas(const specie& sp);
+        inline incompressiblePerfectGas(const Specie& sp);
 
         //- Construct from Istream
         incompressiblePerfectGas(Istream&);
@@ -132,25 +172,25 @@ public:
 
     // Friend operators
 
-        inline friend incompressiblePerfectGas operator+
+        inline friend incompressiblePerfectGas operator+ <Specie>
         (
             const incompressiblePerfectGas&,
             const incompressiblePerfectGas&
         );
 
-        inline friend incompressiblePerfectGas operator-
+        inline friend incompressiblePerfectGas operator- <Specie>
         (
             const incompressiblePerfectGas&,
             const incompressiblePerfectGas&
         );
 
-        inline friend incompressiblePerfectGas operator*
+        inline friend incompressiblePerfectGas operator* <Specie>
         (
             const scalar s,
             const incompressiblePerfectGas&
         );
 
-        inline friend incompressiblePerfectGas operator==
+        inline friend incompressiblePerfectGas operator== <Specie>
         (
             const incompressiblePerfectGas&,
             const incompressiblePerfectGas&
@@ -159,7 +199,11 @@ public:
 
     // Ostream Operator
 
-        friend Ostream& operator<<(Ostream&, const incompressiblePerfectGas&);
+        friend Ostream& operator<< <Specie>
+        (
+            Ostream&,
+            const incompressiblePerfectGas&
+        );
 };
 
 
@@ -171,6 +215,10 @@ public:
 
 #include "incompressiblePerfectGasI.H"
 
+#ifdef NoRepository
+#   include "incompressiblePerfectGas.C"
+#endif
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H
index d179d5e2a181c10d78880e3426295ccbd67fd7c8..340d4b1e4c5f502320282a3256228329ee53abe2 100644
--- a/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H
+++ b/src/thermophysicalModels/specie/equationOfState/incompressiblePerfectGas/incompressiblePerfectGasI.H
@@ -27,85 +27,109 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-inline Foam::incompressiblePerfectGas::incompressiblePerfectGas
+template<class Specie>
+inline Foam::incompressiblePerfectGas<Specie>::incompressiblePerfectGas
 (
-    const specie& sp
+    const Specie& sp
 )
 :
-    specie(sp)
+    Specie(sp)
 {}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-inline Foam::incompressiblePerfectGas::incompressiblePerfectGas
+template<class Specie>
+inline Foam::incompressiblePerfectGas<Specie>::incompressiblePerfectGas
 (
     const word& name,
-    const incompressiblePerfectGas& pg
+    const incompressiblePerfectGas<Specie>& pg
 )
 :
-    specie(name, pg)
+    Specie(name, pg)
 {}
 
 
-inline Foam::autoPtr<Foam::incompressiblePerfectGas>
-Foam::incompressiblePerfectGas::clone() const
+template<class Specie>
+inline Foam::autoPtr<Foam::incompressiblePerfectGas<Specie> >
+Foam::incompressiblePerfectGas<Specie>::clone() const
 {
-    return autoPtr<incompressiblePerfectGas>
+    return autoPtr<incompressiblePerfectGas<Specie> >
     (
-        new incompressiblePerfectGas(*this)
+        new incompressiblePerfectGas<Specie>(*this)
     );
 }
 
 
-inline Foam::autoPtr<Foam::incompressiblePerfectGas>
-Foam::incompressiblePerfectGas::New
+template<class Specie>
+inline Foam::autoPtr<Foam::incompressiblePerfectGas<Specie> >
+Foam::incompressiblePerfectGas<Specie>::New
 (
     Istream& is
 )
 {
-    return autoPtr<incompressiblePerfectGas>(new incompressiblePerfectGas(is));
+    return autoPtr<incompressiblePerfectGas<Specie> >
+    (
+        new incompressiblePerfectGas<Specie>(is)
+    );
 }
 
 
-inline Foam::autoPtr<Foam::incompressiblePerfectGas>
-Foam::incompressiblePerfectGas::New
+template<class Specie>
+inline Foam::autoPtr<Foam::incompressiblePerfectGas<Specie> >
+Foam::incompressiblePerfectGas<Specie>::New
 (
     const dictionary& dict
 )
 {
-    return autoPtr<incompressiblePerfectGas>
+    return autoPtr<incompressiblePerfectGas<Specie> >
     (
-        new incompressiblePerfectGas(dict)
+        new incompressiblePerfectGas<Specie>(dict)
     );
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::scalar Foam::incompressiblePerfectGas::rho
+template<class Specie>
+inline Foam::scalar Foam::incompressiblePerfectGas<Specie>::rho
 (
     scalar p,
     scalar T
 ) const
 {
-    return pRef_/(R()*T);
+    return pRef_/(this->R()*T);
 }
 
 
-inline Foam::scalar Foam::incompressiblePerfectGas::psi(scalar, scalar T) const
+template<class Specie>
+inline Foam::scalar Foam::incompressiblePerfectGas<Specie>::psi
+(
+    scalar,
+    scalar T
+) const
 {
     return 0.0;
 }
 
 
-inline Foam::scalar Foam::incompressiblePerfectGas::Z(scalar, scalar) const
+template<class Specie>
+inline Foam::scalar Foam::incompressiblePerfectGas<Specie>::Z
+(
+    scalar,
+    scalar
+) const
 {
     return 0.0;
 }
 
 
-inline Foam::scalar Foam::incompressiblePerfectGas::cpMcv(scalar, scalar) const
+template<class Specie>
+inline Foam::scalar Foam::incompressiblePerfectGas<Specie>::cpMcv
+(
+    scalar,
+    scalar
+) const
 {
     return this->RR;
 }
@@ -114,74 +138,81 @@ inline Foam::scalar Foam::incompressiblePerfectGas::cpMcv(scalar, scalar) const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-inline void Foam::incompressiblePerfectGas::operator+=
+template<class Specie>
+inline void Foam::incompressiblePerfectGas<Specie>::operator+=
 (
-    const incompressiblePerfectGas& pg
+    const incompressiblePerfectGas<Specie>& pg
 )
 {
-    specie::operator+=(pg);
+    Specie::operator+=(pg);
 }
 
 
-inline void Foam::incompressiblePerfectGas::operator-=
+template<class Specie>
+inline void Foam::incompressiblePerfectGas<Specie>::operator-=
 (
-    const incompressiblePerfectGas& pg
+    const incompressiblePerfectGas<Specie>& pg
 )
 {
-    specie::operator-=(pg);
+    Specie::operator-=(pg);
 }
 
 
-inline void Foam::incompressiblePerfectGas::operator*=(const scalar s)
+template<class Specie>
+inline void Foam::incompressiblePerfectGas<Specie>::operator*=(const scalar s)
 {
-    specie::operator*=(s);
+    Specie::operator*=(s);
 }
 
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-inline Foam::incompressiblePerfectGas Foam::operator+
+template<class Specie>
+inline Foam::incompressiblePerfectGas<Specie> Foam::operator+
 (
-    const incompressiblePerfectGas& pg1,
-    const incompressiblePerfectGas& pg2
+    const incompressiblePerfectGas<Specie>& pg1,
+    const incompressiblePerfectGas<Specie>& pg2
 )
 {
-    return incompressiblePerfectGas
+    return incompressiblePerfectGas<Specie>
     (
-        static_cast<const specie&>(pg1)
-      + static_cast<const specie&>(pg2)
+        static_cast<const Specie&>(pg1)
+      + static_cast<const Specie&>(pg2)
     );
 }
 
 
-inline Foam::incompressiblePerfectGas Foam::operator-
+template<class Specie>
+inline Foam::incompressiblePerfectGas<Specie> Foam::operator-
 (
-    const incompressiblePerfectGas& pg1,
-    const incompressiblePerfectGas& pg2
+    const incompressiblePerfectGas<Specie>& pg1,
+    const incompressiblePerfectGas<Specie>& pg2
 )
 {
-    return incompressiblePerfectGas
+    return incompressiblePerfectGas<Specie>
     (
-        static_cast<const specie&>(pg1)
-      - static_cast<const specie&>(pg2)
+        static_cast<const Specie&>(pg1)
+      - static_cast<const Specie&>(pg2)
     );
 }
 
 
-inline Foam::incompressiblePerfectGas Foam::operator*
+template<class Specie>
+inline Foam::incompressiblePerfectGas<Specie> Foam::operator*
 (
     const scalar s,
-    const incompressiblePerfectGas& pg
+    const incompressiblePerfectGas<Specie>& pg
 )
 {
-    return incompressiblePerfectGas(s*static_cast<const specie&>(pg));
+    return incompressiblePerfectGas<Specie>(s*static_cast<const Specie&>(pg));
 }
 
 
-inline Foam::incompressiblePerfectGas Foam::operator==
+template<class Specie>
+inline Foam::incompressiblePerfectGas<Specie> Foam::operator==
 (
-    const incompressiblePerfectGas& pg1,
-    const incompressiblePerfectGas& pg2
+    const incompressiblePerfectGas<Specie>& pg1,
+    const incompressiblePerfectGas<Specie>& pg2
 )
 {
     return pg2 - pg1;
diff --git a/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.C b/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.C
index d575368ce5ecb2272fd5d40f362408d513f709f8..5e4cf38b7299286ce9f545cbcca1c0b8096f0780 100644
--- a/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.C
+++ b/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.C
@@ -28,35 +28,39 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::perfectGas::perfectGas(Istream& is)
+template<class Specie>
+Foam::perfectGas<Specie>::perfectGas(Istream& is)
 :
-    specie(is)
+    Specie(is)
 {
-    is.check("perfectGas::perfectGas(Istream& is)");
+    is.check("perfectGas<Specie>::perfectGas(Istream& is)");
 }
 
 
-Foam::perfectGas::perfectGas(const dictionary& dict)
+template<class Specie>
+Foam::perfectGas<Specie>::perfectGas(const dictionary& dict)
 :
-    specie(dict)
+    Specie(dict)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::perfectGas::write(Ostream& os) const
+template<class Specie>
+void Foam::perfectGas<Specie>::write(Ostream& os) const
 {
-    specie::write(os);
+    Specie::write(os);
 }
 
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
-Foam::Ostream& Foam::operator<<(Ostream& os, const perfectGas& pg)
+template<class Specie>
+Foam::Ostream& Foam::operator<<(Ostream& os, const perfectGas<Specie>& pg)
 {
-    os  << static_cast<const specie&>(pg);
+    os  << static_cast<const Specie&>(pg);
 
-    os.check("Ostream& operator<<(Ostream& os, const perfectGas& st)");
+    os.check("Ostream& operator<<(Ostream& os, const perfectGas<Specie>& st)");
     return os;
 }
 
diff --git a/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H b/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H
index 74853b3f5a63c82a902c1a03a9bd23142e07b80e..a0cff64517e7a484dffa5e9c179a53897e145fce 100644
--- a/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H
+++ b/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGas.H
@@ -36,7 +36,6 @@ SourceFiles
 #ifndef perfectGas_H
 #define perfectGas_H
 
-#include "specie.H"
 #include "autoPtr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,13 +43,54 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+
+template<class Specie> class perfectGas;
+
+template<class Specie>
+inline perfectGas<Specie> operator+
+(
+    const perfectGas<Specie>&,
+    const perfectGas<Specie>&
+);
+
+template<class Specie>
+inline perfectGas<Specie> operator-
+(
+    const perfectGas<Specie>&,
+    const perfectGas<Specie>&
+);
+
+template<class Specie>
+inline perfectGas<Specie> operator*
+(
+    const scalar,
+    const perfectGas<Specie>&
+);
+
+template<class Specie>
+inline perfectGas<Specie> operator==
+(
+    const perfectGas<Specie>&,
+    const perfectGas<Specie>&
+);
+
+template<class Specie>
+Ostream& operator<<
+(
+    Ostream&,
+    const perfectGas<Specie>&
+);
+
+
 /*---------------------------------------------------------------------------*\
                            Class perfectGas Declaration
 \*---------------------------------------------------------------------------*/
 
+template<class Specie>
 class perfectGas
 :
-    public specie
+    public Specie
 {
 
 public:
@@ -58,7 +98,7 @@ public:
     // Constructors
 
         //- Construct from components
-        inline perfectGas(const specie& sp);
+        inline perfectGas(const Specie& sp);
 
         //- Construct from Istream
         perfectGas(Istream&);
@@ -118,25 +158,25 @@ public:
 
     // Friend operators
 
-        inline friend perfectGas operator+
+        inline friend perfectGas operator+ <Specie>
         (
             const perfectGas&,
             const perfectGas&
         );
 
-        inline friend perfectGas operator-
+        inline friend perfectGas operator- <Specie>
         (
             const perfectGas&,
             const perfectGas&
         );
 
-        inline friend perfectGas operator*
+        inline friend perfectGas operator* <Specie>
         (
             const scalar s,
             const perfectGas&
         );
 
-        inline friend perfectGas operator==
+        inline friend perfectGas operator== <Specie>
         (
             const perfectGas&,
             const perfectGas&
@@ -145,7 +185,11 @@ public:
 
     // Ostream Operator
 
-        friend Ostream& operator<<(Ostream&, const perfectGas&);
+        friend Ostream& operator<< <Specie>
+        (
+            Ostream&,
+            const perfectGas&
+        );
 };
 
 
@@ -157,6 +201,10 @@ public:
 
 #include "perfectGasI.H"
 
+#ifdef NoRepository
+#   include "perfectGas.C"
+#endif
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGasI.H b/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGasI.H
index b43bd6fe1a24e8b2b2ede56af1c8970fa095ca8a..7d5c6832ba8ca98ef77123a74b84895ebf9d1901 100644
--- a/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGasI.H
+++ b/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGasI.H
@@ -27,62 +27,77 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-inline Foam::perfectGas::perfectGas(const specie& sp)
+template<class Specie>
+inline Foam::perfectGas<Specie>::perfectGas(const Specie& sp)
 :
-    specie(sp)
+    Specie(sp)
 {}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-inline Foam::perfectGas::perfectGas(const word& name, const perfectGas& pg)
+template<class Specie>
+inline Foam::perfectGas<Specie>::perfectGas
+(
+    const word& name,
+    const perfectGas<Specie>& pg
+)
 :
-    specie(name, pg)
+    Specie(name, pg)
 {}
 
 
-inline Foam::autoPtr<Foam::perfectGas> Foam::perfectGas::clone() const
+template<class Specie>
+inline Foam::autoPtr<Foam::perfectGas<Specie> >
+Foam::perfectGas<Specie>::clone() const
 {
-    return autoPtr<perfectGas>(new perfectGas(*this));
+    return autoPtr<perfectGas<Specie> >(new perfectGas<Specie>(*this));
 }
 
 
-inline Foam::autoPtr<Foam::perfectGas> Foam::perfectGas::New(Istream& is)
+template<class Specie>
+inline Foam::autoPtr<Foam::perfectGas<Specie> >
+Foam::perfectGas<Specie>::New(Istream& is)
 {
-    return autoPtr<perfectGas>(new perfectGas(is));
+    return autoPtr<perfectGas<Specie> >(new perfectGas<Specie>(is));
 }
 
 
-inline Foam::autoPtr<Foam::perfectGas> Foam::perfectGas::New
+template<class Specie>
+inline Foam::autoPtr<Foam::perfectGas<Specie> > Foam::perfectGas<Specie>::New
 (
     const dictionary& dict
 )
 {
-    return autoPtr<perfectGas>(new perfectGas(dict));
+    return autoPtr<perfectGas<Specie> >(new perfectGas<Specie>(dict));
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::scalar Foam::perfectGas::rho(scalar p, scalar T) const
+template<class Specie>
+inline Foam::scalar Foam::perfectGas<Specie>::rho(scalar p, scalar T) const
 {
-    return p/(R()*T);
+    return p/(this->R()*T);
 }
 
 
-inline Foam::scalar Foam::perfectGas::psi(scalar, scalar T) const
+template<class Specie>
+inline Foam::scalar Foam::perfectGas<Specie>::psi(scalar, scalar T) const
 {
-    return 1.0/(R()*T);
+    return 1.0/(this->R()*T);
 }
 
 
-inline Foam::scalar Foam::perfectGas::Z(scalar, scalar) const
+template<class Specie>
+inline Foam::scalar Foam::perfectGas<Specie>::Z(scalar, scalar) const
 {
     return 1.0;
 }
 
 
-inline Foam::scalar Foam::perfectGas::cpMcv(scalar, scalar) const
+template<class Specie>
+inline Foam::scalar Foam::perfectGas<Specie>::cpMcv(scalar, scalar) const
 {
     return this->RR;
 }
@@ -90,68 +105,75 @@ inline Foam::scalar Foam::perfectGas::cpMcv(scalar, scalar) const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-inline void Foam::perfectGas::operator+=(const perfectGas& pg)
+template<class Specie>
+inline void Foam::perfectGas<Specie>::operator+=(const perfectGas<Specie>& pg)
 {
-    specie::operator+=(pg);
+    Specie::operator+=(pg);
 }
 
 
-inline void Foam::perfectGas::operator-=(const perfectGas& pg)
+template<class Specie>
+inline void Foam::perfectGas<Specie>::operator-=(const perfectGas<Specie>& pg)
 {
-    specie::operator-=(pg);
+    Specie::operator-=(pg);
 }
 
 
-inline void Foam::perfectGas::operator*=(const scalar s)
+template<class Specie>
+inline void Foam::perfectGas<Specie>::operator*=(const scalar s)
 {
-    specie::operator*=(s);
+    Specie::operator*=(s);
 }
 
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-inline Foam::perfectGas Foam::operator+
+template<class Specie>
+inline Foam::perfectGas<Specie> Foam::operator+
 (
-    const perfectGas& pg1,
-    const perfectGas& pg2
+    const perfectGas<Specie>& pg1,
+    const perfectGas<Specie>& pg2
 )
 {
-    return perfectGas
+    return perfectGas<Specie>
     (
-        static_cast<const specie&>(pg1)
-      + static_cast<const specie&>(pg2)
+        static_cast<const Specie&>(pg1)
+      + static_cast<const Specie&>(pg2)
     );
 }
 
 
-inline Foam::perfectGas Foam::operator-
+template<class Specie>
+inline Foam::perfectGas<Specie> Foam::operator-
 (
-    const perfectGas& pg1,
-    const perfectGas& pg2
+    const perfectGas<Specie>& pg1,
+    const perfectGas<Specie>& pg2
 )
 {
-    return perfectGas
+    return perfectGas<Specie>
     (
-        static_cast<const specie&>(pg1)
-      - static_cast<const specie&>(pg2)
+        static_cast<const Specie&>(pg1)
+      - static_cast<const Specie&>(pg2)
     );
 }
 
 
-inline Foam::perfectGas Foam::operator*
+template<class Specie>
+inline Foam::perfectGas<Specie> Foam::operator*
 (
     const scalar s,
-    const perfectGas& pg
+    const perfectGas<Specie>& pg
 )
 {
-    return perfectGas(s*static_cast<const specie&>(pg));
+    return perfectGas<Specie>(s*static_cast<const Specie&>(pg));
 }
 
 
-inline Foam::perfectGas Foam::operator==
+template<class Specie>
+inline Foam::perfectGas<Specie> Foam::operator==
 (
-    const perfectGas& pg1,
-    const perfectGas& pg2
+    const perfectGas<Specie>& pg1,
+    const perfectGas<Specie>& pg2
 )
 {
     return pg2 - pg1;
diff --git a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C
index 97e825ffb3e66ce1ba4940b46198eb217feed841..895ec1859f96f8ab57b9811825347cd278ddcc28 100644
--- a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C
+++ b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.C
@@ -28,27 +28,30 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::rhoConst::rhoConst(Istream& is)
+template<class Specie>
+Foam::rhoConst<Specie>::rhoConst(Istream& is)
 :
-    specie(is),
+    Specie(is),
     rho_(readScalar(is))
 {
-    is.check("rhoConst::rhoConst(Istream& is)");
+    is.check("rhoConst<Specie>::rhoConst(Istream& is)");
 }
 
 
-Foam::rhoConst::rhoConst(const dictionary& dict)
+template<class Specie>
+Foam::rhoConst<Specie>::rhoConst(const dictionary& dict)
 :
-    specie(dict),
+    Specie(dict),
     rho_(readScalar(dict.subDict("equationOfState").lookup("rho")))
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::rhoConst::write(Ostream& os) const
+template<class Specie>
+void Foam::rhoConst<Specie>::write(Ostream& os) const
 {
-    specie::write(os);
+    Specie::write(os);
 
     dictionary dict("equationOfState");
     dict.add("rho", rho_);
@@ -59,12 +62,13 @@ void Foam::rhoConst::write(Ostream& os) const
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
-Foam::Ostream& Foam::operator<<(Ostream& os, const rhoConst& ico)
+template<class Specie>
+Foam::Ostream& Foam::operator<<(Ostream& os, const rhoConst<Specie>& ico)
 {
-    os  << static_cast<const specie&>(ico)
+    os  << static_cast<const Specie&>(ico)
         << token::SPACE << ico.rho_;
 
-    os.check("Ostream& operator<<(Ostream& os, const rhoConst& ico)");
+    os.check("Ostream& operator<<(Ostream& os, const rhoConst<Specie>& ico)");
     return os;
 }
 
diff --git a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H
index 49a776c63ef0cd76d2ba74589619e9ffe4b4035b..af4af4adbd8e767fe0fd206048bead4272b0075b 100644
--- a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H
+++ b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConst.H
@@ -36,7 +36,6 @@ SourceFiles
 #ifndef rhoConst_H
 #define rhoConst_H
 
-#include "specie.H"
 #include "autoPtr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,13 +43,54 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+
+template<class Specie> class rhoConst;
+
+template<class Specie>
+inline rhoConst<Specie> operator+
+(
+    const rhoConst<Specie>&,
+    const rhoConst<Specie>&
+);
+
+template<class Specie>
+inline rhoConst<Specie> operator-
+(
+    const rhoConst<Specie>&,
+    const rhoConst<Specie>&
+);
+
+template<class Specie>
+inline rhoConst<Specie> operator*
+(
+    const scalar,
+    const rhoConst<Specie>&
+);
+
+template<class Specie>
+inline rhoConst<Specie> operator==
+(
+    const rhoConst<Specie>&,
+    const rhoConst<Specie>&
+);
+
+template<class Specie>
+Ostream& operator<<
+(
+    Ostream&,
+    const rhoConst<Specie>&
+);
+
+
 /*---------------------------------------------------------------------------*\
                         Class rhoConst Declaration
 \*---------------------------------------------------------------------------*/
 
+template<class Specie>
 class rhoConst
 :
-    public specie
+    public Specie
 {
     // Private data
 
@@ -63,7 +103,7 @@ public:
     // Constructors
 
         //- Construct from components
-        inline rhoConst(const specie& sp, const scalar rho);
+        inline rhoConst(const Specie& sp, const scalar rho);
 
         //- Construct from Istream
         rhoConst(Istream&);
@@ -120,25 +160,25 @@ public:
 
     // Friend operators
 
-        inline friend rhoConst operator+
+        inline friend rhoConst operator+ <Specie>
         (
             const rhoConst&,
             const rhoConst&
         );
 
-        inline friend rhoConst operator-
+        inline friend rhoConst operator- <Specie>
         (
             const rhoConst&,
             const rhoConst&
         );
 
-        inline friend rhoConst operator*
+        inline friend rhoConst operator* <Specie>
         (
             const scalar s,
             const rhoConst&
         );
 
-        inline friend rhoConst operator==
+        inline friend rhoConst operator== <Specie>
         (
             const rhoConst&,
             const rhoConst&
@@ -147,7 +187,11 @@ public:
 
     // Ostream Operator
 
-        friend Ostream& operator<<(Ostream&, const rhoConst&);
+        friend Ostream& operator<< <Specie>
+        (
+            Ostream&,
+            const rhoConst&
+        );
 };
 
 
@@ -159,6 +203,10 @@ public:
 
 #include "rhoConstI.H"
 
+#ifdef NoRepository
+#   include "rhoConst.C"
+#endif
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H
index 79faf014d4eec85ca4887a570bc18fb6ae68cac9..a3898a4d69ffaac8a741ec10d6d6a8af60269e3f 100644
--- a/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H
+++ b/src/thermophysicalModels/specie/equationOfState/rhoConst/rhoConstI.H
@@ -27,65 +27,73 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-inline Foam::rhoConst::rhoConst
+template<class Specie>
+inline Foam::rhoConst<Specie>::rhoConst
 (
-    const specie& sp,
+    const Specie& sp,
     const scalar rho
 )
 :
-    specie(sp),
+    Specie(sp),
     rho_(rho)
 {}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-inline Foam::rhoConst::rhoConst
+template<class Specie>
+inline Foam::rhoConst<Specie>::rhoConst
 (
     const word& name,
-    const rhoConst& ico
+    const rhoConst<Specie>& ico
 )
 :
-    specie(name, ico),
+    Specie(name, ico),
     rho_(ico.rho_)
 {}
 
 
-inline Foam::autoPtr<Foam::rhoConst>
-Foam::rhoConst::clone() const
+template<class Specie>
+inline Foam::autoPtr<Foam::rhoConst<Specie> >
+Foam::rhoConst<Specie>::clone() const
 {
-    return autoPtr<rhoConst>(new rhoConst(*this));
+    return autoPtr<rhoConst<Specie> >(new rhoConst<Specie>(*this));
 }
 
 
-inline Foam::autoPtr<Foam::rhoConst>
-Foam::rhoConst::New(Istream& is)
+template<class Specie>
+inline Foam::autoPtr<Foam::rhoConst<Specie> >
+Foam::rhoConst<Specie>::New(Istream& is)
 {
-    return autoPtr<rhoConst>(new rhoConst(is));
+    return autoPtr<rhoConst<Specie> >(new rhoConst<Specie>(is));
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::scalar Foam::rhoConst::rho(scalar p, scalar T) const
+template<class Specie>
+inline Foam::scalar Foam::rhoConst<Specie>::rho(scalar p, scalar T) const
 {
     return rho_;
 }
 
 
-inline Foam::scalar Foam::rhoConst::psi(scalar, scalar T) const
+template<class Specie>
+inline Foam::scalar Foam::rhoConst<Specie>::psi(scalar, scalar T) const
 {
     return 0.0;
 }
 
 
-inline Foam::scalar Foam::rhoConst::Z(scalar, scalar) const
+template<class Specie>
+inline Foam::scalar Foam::rhoConst<Specie>::Z(scalar, scalar) const
 {
     return 0.0;
 }
 
 
-inline Foam::scalar Foam::rhoConst::cpMcv(scalar, scalar) const
+template<class Specie>
+inline Foam::scalar Foam::rhoConst<Specie>::cpMcv(scalar, scalar) const
 {
     return 0.0;
 }
@@ -93,11 +101,12 @@ inline Foam::scalar Foam::rhoConst::cpMcv(scalar, scalar) const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-inline void Foam::rhoConst::operator+=(const rhoConst& ico)
+template<class Specie>
+inline void Foam::rhoConst<Specie>::operator+=(const rhoConst<Specie>& ico)
 {
     scalar molr1 = this->nMoles();
 
-    specie::operator+=(ico);
+    Specie::operator+=(ico);
 
     molr1 /= this->nMoles();
     scalar molr2 = ico.nMoles()/this->nMoles();
@@ -106,11 +115,12 @@ inline void Foam::rhoConst::operator+=(const rhoConst& ico)
 }
 
 
-inline void Foam::rhoConst::operator-=(const rhoConst& ico)
+template<class Specie>
+inline void Foam::rhoConst<Specie>::operator-=(const rhoConst<Specie>& ico)
 {
     scalar molr1 = this->nMoles();
 
-    specie::operator-=(ico);
+    Specie::operator-=(ico);
 
     molr1 /= this->nMoles();
     scalar molr2 = ico.nMoles()/this->nMoles();
@@ -119,66 +129,71 @@ inline void Foam::rhoConst::operator-=(const rhoConst& ico)
 }
 
 
-inline void Foam::rhoConst::operator*=(const scalar s)
+template<class Specie>
+inline void Foam::rhoConst<Specie>::operator*=(const scalar s)
 {
-    specie::operator*=(s);
+    Specie::operator*=(s);
 }
 
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-inline Foam::rhoConst Foam::operator+
+template<class Specie>
+inline Foam::rhoConst<Specie> Foam::operator+
 (
-    const rhoConst& ico1,
-    const rhoConst& ico2
+    const rhoConst<Specie>& ico1,
+    const rhoConst<Specie>& ico2
 )
 {
     scalar nMoles = ico1.nMoles() + ico2.nMoles();
     scalar molr1 = ico1.nMoles()/nMoles;
     scalar molr2 = ico2.nMoles()/nMoles;
 
-    return rhoConst
+    return rhoConst<Specie>
     (
-        static_cast<const specie&>(ico1)
-      + static_cast<const specie&>(ico2),
+        static_cast<const Specie&>(ico1)
+      + static_cast<const Specie&>(ico2),
         molr1*ico1.rho_ + molr2*ico2.rho_
     );
 }
 
 
-inline Foam::rhoConst Foam::operator-
+template<class Specie>
+inline Foam::rhoConst<Specie> Foam::operator-
 (
-    const rhoConst& ico1,
-    const rhoConst& ico2
+    const rhoConst<Specie>& ico1,
+    const rhoConst<Specie>& ico2
 )
 {
     scalar nMoles = ico1.nMoles() + ico2.nMoles();
     scalar molr1 = ico1.nMoles()/nMoles;
     scalar molr2 = ico2.nMoles()/nMoles;
 
-    return rhoConst
+    return rhoConst<Specie>
     (
-        static_cast<const specie&>(ico1)
-      - static_cast<const specie&>(ico2),
+        static_cast<const Specie&>(ico1)
+      - static_cast<const Specie&>(ico2),
         molr1*ico1.rho_ - molr2*ico2.rho_
     );
 }
 
 
-inline Foam::rhoConst Foam::operator*
+template<class Specie>
+inline Foam::rhoConst<Specie> Foam::operator*
 (
     const scalar s,
-    const rhoConst& ico
+    const rhoConst<Specie>& ico
 )
 {
-    return rhoConst(s*static_cast<const specie&>(ico), ico.rho_);
+    return rhoConst<Specie>(s*static_cast<const Specie&>(ico), ico.rho_);
 }
 
 
-inline Foam::rhoConst Foam::operator==
+template<class Specie>
+inline Foam::rhoConst<Specie> Foam::operator==
 (
-    const rhoConst& ico1,
-    const rhoConst& ico2
+    const rhoConst<Specie>& ico1,
+    const rhoConst<Specie>& ico2
 )
 {
     return ico2 - ico1;
diff --git a/src/thermophysicalModels/specie/include/thermoPhysicsTypes.H b/src/thermophysicalModels/specie/include/thermoPhysicsTypes.H
index d4cefcc05f5d69dabe63b5a74e2af3fefae9e0d3..a3a276b926a44c307671aaeecfcd2a809fa72601 100644
--- a/src/thermophysicalModels/specie/include/thermoPhysicsTypes.H
+++ b/src/thermophysicalModels/specie/include/thermoPhysicsTypes.H
@@ -32,6 +32,7 @@ Description
 #ifndef thermoPhysicsTypes_H
 #define thermoPhysicsTypes_H
 
+#include "specie.H"
 #include "perfectGas.H"
 #include "incompressiblePerfectGas.H"
 #include "hConstThermo.H"
@@ -56,7 +57,7 @@ namespace Foam
         <
             hConstThermo
             <
-                perfectGas
+                perfectGas<specie>
             >,
             sensibleEnthalpy
         >
@@ -69,7 +70,7 @@ namespace Foam
         <
             janafThermo
             <
-                perfectGas
+                perfectGas<specie>
             >,
             sensibleEnthalpy
         >
@@ -82,7 +83,7 @@ namespace Foam
         <
             hConstThermo
             <
-                incompressiblePerfectGas
+                incompressiblePerfectGas<specie>
             >,
             sensibleEnthalpy
         >
@@ -95,7 +96,7 @@ namespace Foam
         <
             janafThermo
             <
-                incompressiblePerfectGas
+                incompressiblePerfectGas<specie>
             >,
             sensibleEnthalpy
         >
@@ -108,7 +109,7 @@ namespace Foam
         <
             hPolynomialThermo
             <
-                icoPolynomial<8>,
+                icoPolynomial<specie, 8>,
                 8
             >,
             sensibleEnthalpy