diff --git a/applications/solvers/combustion/PDRFoam/PDRFoam.C b/applications/solvers/combustion/PDRFoam/PDRFoam.C
index 89a7b7f956d2d4b5c03e06f28bb4b0ef4dd33cbc..daeed4ebf086e874af5251c7101bcba8e20de81b 100644
--- a/applications/solvers/combustion/PDRFoam/PDRFoam.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoam.C
@@ -36,7 +36,7 @@ Description
     to be appropriate by comparison with the results from the
     spectral model.
 
-    Strain effects are encorporated directly into the Xi equation
+    Strain effects are incorporated directly into the Xi equation
     but not in the algebraic approximation.  Further work need to be
     done on this issue, particularly regarding the enhanced removal rate
     caused by flame compression.  Analysis using results of the spectral
@@ -70,53 +70,53 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
+    #include "setRootCase.H"
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readCombustionProperties.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readCombustionProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
-scalar StCoNum = 0.0;
+    scalar StCoNum = 0.0;
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
-#       include "readPISOControls.H"
-#       include "CourantNo.H"
-#       include "setDeltaT.H"
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "\n\nTime = " << runTime.timeName() << endl;
 
-#       include "rhoEqn.H"
-#       include "UEqn.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=1; corr<=nCorr; corr++)
         {
-#           include "bEqn.H"
-#           include "ftEqn.H"
-#           include "huEqn.H"
-#           include "hEqn.H"
+            #include "bEqn.H"
+            #include "ftEqn.H"
+            #include "huEqn.H"
+            #include "hEqn.H"
 
             if (!ign.ignited())
             {
                 hu == h;
             }
 
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();
diff --git a/applications/solvers/combustion/PDRFoam/StCourantNo.H b/applications/solvers/combustion/PDRFoam/StCourantNo.H
index e529409475414941805f21949257eaac574a0f99..18d6e86a941045877155f7d01b3c844984c76965 100644
--- a/applications/solvers/combustion/PDRFoam/StCourantNo.H
+++ b/applications/solvers/combustion/PDRFoam/StCourantNo.H
@@ -31,23 +31,25 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-scalar meanStCoNum = 0.0;
+    scalar meanStCoNum = 0.0;
 
-if (mesh.nInternalFaces())
-{
-    surfaceScalarField SfUfbyDelta = 
-        mesh.surfaceInterpolation::deltaCoeffs()
-       *mag(phiSt/fvc::interpolate(rho));
+    if (mesh.nInternalFaces())
+    {
+        surfaceScalarField SfUfbyDelta =
+            mesh.surfaceInterpolation::deltaCoeffs()
+        *mag(phiSt/fvc::interpolate(rho));
 
-    StCoNum = max(SfUfbyDelta/mesh.magSf())
-        .value()*runTime.deltaT().value();
+        StCoNum =
+            max(SfUfbyDelta/mesh.magSf()).value()
+           *runTime.deltaT().value();
 
-    meanStCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
-        .value()*runTime.deltaT().value();
-}
+        meanStCoNum =
+            (sum(SfUfbyDelta)/sum(mesh.magSf())).value()
+        *runTime.deltaT().value();
+    }
 
-Info<< "St courant Number mean: " << meanStCoNum
-    << " max: " << StCoNum << endl;
+    Info<< "St courant Number mean: " << meanStCoNum
+        << " max: " << StCoNum << endl;
 }
 
 // ************************************************************************* //
diff --git a/applications/solvers/combustion/PDRFoam/bEqn.H b/applications/solvers/combustion/PDRFoam/bEqn.H
index 3fd99b65eaeaed0e25ce2c7e90f19fde9ff7ccf1..cb4493154030fc42af097a77dfd9258930ac62ec 100644
--- a/applications/solvers/combustion/PDRFoam/bEqn.H
+++ b/applications/solvers/combustion/PDRFoam/bEqn.H
@@ -1,7 +1,7 @@
 tmp<fv::convectionScheme<scalar> > mvConvection
 (
     fv::convectionScheme<scalar>::New
-    (    
+    (
         mesh,
         fields,
         phi,
@@ -25,7 +25,7 @@ if (ign.ignited())
 
     // Unburnt gas density
     // ~~~~~~~~~~~~~~~~~~~
-    volScalarField rhou = thermo->rhou();
+    volScalarField rhou = thermo.rhou();
 
     // Calculate flame normal etc.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/applications/solvers/combustion/PDRFoam/createFields.H b/applications/solvers/combustion/PDRFoam/createFields.H
index ba07a6dd4659b3c5e127efa8b96abd634611d42a..74a18ab6f5affe39d56ff90d41f084cd6e09e4f3 100644
--- a/applications/solvers/combustion/PDRFoam/createFields.H
+++ b/applications/solvers/combustion/PDRFoam/createFields.H
@@ -1,10 +1,11 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<hhuCombustionThermo> thermo
+    autoPtr<hhuCombustionThermo> pThermo
     (
         hhuCombustionThermo::New(mesh)
     );
-    combustionMixture& composition = thermo->composition();
+    hhuCombustionThermo& thermo = pThermo();
+    basicMultiComponentMixture& composition = thermo.composition();
 
     volScalarField rho
     (
@@ -16,13 +17,13 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
-    const volScalarField& psi = thermo->psi();
-    volScalarField& h = thermo->h();
-    volScalarField& hu = thermo->hu();
+    volScalarField& p = thermo.p();
+    const volScalarField& psi = thermo.psi();
+    volScalarField& h = thermo.h();
+    volScalarField& hu = thermo.hu();
 
     volScalarField& b = composition.Y("b");
     Info<< "min(b) = " << min(b).value() << endl;
@@ -54,7 +55,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
diff --git a/applications/solvers/combustion/PDRFoam/hEqn.H b/applications/solvers/combustion/PDRFoam/hEqn.H
index 807792948793710b37dea4320ffcfb34ec933db2..7f5292d01a903f7931aa0015e81cb26c9250253a 100644
--- a/applications/solvers/combustion/PDRFoam/hEqn.H
+++ b/applications/solvers/combustion/PDRFoam/hEqn.H
@@ -8,5 +8,5 @@
         betav*DpDt
     );
 
-    thermo->correct();
+    thermo.correct();
 }
diff --git a/applications/solvers/combustion/PDRFoam/huEqn.H b/applications/solvers/combustion/PDRFoam/huEqn.H
index 462f271f4b75dd60c6efc390b87b9ffc959e32e6..3467bc6b751177800f70b927250e2b61926ff04b 100644
--- a/applications/solvers/combustion/PDRFoam/huEqn.H
+++ b/applications/solvers/combustion/PDRFoam/huEqn.H
@@ -13,6 +13,6 @@ if (ign.ignited())
     //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hu)
 
      ==
-        betav*DpDt*rho/thermo->rhou()
+        betav*DpDt*rho/thermo.rhou()
     );
 }
diff --git a/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H b/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H
index 1ff6f981bbae09ec27284867cec07f78c0292a97..1931dbbb3554e8b9747fd107b799c6e9f0e4ce4c 100644
--- a/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H
+++ b/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H
@@ -196,8 +196,7 @@ public:
 
 
     // Destructor
-
-        ~SCOPE();
+    ~SCOPE();
 
 
     // Member functions
diff --git a/applications/solvers/combustion/PDRFoam/pEqn.H b/applications/solvers/combustion/PDRFoam/pEqn.H
index 312de8301a25a9f843956a817d4f1beb5c27c475..524c8eac74935856805699ace15ed3bc019fb29b 100644
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/applications/solvers/combustion/PDRFoam/pEqn.H
@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = invA & UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -34,7 +34,7 @@ if (transonic)
 }
 else
 {
-    phi = 
+    phi =
         fvc::interpolate(rho)*
         (
             (fvc::interpolate(U) & mesh.Sf())
diff --git a/applications/solvers/combustion/XiFoam/XiFoam.C b/applications/solvers/combustion/XiFoam/XiFoam.C
index d18ef0d2f4b91086fcd6bcdd93ad12a5942eb523..cb7f7d2f959a4542ea5462e3943b05c52bbe37bc 100644
--- a/applications/solvers/combustion/XiFoam/XiFoam.C
+++ b/applications/solvers/combustion/XiFoam/XiFoam.C
@@ -109,7 +109,7 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 
diff --git a/applications/solvers/combustion/XiFoam/bEqn.H b/applications/solvers/combustion/XiFoam/bEqn.H
index 739d6987e833f64c06880d7b31d85351a1c77933..33ef24bfe655b2d587b7e403fbd8ee1f5d25b736 100644
--- a/applications/solvers/combustion/XiFoam/bEqn.H
+++ b/applications/solvers/combustion/XiFoam/bEqn.H
@@ -6,7 +6,7 @@ if (ign.ignited())
 
     // Unburnt gas density
     // ~~~~~~~~~~~~~~~~~~~
-    volScalarField rhou = thermo->rhou();
+    volScalarField rhou = thermo.rhou();
 
     // Calculate flame normal etc.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -76,7 +76,7 @@ if (ign.ignited())
 
     volScalarField epsilon = pow(uPrimeCoef, 3)*turbulence->epsilon();
 
-    volScalarField tauEta = sqrt(thermo->muu()/(rhou*epsilon));
+    volScalarField tauEta = sqrt(thermo.muu()/(rhou*epsilon));
 
     volScalarField Reta = up/
     (
@@ -180,7 +180,7 @@ if (ign.ignited())
         // with a linear correction function to give a plausible profile for Xi
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        volScalarField XiEqStar = 
+        volScalarField XiEqStar =
             scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta;
 
         volScalarField XiEq =
diff --git a/applications/solvers/combustion/XiFoam/createFields.H b/applications/solvers/combustion/XiFoam/createFields.H
index 812188640279377c734494a62d05365bf982f28c..ef16bd615c9bef3cfd86788e8fbd5a8a2bb879f3 100644
--- a/applications/solvers/combustion/XiFoam/createFields.H
+++ b/applications/solvers/combustion/XiFoam/createFields.H
@@ -1,10 +1,11 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<hhuCombustionThermo> thermo
+    autoPtr<hhuCombustionThermo> pThermo
     (
         hhuCombustionThermo::New(mesh)
     );
-    combustionMixture& composition = thermo->composition();
+    hhuCombustionThermo& thermo = pThermo();
+    basicMultiComponentMixture& composition = thermo.composition();
 
     volScalarField rho
     (
@@ -16,18 +17,18 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
-    const volScalarField& psi = thermo->psi();
-    volScalarField& h = thermo->h();
-    volScalarField& hu = thermo->hu();
+    volScalarField& p = thermo.p();
+    const volScalarField& psi = thermo.psi();
+    volScalarField& h = thermo.h();
+    volScalarField& hu = thermo.hu();
 
     volScalarField& b = composition.Y("b");
     Info<< "min(b) = " << min(b).value() << endl;
 
-    const volScalarField& T = thermo->T();
+    const volScalarField& T = thermo.T();
 
 
     Info<< "\nReading field U\n" << endl;
@@ -55,7 +56,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
diff --git a/applications/solvers/combustion/XiFoam/hEqn.H b/applications/solvers/combustion/XiFoam/hEqn.H
index 2a5204aa8d21bbc24d91eb0f6d9c29d1f43149a0..ebce30e24e50dee92824de875abfa138f0abf3b1 100644
--- a/applications/solvers/combustion/XiFoam/hEqn.H
+++ b/applications/solvers/combustion/XiFoam/hEqn.H
@@ -8,5 +8,5 @@
         DpDt
     );
 
-    thermo->correct();
+    thermo.correct();
 }
diff --git a/applications/solvers/combustion/XiFoam/huEqn.H b/applications/solvers/combustion/XiFoam/huEqn.H
index 2f001b003265fca1aeeb2be7dd223546bf845bb3..0b4068344bc8b74a46691ba7ed205c6cd4efc9ba 100644
--- a/applications/solvers/combustion/XiFoam/huEqn.H
+++ b/applications/solvers/combustion/XiFoam/huEqn.H
@@ -13,6 +13,6 @@ if (ign.ignited())
     //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hu)
 
      ==
-        DpDt*rho/thermo->rhou()
+        DpDt*rho/thermo.rhou()
     );
 }
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index 05db89627dabf5400a98c2e6a37784dcb9cc53cd..9443f909a356699185bfbc8497cf46e26fd1ed3b 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -34,9 +34,9 @@ if (transonic)
 }
 else
 {
-    phi = 
-        fvc::interpolate(rho)*
-        (
+    phi =
+        fvc::interpolate(rho)
+       *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
diff --git a/applications/solvers/combustion/coalChemistryFoam/createFields.H b/applications/solvers/combustion/coalChemistryFoam/createFields.H
index b8f04a22d8c25b89a7437f6977a838a18ee7d39e..faf324fee108df509df0c88a54515afdd1468108 100644
--- a/applications/solvers/combustion/coalChemistryFoam/createFields.H
+++ b/applications/solvers/combustion/coalChemistryFoam/createFields.H
@@ -84,7 +84,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
     DimensionedField<scalar, volMesh> kappa
     (
diff --git a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
index ddb4c3553009eb6acb7284c01e7c1f6517dfb89e..67def3fb5901fc279f6e51754b2c7070d658b185 100644
--- a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
+++ b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
@@ -33,7 +33,7 @@ Description
 #include "fvCFD.H"
 #include "engineTime.H"
 #include "engineMesh.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "OFstream.H"
 
@@ -41,27 +41,27 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
+    #include "setRootCase.H"
 
-#   include "createEngineTime.H"
-#   include "createEngineMesh.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-#   include "readEngineTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-#   include "startSummary.H"
+    #include "createEngineTime.H"
+    #include "createEngineMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readEngineTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+    #include "startSummary.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readEngineTimeControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readEngineTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -70,22 +70,22 @@ int main(int argc, char *argv[])
 
         mesh.move();
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=1; corr<=nCorr; corr++)
         {
-#           include "hEqn.H"
-#           include "pEqn.H"
+            #include "hEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();
 
         runTime.write();
 
-#       include "logSummary.H"
+        #include "logSummary.H"
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
diff --git a/applications/solvers/combustion/coldEngineFoam/createFields.H b/applications/solvers/combustion/coldEngineFoam/createFields.H
index 7a369df4721e96dd98cf3d86d3d581d86ba59ef2..6bc3139ce6f070445cc5e55ba3b3df839522ab50 100644
--- a/applications/solvers/combustion/coldEngineFoam/createFields.H
+++ b/applications/solvers/combustion/coldEngineFoam/createFields.H
@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -15,13 +16,13 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
-    const volScalarField& psi = thermo->psi();
-    volScalarField& h = thermo->h();
-    const volScalarField& T = thermo->T();
+    volScalarField& p = thermo.p();
+    const volScalarField& psi = thermo.psi();
+    volScalarField& h = thermo.h();
+    const volScalarField& T = thermo.T();
 
 
     Info<< "\nReading field U\n" << endl;
@@ -38,7 +39,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
 
     Info<< "Creating turbulence model\n" << endl;
@@ -49,7 +50,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
diff --git a/applications/solvers/combustion/coldEngineFoam/hEqn.H b/applications/solvers/combustion/coldEngineFoam/hEqn.H
index f72ef0c89cb973b0150b3ec889b51c5be403b872..ae60d3316ec77061804e629360ed13f6cd891f68 100644
--- a/applications/solvers/combustion/coldEngineFoam/hEqn.H
+++ b/applications/solvers/combustion/coldEngineFoam/hEqn.H
@@ -8,5 +8,5 @@
         DpDt
     );
 
-    thermo->correct();
+    thermo.correct();
 }
diff --git a/applications/solvers/combustion/dieselEngineFoam/createFields.H b/applications/solvers/combustion/dieselEngineFoam/createFields.H
index cf0532b0f3a9079454779e4cb04a1054277129b3..9d9229cc3c1dbed27219d4e53ff49af24f78767d 100644
--- a/applications/solvers/combustion/dieselEngineFoam/createFields.H
+++ b/applications/solvers/combustion/dieselEngineFoam/createFields.H
@@ -1,13 +1,17 @@
 Info<< nl << "Reading thermophysicalProperties" << endl;
-autoPtr<hCombustionThermo> thermo
+
+autoPtr<psiChemistryModel> pChemistry
 (
-    hCombustionThermo::New(mesh)
+    psiChemistryModel::New(mesh)
 );
+psiChemistryModel& chemistry = pChemistry();
+
+hCombustionThermo& thermo = chemistry.thermo();
 
-combustionMixture& composition = thermo->composition();
+basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 
-word inertSpecie(thermo->lookup("inertSpecie"));
+word inertSpecie(thermo.lookup("inertSpecie"));
 
 volScalarField rho
 (
@@ -17,7 +21,7 @@ volScalarField rho
         runTime.timeName(),
         mesh
     ),
-    thermo->rho()
+    thermo.rho()
 );
 
 Info<< "Reading field U\n" << endl;
@@ -35,10 +39,10 @@ volVectorField U
 );
 
 
-volScalarField& p = thermo->p();
-const volScalarField& psi = thermo->psi();
-const volScalarField& T = thermo->T();
-volScalarField& h = thermo->h();
+volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
+const volScalarField& T = thermo.T();
+volScalarField& h = thermo.h();
 
 
 #include "compressibleCreatePhi.H"
@@ -65,7 +69,7 @@ autoPtr<compressible::turbulenceModel> turbulence
         rho,
         U,
         phi,
-        thermo()
+        thermo
     )
 );
 
@@ -73,31 +77,11 @@ Info<< "Creating field DpDt\n" << endl;
 volScalarField DpDt =
     fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
-Info << "Constructing chemical mechanism" << endl;
-chemistryModel chemistry
-(
-    thermo(),
-    rho
-);
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-for(label i=0; i<Y.size(); i++)
+forAll (Y, i)
 {
     fields.add(Y[i]);
 }
 fields.add(h);
-
-volScalarField dQ
-(
-    IOobject
-    (
-        "dQ",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("zero", dimensionSet(1,-3,-1,0,0,0,0), 0.0)
-);
diff --git a/applications/solvers/combustion/dieselEngineFoam/createSpray.H b/applications/solvers/combustion/dieselEngineFoam/createSpray.H
index ac473957e5d1e86395509f79c2175c4f8039bc4d..4bc1a32f11becdd7e88f040685a97ae98be7e4e3 100644
--- a/applications/solvers/combustion/dieselEngineFoam/createSpray.H
+++ b/applications/solvers/combustion/dieselEngineFoam/createSpray.H
@@ -1,14 +1,15 @@
 Info << "Constructing Spray" << endl;
 
-PtrList<specieProperties> gasProperties(Y.size());
+PtrList<gasThermoPhysics> gasProperties(Y.size());
 forAll(gasProperties, i)
 {
     gasProperties.set
     (
         i,
-        new specieProperties
+        new gasThermoPhysics
         (
-            dynamic_cast<const reactingMixture&>(thermo()).speciesData()[i]
+            dynamic_cast<const reactingMixture<gasThermoPhysics>&>
+                (thermo).speciesData()[i]
         )
     );
 }
diff --git a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
index 927afe9bc28e486ad50704f7b9ad2daecbe88ba8..ba910bcf56e2bd70d03251deb46429622fb3a0af 100644
--- a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
+++ b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
@@ -36,40 +36,41 @@ Description
 #include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "spray.H"
-#include "chemistryModel.H"
+#include "psiChemistryModel.H"
 #include "chemistrySolver.H"
 #include "multivariateScheme.H"
 #include "Switch.H"
 #include "OFstream.H"
 #include "volPointInterpolation.H"
+#include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
-#   include "createEngineTime.H"
-#   include "createEngineMesh.H"
-#   include "createFields.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readCombustionProperties.H"
-#   include "createSpray.H"
-#   include "initContinuityErrs.H"
-#   include "readEngineTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-#   include "startSummary.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    #include "setRootCase.H"
+    #include "createEngineTime.H"
+    #include "createEngineMesh.H"
+    #include "createFields.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readCombustionProperties.H"
+    #include "createSpray.H"
+    #include "initContinuityErrs.H"
+    #include "readEngineTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+    #include "startSummary.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readEngineTimeControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readEngineTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -101,27 +102,27 @@ int main(int argc, char *argv[])
             kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
         }
 
-#       include "rhoEqn.H"
-#       include "UEqn.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
-#           include "YEqn.H"
-#           include "hEqn.H"
+            #include "YEqn.H"
+            #include "hEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-#               include "pEqn.H"
+                #include "pEqn.H"
             }
         }
 
         turbulence->correct();
 
-#       include "logSummary.H"
-#       include "spraySummary.H"
+        #include "logSummary.H"
+        #include "spraySummary.H"
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 
diff --git a/applications/solvers/combustion/dieselEngineFoam/hEqn.H b/applications/solvers/combustion/dieselEngineFoam/hEqn.H
index 347fef1a9c33231900aeecc5703f7e120f26fc4c..0406f7fbd492100a3814d0531ef25487636df4c1 100644
--- a/applications/solvers/combustion/dieselEngineFoam/hEqn.H
+++ b/applications/solvers/combustion/dieselEngineFoam/hEqn.H
@@ -9,32 +9,5 @@
      + dieselSpray.heatTransferSource()
     );
 
-    thermo->correct();
-
-    forAll(dQ, i)
-    {
-        dQ[i] = 0.0;
-    }
-
-    scalarField cp(dQ.size(), 0.0);
-
-    forAll(Y, i)
-    {
-        volScalarField RRi = chemistry.RR(i);
-
-        forAll(h, celli)
-        {
-            scalar Ti = T[celli];
-            cp[celli] += Y[i][celli]*chemistry.specieThermo()[i].Cp(Ti);
-            scalar hi = chemistry.specieThermo()[i].h(Ti);
-            scalar RR = RRi[celli];
-            dQ[celli] -= hi*RR;
-        }
-
-    }
-
-    forAll(dQ, celli)
-    {
-        dQ[celli] /= cp[celli];
-    }
+    thermo.correct();
 }
diff --git a/applications/solvers/combustion/dieselEngineFoam/pEqn.H b/applications/solvers/combustion/dieselEngineFoam/pEqn.H
index 0324a47ce01ecdc7496cdacc24b3a137b01fb703..b68ae732962c2e1f26c46876b69efae206eaaa45 100644
--- a/applications/solvers/combustion/dieselEngineFoam/pEqn.H
+++ b/applications/solvers/combustion/dieselEngineFoam/pEqn.H
@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField A = UEqn.A();
 U = UEqn.H()/A;
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
diff --git a/applications/solvers/combustion/dieselEngineFoam/rhoEqn.H b/applications/solvers/combustion/dieselEngineFoam/rhoEqn.H
index 396bf33f300d9549874fafa6aa92ae9b28d3a895..dd396486d35d57a7b96baeb69c944e461ab7c11f 100644
--- a/applications/solvers/combustion/dieselEngineFoam/rhoEqn.H
+++ b/applications/solvers/combustion/dieselEngineFoam/rhoEqn.H
@@ -44,7 +44,7 @@ volScalarField Sevap
     dimensionedScalar("zero", dimensionSet(1, -3, -1, 0, 0), 0.0)
 );
 
-for(label i=0; i<Y.size(); i++)
+for (label i=0; i<Y.size(); i++)
 {
     if (dieselSpray.isLiquidFuel()[i])
     {
diff --git a/applications/solvers/combustion/dieselEngineFoam/spraySummary.H b/applications/solvers/combustion/dieselEngineFoam/spraySummary.H
index 5b251e22902437f4b4bf964ee25c3d9c3169e90e..31337d0a77ee0dea6b686bd8f1c2b551d317329a 100644
--- a/applications/solvers/combustion/dieselEngineFoam/spraySummary.H
+++ b/applications/solvers/combustion/dieselEngineFoam/spraySummary.H
@@ -1,30 +1,30 @@
-        label Nparcels = dieselSpray.size();
-        reduce(Nparcels, sumOp<label>());
+    label Nparcels = dieselSpray.size();
+    reduce(Nparcels, sumOp<label>());
 
-        Info<< "\nNumber of parcels in system.... | "
-            << Nparcels << endl
-            << "Injected liquid mass........... | "
-            << 1e6*dieselSpray.injectedMass(runTime.value()) << " mg" << endl
-            << "Liquid Mass in system.......... | "
-            << 1e6*dieselSpray.liquidMass() << " mg" << endl
-            << "SMD, Dmax...................... | "
-            << dieselSpray.smd()*1e6 << " mu, "
-            << dieselSpray.maxD()*1e6 << " mu"
-            << endl;
+    Info<< "\nNumber of parcels in system.... | "
+        << Nparcels << endl
+        << "Injected liquid mass........... | "
+        << 1e6*dieselSpray.injectedMass(runTime.value()) << " mg" << endl
+        << "Liquid Mass in system.......... | "
+        << 1e6*dieselSpray.liquidMass() << " mg" << endl
+        << "SMD, Dmax...................... | "
+        << dieselSpray.smd()*1e6 << " mu, "
+        << dieselSpray.maxD()*1e6 << " mu"
+        << endl;
 
-        scalar evapMass =
-            dieselSpray.injectedMass(runTime.value())
-          - dieselSpray.liquidMass();
+    scalar evapMass =
+        dieselSpray.injectedMass(runTime.value())
+    - dieselSpray.liquidMass();
 
-        scalar gasMass = fvc::domainIntegrate(rho).value();
+    scalar gasMass = fvc::domainIntegrate(rho).value();
 
-        if (dieselSpray.twoD())
-        {
-            gasMass *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge();
-        }
+    if (dieselSpray.twoD())
+    {
+        gasMass *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge();
+    }
 
-        scalar addedMass = gasMass - gasMass0;
+    scalar addedMass = gasMass - gasMass0;
 
-        Info<< "Added gas mass................. | " << 1e6*addedMass << " mg"
-            << nl << "Evaporation Continuity Error... | "
-            << 1e6*(addedMass - evapMass) << " mg" << endl;
+    Info<< "Added gas mass................. | " << 1e6*addedMass << " mg"
+        << nl << "Evaporation Continuity Error... | "
+        << 1e6*(addedMass - evapMass) << " mg" << endl;
diff --git a/applications/solvers/combustion/dieselFoam/dieselFoam.C b/applications/solvers/combustion/dieselFoam/dieselFoam.C
index 31d034bc5734e6807843059dff8eea02e17474a8..11cd91a488e349464b1b7eb8f9dce31c26e3a919 100644
--- a/applications/solvers/combustion/dieselFoam/dieselFoam.C
+++ b/applications/solvers/combustion/dieselFoam/dieselFoam.C
@@ -34,7 +34,7 @@ Description
 #include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "spray.H"
-#include "chemistryModel.H"
+#include "psiChemistryModel.H"
 #include "chemistrySolver.H"
 
 #include "multivariateScheme.H"
@@ -46,28 +46,27 @@ Description
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "createFields.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readCombustionProperties.H"
-#   include "createSpray.H"
-#   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readCombustionProperties.H"
+    #include "createSpray.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -94,26 +93,26 @@ int main(int argc, char *argv[])
             kappa = (runTime.deltaT() + tc)/(runTime.deltaT()+tc+tk);
         }
 
-#       include "rhoEqn.H"
-#       include "UEqn.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
-#           include "YEqn.H"
-#           include "hEqn.H"
+            #include "YEqn.H"
+            #include "hEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-#               include "pEqn.H"
+                #include "pEqn.H"
             }
         }
 
         turbulence->correct();
 
-#       include "spraySummary.H"
+        #include "spraySummary.H"
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 
diff --git a/applications/solvers/combustion/dieselFoam/pEqn.H b/applications/solvers/combustion/dieselFoam/pEqn.H
index 37e63e8635cfa1bdef7ec2021711add2c9b7b464..d74947fe5305882c8958f9237c17709ce530dbb6 100644
--- a/applications/solvers/combustion/dieselFoam/pEqn.H
+++ b/applications/solvers/combustion/dieselFoam/pEqn.H
@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -36,9 +36,9 @@ if (transonic)
 }
 else
 {
-    phi = 
-        fvc::interpolate(rho)*
-        (
+    phi =
+        fvc::interpolate(rho)
+       *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
diff --git a/applications/solvers/combustion/engineFoam/engineFoam.C b/applications/solvers/combustion/engineFoam/engineFoam.C
index a835e8ea4295fc3d04efc3665c751c3ec504020c..33a8ea6f5d64deeca30d0c2805a81341e35175fa 100644
--- a/applications/solvers/combustion/engineFoam/engineFoam.C
+++ b/applications/solvers/combustion/engineFoam/engineFoam.C
@@ -63,29 +63,29 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
-
-#   include "createEngineTime.H"
-#   include "createEngineMesh.H"
-#   include "readPISOControls.H"
-#   include "readCombustionProperties.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-#   include "readEngineTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-#   include "startSummary.H"
+    #include "setRootCase.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    #include "createEngineTime.H"
+    #include "createEngineMesh.H"
+    #include "readPISOControls.H"
+    #include "readCombustionProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readEngineTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+    #include "startSummary.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readEngineTimeControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readEngineTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -93,31 +93,31 @@ int main(int argc, char *argv[])
 
         mesh.move();
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=1; corr<=nCorr; corr++)
         {
-#           include "ftEqn.H"
-#           include "bEqn.H"
-#           include "huEqn.H"
-#           include "hEqn.H"
+            #include "ftEqn.H"
+            #include "bEqn.H"
+            #include "huEqn.H"
+            #include "hEqn.H"
 
             if (!ign.ignited())
             {
                 hu == h;
             }
 
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();
 
-#       include "logSummary.H"
+        #include "logSummary.H"
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 638d952f1753a9ea38a5dd317c23ea543f7cc4d1..39b4967312055a320025d3ca3ba70b809fafe762 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
@@ -8,8 +8,8 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
-        *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
+        fvc::interpolate(psi)
+       *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
diff --git a/applications/solvers/combustion/engineFoam/readEngineTimeControls.H b/applications/solvers/combustion/engineFoam/readEngineTimeControls.H
index 2aea3445eb2705ed41469074bcae777250a7f378..8d6d5f28b3b50a9d5fd54f33fd5c0a565e187b84 100644
--- a/applications/solvers/combustion/engineFoam/readEngineTimeControls.H
+++ b/applications/solvers/combustion/engineFoam/readEngineTimeControls.H
@@ -1,3 +1,3 @@
-#include "readTimeControls.H"
+    #include "readTimeControls.H"
 
     maxDeltaT = runTime.userTimeToTime(maxDeltaT);
diff --git a/applications/solvers/combustion/reactingFoam/YEqn.H b/applications/solvers/combustion/reactingFoam/YEqn.H
index 873351ac9b4d4c03c6da75e907ab1eea09fdaac6..cda24ec2f72efc49c797a83bbcbb4ba2be6007b5 100644
--- a/applications/solvers/combustion/reactingFoam/YEqn.H
+++ b/applications/solvers/combustion/reactingFoam/YEqn.H
@@ -13,7 +13,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     label inertIndex = -1;
     volScalarField Yt = 0.0*Y[0];
 
-    for(label i=0; i<Y.size(); i++)
+    for (label i=0; i<Y.size(); i++)
     {
         if (Y[i].name() != inertSpecie)
         {
@@ -37,7 +37,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
             inertIndex = i;
         }
     }
-    
+
     Y[inertIndex] = scalar(1) - Yt;
     Y[inertIndex].max(0.0);
 }
diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H
index 510f7d4236631404dea74eecb2886c3600d2cb4c..d92766db5498508a45847a41101c7a54f5960701 100644
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/createFields.H
@@ -1,13 +1,16 @@
 Info<< nl << "Reading thermophysicalProperties" << endl;
-autoPtr<hCombustionThermo> thermo
+autoPtr<psiChemistryModel> pChemistry
 (
-    hCombustionThermo::New(mesh)
+    psiChemistryModel::New(mesh)
 );
+psiChemistryModel& chemistry = pChemistry();
 
-combustionMixture& composition = thermo->composition();
+hCombustionThermo& thermo = chemistry.thermo();
+
+basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 
-word inertSpecie(thermo->lookup("inertSpecie"));
+word inertSpecie(thermo.lookup("inertSpecie"));
 
 volScalarField rho
 (
@@ -17,7 +20,7 @@ volScalarField rho
         runTime.timeName(),
         mesh
     ),
-    thermo->rho()
+    thermo.rho()
 );
 
 Info<< "Reading field U\n" << endl;
@@ -35,9 +38,9 @@ volVectorField U
 );
 
 
-volScalarField& p = thermo->p();
-const volScalarField& psi = thermo->psi();
-volScalarField& h = thermo->h();
+volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
+volScalarField& h = thermo.h();
 
 
 #include "compressibleCreatePhi.H"
@@ -64,7 +67,7 @@ autoPtr<compressible::turbulenceModel> turbulence
         rho,
         U,
         phi,
-        thermo()
+        thermo
     )
 );
 
@@ -72,31 +75,11 @@ Info<< "Creating field DpDt\n" << endl;
 volScalarField DpDt =
     fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
-Info << "Constructing chemical mechanism" << endl;
-chemistryModel chemistry
-(
-    thermo(),
-    rho
-);
-
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-for(label i=0; i<Y.size(); i++)
+forAll (Y, i)
 {
     fields.add(Y[i]);
 }
 fields.add(h);
 
-volScalarField dQ
-(
-    IOobject
-    (
-        "dQ",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("zero", dimensionSet(1, -3, -1, 0, 0, 0, 0), 0.0)
-);
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index e915818bdc14525026340f04a7e5339bcd08e740..2a9e34e577395d117645f4f88dd19c90f9add5c2 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C
@@ -33,7 +33,7 @@ Description
 #include "fvCFD.H"
 #include "hCombustionThermo.H"
 #include "turbulenceModel.H"
-#include "chemistryModel.H"
+#include "psiChemistryModel.H"
 #include "chemistrySolver.H"
 #include "multivariateScheme.H"
 
@@ -41,52 +41,52 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readChemistryProperties.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readChemistryProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "chemistry.H"
-#       include "rhoEqn.H"
-#       include "UEqn.H"
+        #include "chemistry.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
-#           include "YEqn.H"
+            #include "YEqn.H"
 
-#           define Db turbulence->alphaEff()
-#           include "hEqn.H"
+            #define Db turbulence->alphaEff()
+            #include "hEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-#               include "pEqn.H"
+                #include "pEqn.H"
             }
         }
 
         turbulence->correct();
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 
diff --git a/applications/solvers/compressible/rhoCentralFoam/compressibleCourantNo.H b/applications/solvers/compressible/rhoCentralFoam/compressibleCourantNo.H
index ebec931a54fcda48dd8777c8fc77fac5a42400f8..33b4edc8a7d31da894068a3c7016dc38bf21137e 100644
--- a/applications/solvers/compressible/rhoCentralFoam/compressibleCourantNo.H
+++ b/applications/solvers/compressible/rhoCentralFoam/compressibleCourantNo.H
@@ -38,11 +38,11 @@ if (mesh.nInternalFaces())
     surfaceScalarField amaxSfbyDelta =
         mesh.surfaceInterpolation::deltaCoeffs()*amaxSf;
 
-    CoNum = max(amaxSfbyDelta/mesh.magSf())
-        .value()*runTime.deltaT().value();
+    CoNum = max(amaxSfbyDelta/mesh.magSf()).value()*runTime.deltaT().value();
 
-    meanCoNum = (sum(amaxSfbyDelta)/sum(mesh.magSf()))
-        .value()*runTime.deltaT().value();
+    meanCoNum =
+        (sum(amaxSfbyDelta)/sum(mesh.magSf())).value()
+       *runTime.deltaT().value();
 }
 
 Info<< "Mean and max Courant Numbers = "
diff --git a/applications/solvers/compressible/rhoCentralFoam/createFields.H b/applications/solvers/compressible/rhoCentralFoam/createFields.H
index 3b273bb3dda1ed515b346870e932faa9a73c3123..daf3a88435cc8b9bfbae14ab33085c9dcd72d55e 100644
--- a/applications/solvers/compressible/rhoCentralFoam/createFields.H
+++ b/applications/solvers/compressible/rhoCentralFoam/createFields.H
@@ -1,15 +1,16 @@
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<basicThermo> thermo
+autoPtr<basicPsiThermo> pThermo
 (
-    basicThermo::New(mesh)
+    basicPsiThermo::New(mesh)
 );
+basicPsiThermo& thermo = pThermo();
 
-volScalarField& p = thermo->p();
-volScalarField& h = thermo->h();
-const volScalarField& T = thermo->T();
-const volScalarField& psi = thermo->psi();
-const volScalarField& mu = thermo->mu();
+volScalarField& p = thermo.p();
+volScalarField& h = thermo.h();
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();
+const volScalarField& mu = thermo.mu();
 
 bool inviscid(true);
 if (max(mu.internalField()) > 0.0)
@@ -42,7 +43,7 @@ volScalarField rho
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
-    thermo->rho(),
+    thermo.rho(),
     rhoBoundaryTypes
 );
 
diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoBoundaryTypes.H b/applications/solvers/compressible/rhoCentralFoam/rhoBoundaryTypes.H
index c85f49198529d7ae1c3e400a5c1ae93afef47f1a..7118d24fc9247c5ed8f0301d9a7bb837ba278b0b 100644
--- a/applications/solvers/compressible/rhoCentralFoam/rhoBoundaryTypes.H
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoBoundaryTypes.H
@@ -3,10 +3,7 @@ wordList rhoBoundaryTypes = pbf.types();
 
 forAll(rhoBoundaryTypes, patchi)
 {
-    if
-    (
-        rhoBoundaryTypes[patchi] == "waveTransmissive"
-    )
+    if (rhoBoundaryTypes[patchi] == "waveTransmissive")
     {
         rhoBoundaryTypes[patchi] = zeroGradientFvPatchScalarField::typeName;
     }
diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
index 502406d1322f64edbdd77550ffa90b1f040e1ce1..f418f176374b384112da5f8cd760be97d61efc02 100644
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "zeroGradientFvPatchFields.H"
 #include "fixedRhoFvPatchScalarField.H"
 
@@ -40,18 +40,17 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "readThermophysicalProperties.H"
+    #include "readTimeControls.H"
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "createFields.H"
-#   include "readThermophysicalProperties.H"
-#   include "readTimeControls.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#   include "readFluxScheme.H"
+    #include "readFluxScheme.H"
 
     dimensionedScalar v_zero("v_zero",dimVolume/dimTime, 0.0);
 
@@ -91,7 +90,7 @@ int main(int argc, char *argv[])
         surfaceScalarField phiv_pos = U_pos & mesh.Sf();
         surfaceScalarField phiv_neg = U_neg & mesh.Sf();
 
-        volScalarField c = sqrt(thermo->Cp()/thermo->Cv()*rPsi);
+        volScalarField c = sqrt(thermo.Cp()/thermo.Cv()*rPsi);
         surfaceScalarField cSf_pos = fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf();
         surfaceScalarField cSf_neg = fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf();
 
@@ -102,9 +101,9 @@ int main(int argc, char *argv[])
 
         surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
 
-#       include "compressibleCourantNo.H"
-#       include "readTimeControls.H"
-#       include "setDeltaT.H"
+        #include "compressibleCourantNo.H"
+        #include "readTimeControls.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -183,7 +182,7 @@ int main(int argc, char *argv[])
 
         h = (rhoE + p)/rho - 0.5*magSqr(U);
         h.correctBoundaryConditions();
-        thermo->correct();
+        thermo.correct();
         rhoE.boundaryField() =
             rho.boundaryField()*
             (
@@ -193,15 +192,15 @@ int main(int argc, char *argv[])
 
         if (!inviscid)
         {
-            volScalarField k("k", thermo->Cp()*mu/Pr);
+            volScalarField k("k", thermo.Cp()*mu/Pr);
             solve
             (
                 fvm::ddt(rho, h) - fvc::ddt(rho, h)
-              - fvm::laplacian(thermo->alpha(), h)
-              + fvc::laplacian(thermo->alpha(), h)
+              - fvm::laplacian(thermo.alpha(), h)
+              + fvc::laplacian(thermo.alpha(), h)
               - fvc::laplacian(k, T)
             );
-            thermo->correct();
+            thermo.correct();
             rhoE = rho*(h + 0.5*magSqr(U)) - p;
         }
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index e836704cb1fd79e93dd0dffec96d15d0001830e0..1433ea5cc3d7180c4310f7b0710b913025030977 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -1,13 +1,14 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
-    const volScalarField& psi = thermo->psi();
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -19,7 +20,7 @@
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "Reading field U\n" << endl;
@@ -51,7 +52,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
diff --git a/applications/solvers/compressible/rhoPimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPimpleFoam/hEqn.H
index dbc76740f51da87673524ba4bda189998d837864..e66b99442b87bd6e3a4fdeb4025e94f5f88e09c3 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/hEqn.H
@@ -19,5 +19,5 @@
         hEqn.solve();
     }
 
-    thermo->correct();
+    thermo.correct();
 }
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 7918618404a9653761dff2f60d6b67ce73b91bfb..9c4e9658cc703bb0234c79a797fdfee4d2c8b22d 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn().A();
 U = rUA*UEqn().H();
@@ -13,7 +13,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -99,7 +99,7 @@ else
     // Explicitly relax pressure for momentum corrector
     p.relax();
 
-    rho = thermo->rho();
+    rho = thermo.rho();
     rho.relax();
     Info<< "rho max/min : " << max(rho).value()
         << " " << min(rho).value() << endl;
@@ -117,7 +117,7 @@ bound(p, pMin);
 /*
 if (closedVolume)
 {
-    p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
-        /fvc::domainIntegrate(thermo->psi());
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
 }
 */
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index 8b916b308ef081bc5c90d7695a5bc52df53d4fad..19a26875db50c7d5f224371c229db652adc42ca3 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -35,7 +35,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "bound.H"
 
diff --git a/applications/solvers/compressible/rhoPisoFoam/createFields.H b/applications/solvers/compressible/rhoPisoFoam/createFields.H
index 614f01702df40cf576ae7b94f5a82cd1405ded64..fdb706a14ce0651bf3722244ea51dae8c0df2d6b 100644
--- a/applications/solvers/compressible/rhoPisoFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPisoFoam/createFields.H
@@ -1,13 +1,14 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
-    const volScalarField& psi = thermo->psi();
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -19,7 +20,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "\nReading field U\n" << endl;
@@ -47,7 +48,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
diff --git a/applications/solvers/compressible/rhoPisoFoam/hEqn.H b/applications/solvers/compressible/rhoPisoFoam/hEqn.H
index f72ef0c89cb973b0150b3ec889b51c5be403b872..ae60d3316ec77061804e629360ed13f6cd891f68 100644
--- a/applications/solvers/compressible/rhoPisoFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoPisoFoam/hEqn.H
@@ -8,5 +8,5 @@
         DpDt
     );
 
-    thermo->correct();
+    thermo.correct();
 }
diff --git a/applications/solvers/compressible/rhoPisoFoam/pEqn.H b/applications/solvers/compressible/rhoPisoFoam/pEqn.H
index 05db89627dabf5400a98c2e6a37784dcb9cc53cd..280842ecc3874e9ba634fd997b194fd940635a57 100644
--- a/applications/solvers/compressible/rhoPisoFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPisoFoam/pEqn.H
@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -34,7 +34,7 @@ if (transonic)
 }
 else
 {
-    phi = 
+    phi =
         fvc::interpolate(rho)*
         (
             (fvc::interpolate(U) & mesh.Sf())
diff --git a/applications/solvers/compressible/rhoPisoFoam/rhoPisoFoam.C b/applications/solvers/compressible/rhoPisoFoam/rhoPisoFoam.C
index c349e7520e6ad34fd55550c77574107a6ff9fd13..2d9cccd29e6515482f043586d98217c606a8830c 100644
--- a/applications/solvers/compressible/rhoPisoFoam/rhoPisoFoam.C
+++ b/applications/solvers/compressible/rhoPisoFoam/rhoPisoFoam.C
@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -77,7 +77,7 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
index eee3959a4e0850702aa4ef1a47052aece69971fb..1177cba2a4b86adcdc0081e5305b9095214976a9 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -15,11 +16,12 @@
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
 
 
     Info<< "Reading field U\n" << endl;
@@ -36,7 +38,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
 
     label pRefCell = 0;
@@ -56,7 +58,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
index 8eb03f95eb46b687a79f9f68c09152bbd269224a..605b8820d1816daeaa6a6b4a2f20965e388a1c8f 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
@@ -14,5 +14,5 @@
     eqnResidual = hEqn.solve().initialResidual();
     maxResidual = max(eqnResidual, maxResidual);
 
-    thermo->correct();
+    thermo.correct();
 }
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
index ae41da9e36fd7af8dd172b163b16dd8829a93b8b..4d8e010f7e352464dd92558f23342c10774f6575 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
@@ -65,10 +65,10 @@ bound(p, pMin);
 // to obey overall mass continuity
 if (closedVolume)
 {
-    p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
-        /fvc::domainIntegrate(thermo->psi());
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
 }
 
-rho = thermo->rho();
+rho = thermo.rho();
 rho.relax();
 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C b/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
index 0cd144fce23508b23d803c7abf25510f67f17f99..fbc27d884603ac2c9c80a702d0c7dc9286bd6211 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
@@ -27,12 +27,12 @@ Application
 
 Description
     Steady-state solver for turbulent flow of compressible fluids with
-    implicit or explicit porosity treatment
+    RANS turbulence modelling, and implicit or explicit porosity treatment
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "RASModel.H"
 #include "porousZones.H"
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index 9c4f446c37ff45885232c1219825f8840fe9cbdb..690b53760d75b22618fc9cd0c1358c015d867854 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -15,12 +16,12 @@
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
-
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
 
     Info<< "Reading field U\n" << endl;
     volVectorField U
@@ -56,7 +57,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
index 6a87bbdf11941f36afe9fb7a43186f92cf554e0f..e299d99f83c9b1e0c0ab95b35d015300f4f18d7b 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
@@ -14,5 +14,5 @@
     eqnResidual = hEqn.solve().initialResidual();
     maxResidual = max(eqnResidual, maxResidual);
 
-    thermo->correct();
+    thermo.correct();
 }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index e4fbe15aa485869258c231d8aef86f4e49a5c096..f6a433fd6164d0f803098051d4a9c44f4fafa38c 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn().A();
 U = rUA*UEqn().H();
@@ -11,7 +11,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())*(fvc::interpolate(U) & mesh.Sf())
+        fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
     );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -82,7 +82,7 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-rho = thermo->rho();
+rho = thermo.rho();
 rho.relax();
 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
 
@@ -95,6 +95,6 @@ bound(p, pMin);
 // to obey overall mass continuity
 if (closedVolume)
 {
-    p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
-        /fvc::domainIntegrate(thermo->psi());
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
 }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
index bf6a9f24f3386c7c1f93237fc244c30f85e87951..099257e5ae40a7deec364106327f1a82205a935a 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
@@ -26,13 +26,13 @@ Application
     rhoSimpleFoam
 
 Description
-    Steady-state SIMPLE solver for laminar or turbulent flow of
+    Steady-state SIMPLE solver for laminar or turbulent RANS flow of
     compressible fluids.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "RASModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/compressible/sonicDyMFoam/compressibleContinuityErrs.H b/applications/solvers/compressible/sonicDyMFoam/compressibleContinuityErrs.H
index eec8485329cfcb0956c0b49b9a6108ead3cc797e..1d869cc58dae738d345b2448f7dfdcf7a63c36cb 100644
--- a/applications/solvers/compressible/sonicDyMFoam/compressibleContinuityErrs.H
+++ b/applications/solvers/compressible/sonicDyMFoam/compressibleContinuityErrs.H
@@ -1,5 +1,5 @@
 {
-#   include "rhoEqn.H"
+    #include "rhoEqn.H"
 }
 {
     scalar sumLocalContErr =
diff --git a/applications/solvers/compressible/sonicDyMFoam/createFields.H b/applications/solvers/compressible/sonicDyMFoam/createFields.H
index eb45a7c23af7826cd0843adead9e5a3ea74ccfbf..bbb5d105269512474d7753bb93ef97affda1e8cc 100644
--- a/applications/solvers/compressible/sonicDyMFoam/createFields.H
+++ b/applications/solvers/compressible/sonicDyMFoam/createFields.H
@@ -1,13 +1,14 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
-    const volScalarField& psi = thermo->psi();
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -17,7 +18,7 @@
             runTime.timeName(),
             mesh
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "Reading field U\n" << endl;
@@ -45,7 +46,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
diff --git a/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C b/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C
index 57b15a0f9ba4e546d59ad17eb9a7129989ca4ad1..d86e41913e94044061c9084c219c934c2df3c63a 100644
--- a/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C
+++ b/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C
@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "motionSolver.H"
 
@@ -84,8 +84,8 @@ int main(int argc, char *argv[])
             surfaceScalarField phid
             (
                 "phid",
-                fvc::interpolate(psi)*
-                (
+                fvc::interpolate(psi)
+               *(
                     (fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)
                 )
             );
diff --git a/applications/solvers/compressible/sonicFoam/createFields.H b/applications/solvers/compressible/sonicFoam/createFields.H
index eb45a7c23af7826cd0843adead9e5a3ea74ccfbf..bbb5d105269512474d7753bb93ef97affda1e8cc 100644
--- a/applications/solvers/compressible/sonicFoam/createFields.H
+++ b/applications/solvers/compressible/sonicFoam/createFields.H
@@ -1,13 +1,14 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
-    const volScalarField& psi = thermo->psi();
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -17,7 +18,7 @@
             runTime.timeName(),
             mesh
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "Reading field U\n" << endl;
@@ -45,7 +46,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
diff --git a/applications/solvers/compressible/sonicFoam/hEqn.H b/applications/solvers/compressible/sonicFoam/hEqn.H
index 5cb4c48c0883e6fce99b9e9a2e76b277e9cbfef7..baa2dab34334e1d942f17c29e8e3675747def570 100644
--- a/applications/solvers/compressible/sonicFoam/hEqn.H
+++ b/applications/solvers/compressible/sonicFoam/hEqn.H
@@ -8,5 +8,5 @@
             DpDt
         );
 
-        thermo->correct();
+        thermo.correct();
 }
diff --git a/applications/solvers/compressible/sonicFoam/sonicFoam.C b/applications/solvers/compressible/sonicFoam/sonicFoam.C
index 6ff79607986802a024e450042be2ad286d6d5b12..c9eda00fa34e1736a677d0d36aa4092a1a3b79a6 100644
--- a/applications/solvers/compressible/sonicFoam/sonicFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicFoam.C
@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -80,7 +80,7 @@ int main(int argc, char *argv[])
             surfaceScalarField phid
             (
                 "phid",
-                fvc::interpolate(thermo->psi())
+                fvc::interpolate(psi)
                *(
                    (fvc::interpolate(U) & mesh.Sf())
                  + fvc::ddtPhiCorr(rUA, rho, U, phi)
diff --git a/applications/solvers/compressible/sonicLiquidFoam/createFields.H b/applications/solvers/compressible/sonicLiquidFoam/createFields.H
index f419234c7dbf28ee9231ce36837df22d996912ae..1268cbb7b0657b70f8ffa4131204ad516b2e1b79 100644
--- a/applications/solvers/compressible/sonicLiquidFoam/createFields.H
+++ b/applications/solvers/compressible/sonicLiquidFoam/createFields.H
@@ -41,4 +41,4 @@
     );
 
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
diff --git a/applications/solvers/compressible/sonicLiquidFoam/sonicLiquidFoam.C b/applications/solvers/compressible/sonicLiquidFoam/sonicLiquidFoam.C
index 0324da8299896476159dca790e388716a8ca8ac3..763e03a4b9759e3210379affa4bc44341e7ea848 100644
--- a/applications/solvers/compressible/sonicLiquidFoam/sonicLiquidFoam.C
+++ b/applications/solvers/compressible/sonicLiquidFoam/sonicLiquidFoam.C
@@ -37,16 +37,15 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readThermodynamicProperties.H"
+    #include "readTransportProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
 
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readThermodynamicProperties.H"
-#   include "readTransportProperties.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -54,10 +53,10 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
 
         fvVectorMatrix UEqn
         (
@@ -79,8 +78,8 @@ int main(int argc, char *argv[])
             surfaceScalarField phid
             (
                 "phid",
-                psi*
-                (
+                psi
+               *(
                     (fvc::interpolate(U) & mesh.Sf())
                   + fvc::ddtPhiCorr(rUA, rho, U, phi)
                 )
@@ -100,7 +99,7 @@ int main(int argc, char *argv[])
 
             phi += pEqn.flux();
 
-#           include "compressibleContinuityErrs.H"
+            #include "compressibleContinuityErrs.H"
 
             U -= rUA*fvc::grad(p);
             U.correctBoundaryConditions();
diff --git a/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C b/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C
index b7efa050deee79004a5ec3076695c172fdc7d91a..b5055812c40e362d420c08ab862c9db91e2469c5 100644
--- a/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C
@@ -81,7 +81,7 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 
diff --git a/applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H
index bd0d71dc57105f897441cde7bfdfc5e16842ce16..b5c84e435e3b2a95ad2cb8b67e215559d49758a0 100644
--- a/applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H
@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicRhoThermo> thermo
+    autoPtr<basicRhoThermo> pThermo
     (
         basicRhoThermo::New(mesh)
     );
+    basicRhoThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -15,12 +16,12 @@
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
-    const volScalarField& psi = thermo->psi();
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
 
 
     Info<< "Reading field U\n" << endl;
@@ -37,7 +38,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
 
     Info<< "Creating turbulence model\n" << endl;
@@ -48,7 +49,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
@@ -59,7 +60,7 @@
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
     );
 
-    thermo->correct();
+    thermo.correct();
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
diff --git a/applications/solvers/heatTransfer/buoyantPisoFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantPisoFoam/hEqn.H
index f1a87c6a798229100ff3493a194d73ccfdd601b9..3125cc3ffa86ce120e7dbbf774c9b46941105418 100644
--- a/applications/solvers/heatTransfer/buoyantPisoFoam/hEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPisoFoam/hEqn.H
@@ -11,5 +11,5 @@
     hEqn.relax();
     hEqn.solve();
 
-    thermo->correct();
+    thermo.correct();
 }
diff --git a/applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H
index cfc8287cadd342297c65423a92f6571f3595f3ae..c954c0ecb193a86033981bc7123725bc9335a9cc 100644
--- a/applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H
@@ -1,11 +1,11 @@
 {
     bool closedVolume = p.needReference();
 
-    rho = thermo->rho();
+    rho = thermo.rho();
 
     // Thermodynamic density needs to be updated by psi*d(p) after the
     // pressure solution - done in 2 parts. Part 1:
-    thermo->rho() -= psi*p;
+    thermo.rho() -= psi*p;
 
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
@@ -48,7 +48,7 @@
     }
 
     // Second part of thermodynamic density update
-    thermo->rho() += psi*p;
+    thermo.rho() += psi*p;
 
     U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
     U.correctBoundaryConditions();
@@ -65,7 +65,7 @@
         p +=
             (initialMass - fvc::domainIntegrate(psi*p))
             /fvc::domainIntegrate(psi);
-        thermo->rho() = psi*p;
+        thermo.rho() = psi*p;
         rho += (initialMass - fvc::domainIntegrate(rho))/totalVolume;
     }
 }
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
index b2a8c5c374211719be84f7e03b6340f55a29e2e9..52f915bd7abe4b59a89aa7b30f1e40de0279e07c 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "RASModel.H"
 #include "fixedGradientFvPatchFields.H"
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index 0d9d9d5ec7fc667216c1af6937a24f50d7905986..b304ace1046bbb82fff0318762b17848835cce1d 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -15,11 +16,12 @@
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
 
 
     Info<< "Reading field U\n" << endl;
@@ -36,7 +38,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
 
     Info<< "Creating turbulence model\n" << endl;
@@ -47,11 +49,11 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
-    thermo->correct();
+    thermo.correct();
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H
index dad078c2193a6797d2f08ed9751cdac4ca4ff157..0c507ce3d18bff29619626d62d5570c67e415cc2 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H
@@ -14,5 +14,5 @@
     eqnResidual = hEqn.solve().initialResidual();
     maxResidual = max(eqnResidual, maxResidual);
 
-    thermo->correct();
+    thermo.correct();
 }
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/initConvergenceCheck.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/initConvergenceCheck.H
index b56197f22a50cfd07b04fc14d40b9a5454da8c5b..c920b6708d044d177227d79973ed61bc7f2658fc 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/initConvergenceCheck.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/initConvergenceCheck.H
@@ -1,4 +1,4 @@
-// initialize values for convergence checks
+    // initialize values for convergence checks
 
     scalar eqnResidual = 1, maxResidual = 0;
     scalar convergenceCriterion = 0;
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index d592bc3c43eb47f568b79f123f27dd95654fe6ff..e680e9421aee62a3dac1d09b8174005925399050 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -1,5 +1,5 @@
 {
-    rho = thermo->rho();
+    rho = thermo.rho();
 
     volScalarField rUA = 1.0/UEqn().A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
@@ -39,8 +39,8 @@
             // to obey overall mass continuity
             if (closedVolume)
             {
-                p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
-                    /fvc::domainIntegrate(thermo->psi());
+                p += (initialMass - fvc::domainIntegrate(psi*p))
+                    /fvc::domainIntegrate(psi);
             }
 
             // Calculate the conservative fluxes
@@ -58,7 +58,7 @@
 
     #include "continuityErrs.H"
 
-    rho = thermo->rho();
+    rho = thermo.rho();
     rho.relax();
     Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
         << endl;
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
index 5996e6be23b5fe4d1845f46ea34965ada35422dd..0862d20bb2ce9fe649be9b286df8b149ae3e024f 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "RASModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "radiationModel.H"
@@ -61,7 +61,7 @@ int main(int argc, char *argv[])
 #       include "readSIMPLEControls.H"
 #       include "initConvergenceCheck.H"
 
-        pd.storePrevIter();
+        p.storePrevIter();
         rho.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H
index 31462927b94bdd60bba71db3718e6d98a7e9defb..24b17645d2c3f62b38cdca9a4c7b908a8e08931f 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H
@@ -7,7 +7,7 @@
      ==
         fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
       - p*fvc::div(phi/fvc::interpolate(rho))
-      + radiation->Sh(thermo())
+      + radiation->Sh(thermo)
     );
 
     hEqn.relax();
@@ -15,7 +15,7 @@
     eqnResidual = hEqn.solve().initialResidual();
     maxResidual = max(eqnResidual, maxResidual);
 
-    thermo->correct();
+    thermo.correct();
 
     radiation->correct();
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
index 7ad940bcff1bf5910c2c78253d5068b78a80f768..54e1c8cb1840d4719b1f2f360a7ab3d4ba81b052 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "regionProperties.H"
@@ -42,37 +42,36 @@ Description
 
 int main(int argc, char *argv[])
 {
-
-#   include "setRootCase.H"
-#   include "createTime.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
 
     regionProperties rp(runTime);
 
-#   include "createFluidMeshes.H"
-#   include "createSolidMeshes.H"
+    #include "createFluidMeshes.H"
+    #include "createSolidMeshes.H"
 
-#   include "createFluidFields.H"
-#   include "createSolidFields.H"
+    #include "createFluidFields.H"
+    #include "createSolidFields.H"
 
-#   include "initContinuityErrs.H"
+    #include "initContinuityErrs.H"
 
-#   include "readTimeControls.H"
+    #include "readTimeControls.H"
 
     if (fluidRegions.size())
     {
-#       include "compressibleMultiRegionCourantNo.H"
-#       include "setInitialDeltaT.H"
+        #include "compressibleMultiRegionCourantNo.H"
+        #include "setInitialDeltaT.H"
     }
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
-#       include "readPIMPLEControls.H"
+        #include "readTimeControls.H"
+        #include "readPIMPLEControls.H"
 
         if (fluidRegions.size())
         {
-#           include "compressibleMultiRegionCourantNo.H"
-#           include "setDeltaT.H"
+            #include "compressibleMultiRegionCourantNo.H"
+            #include "setDeltaT.H"
         }
 
         runTime++;
@@ -83,8 +82,8 @@ int main(int argc, char *argv[])
         {
             forAll(fluidRegions, i)
             {
-#               include "setRegionFluidFields.H"
-#               include "storeOldFluidFields.H"
+                #include "setRegionFluidFields.H"
+                #include "storeOldFluidFields.H"
             }
         }
 
@@ -96,18 +95,18 @@ int main(int argc, char *argv[])
             {
                 Info<< "\nSolving for fluid region "
                     << fluidRegions[i].name() << endl;
-#               include "setRegionFluidFields.H"
-#               include "readFluidMultiRegionPIMPLEControls.H"
-#               include "solveFluid.H"
+                #include "setRegionFluidFields.H"
+                #include "readFluidMultiRegionPIMPLEControls.H"
+                #include "solveFluid.H"
             }
 
             forAll(solidRegions, i)
             {
                 Info<< "\nSolving for solid region "
                     << solidRegions[i].name() << endl;
-#               include "setRegionSolidFields.H"
-#               include "readSolidMultiRegionPIMPLEControls.H"
-#               include "solveSolid.H"
+                #include "setRegionSolidFields.H"
+                #include "readSolidMultiRegionPIMPLEControls.H"
+                #include "solveSolid.H"
             }
         }
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index fd9624685f8fcc07cc4f82d014676fb469bbd064..7aa01dee7939de0bf5ce3595f0cefc73a5d7653f 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -1,5 +1,5 @@
     // Initialise fluid field pointer lists
-    PtrList<basicThermo> thermoFluid(fluidRegions.size());
+    PtrList<basicPsiThermo> thermoFluid(fluidRegions.size());
     PtrList<volScalarField> rhoFluid(fluidRegions.size());
     PtrList<volScalarField> KFluid(fluidRegions.size());
     PtrList<volVectorField> UFluid(fluidRegions.size());
@@ -19,7 +19,7 @@
         thermoFluid.set
         (
             i,
-            basicThermo::New(fluidRegions[i]).ptr()
+            basicPsiThermo::New(fluidRegions[i]).ptr()
         );
 
         Info<< "    Adding to rhoFluid\n" << endl;
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 72c9bbc4faf7e81c0fd759466c210a536929ff2c..29aea044664c8a7c2045525c5e94d94e0c0b1e48 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -1,6 +1,6 @@
     const fvMesh& mesh = fluidRegions[i];
 
-    basicThermo& thermo = thermoFluid[i];
+    basicPsiThermo& thermo = thermoFluid[i];
     volScalarField& rho = rhoFluid[i];
     volScalarField& K = KFluid[i];
     volVectorField& U = UFluid[i];
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFoam/Make/options
index 7d9fba61432e98b3f2cf2b390e4dcc547d488f8a..d10253a8fa3615ecf3f0d435390af869fb16a743 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/Make/options
@@ -12,7 +12,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude
@@ -31,7 +31,7 @@ EXE_LIBS = \
     -lsolids \
     -lsolidMixture \
     -lthermophysicalFunctions \
-    -lcombustionThermophysicalModels \
+    -lreactionThermophysicalModels \
     -lchemistryModel \
     -lradiation \
     -lODE
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/additionalOutput.H b/applications/solvers/lagrangian/reactingParcelFoam/additionalOutput.H
deleted file mode 100644
index 9edd35eb7cab9acd294ba6a9e73b6ee8087b1c0d..0000000000000000000000000000000000000000
--- a/applications/solvers/lagrangian/reactingParcelFoam/additionalOutput.H
+++ /dev/null
@@ -1,48 +0,0 @@
-{
-    tmp<volScalarField> tdQ
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "dQ",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            mesh,
-            dimensionedScalar
-            (
-                "zero",
-                dimensionSet(1, -3, -1, 0, 0, 0, 0),
-                0.0
-            )
-        )
-    );
-
-    scalarField& dQ = tdQ();
-
-    scalarField cp(dQ.size(), 0.0);
-
-    forAll(Y, i)
-    {
-        volScalarField RRi = chemistry.RR(i);
-
-        forAll(h, celli)
-        {
-            scalar Ti = T[celli];
-            cp[celli] += Y[i][celli]*chemistry.specieThermo()[i].Cp(Ti);
-            scalar hi = chemistry.specieThermo()[i].h(Ti);
-            scalar RR = RRi[celli];
-            dQ[celli] -= hi*RR;
-        }
-    }
-
-    forAll(dQ, celli)
-    {
-        dQ[celli] /= cp[celli];
-    }
-
-    tdQ().write();
-}
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
index fcab29c923c274b45ec86647cb3b5eb3a54ecc7b..4ae0633ab7bc74f993286571c10f3c52f645e7e4 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
@@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-BasicReactingCloud<specieReactingProperties> parcels
+BasicReactingCloud<gasThermoPhysics> parcels
 (
     "reactingCloud1",
     rho,
     U,
     g,
-    thermo()
+    thermo
 );
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
index 8f0d59e77d4d2e88e6080b7b029dcb1facefa450..cad1330c68da749ecefa26e80375ef7d43e6ca25 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
@@ -1,19 +1,22 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<hCombustionThermo> thermo
+    autoPtr<psiChemistryModel> pChemistry
     (
-        hCombustionThermo::New(mesh)
+        psiChemistryModel::New(mesh)
     );
+    psiChemistryModel& chemistry = pChemistry();
 
-    combustionMixture& composition = thermo->composition();
+    hCombustionThermo& thermo = chemistry.thermo();
+
+    basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
 
-    word inertSpecie(thermo->lookup("inertSpecie"));
+    word inertSpecie(thermo.lookup("inertSpecie"));
 
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
-    const volScalarField& T = thermo->T();
-    const volScalarField& psi = thermo->psi();
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& T = thermo.T();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -25,7 +28,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "\nReading field U\n" << endl;
@@ -42,7 +45,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
     DimensionedField<scalar, volMesh> kappa
     (
@@ -66,7 +69,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
@@ -74,13 +77,6 @@
     volScalarField DpDt =
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
-    Info << "Constructing chemical mechanism" << endl;
-    chemistryModel chemistry
-    (
-        thermo(),
-        rho
-    );
-
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
     forAll (Y, i)
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/hEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/hEqn.H
index 7909d8c67d2b3962b6da9a75da8637ccd4673ed0..020c9bdb118bb15849a81c4febf0f712b16cf0f7 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/hEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/hEqn.H
@@ -7,14 +7,14 @@
      ==
         DpDt
      +  parcels.Sh()
-     +  radiation->Sh(thermo())
+     +  radiation->Sh(thermo)
     );
 
     hEqn.relax();
 
     hEqn.solve();
 
-    thermo->correct();
+    thermo.correct();
 
     radiation->correct();
 }
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index 63370f0cd41ea582ddaa646e9861f7d362293f00..d27afcce8ac2464b63b09e05645bcb5de6aa1282 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -37,8 +37,8 @@ if (transonic)
 else
 {
     phi =
-        fvc::interpolate(rho)*
-        (
+        fvc::interpolate(rho)
+       *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index d93b8236612834565da42ca395fc2cbf8ac9bf37..f60de4a3690648beffdbd72e14e1077c1902180b 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -35,9 +35,9 @@ Description
 #include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "BasicReactingCloud.H"
-#include "chemistryModel.H"
+#include "psiChemistryModel.H"
 #include "chemistrySolver.H"
-#include "thermoPhsyicsTypes.H"
+#include "thermoPhysicsTypes.H"
 #include "radiationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -100,12 +100,9 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
-        if (runTime.write())
-        {
-            #include "additionalOutput.H"
-        }
+        runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H
index 1b2882a2b7b7e6bd09521babeb2cda56feb5f0aa..131c144114eb2b6f04bfd0798ec2facac0516ef0 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/Make/options b/applications/solvers/lagrangian/trackedReactingParcelFoam/Make/options
index 76bc8b084f460664af2e459e21f691c45be205e4..76ae4faa08382045bd9c93b3a3c908bb1a404a55 100644
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/trackedReactingParcelFoam/Make/options
@@ -1,4 +1,7 @@
+REACTINGPARCELFOAM = ../reactingParcelFoam
+
 EXE_INC = \
+    -I$(REACTINGPARCELFOAM) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/UEqn.H
index c82307305b8cba5e0877a7f0a13784ab2f66978b..c3716a63818db9de23ba724c8dbfed1f3eff5be5 100644
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/trackedReactingParcelFoam/UEqn.H
@@ -5,7 +5,7 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho.dimensionedInternalField()*g
-      + reactingParcels.SU()
+      + parcels.SU()
     );
 
     UEqn.relax();
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/YEqn.H
index 8f77f23bc2f3f3799ed84830c75eab4bbd5b438b..2fc0f375664ed2f6c0b80f0ee7910881134020c8 100644
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/trackedReactingParcelFoam/YEqn.H
@@ -26,7 +26,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->muEff(), Yi)
               ==
-                reactingParcels.Srho(i)
+                parcels.Srho(i)
               + kappa*chemistry.RR(i)().dimensionedInternalField()
               + pointMassSources.Su(i)
             );
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/additionalOutput.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/additionalOutput.H
deleted file mode 100644
index 9edd35eb7cab9acd294ba6a9e73b6ee8087b1c0d..0000000000000000000000000000000000000000
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/additionalOutput.H
+++ /dev/null
@@ -1,48 +0,0 @@
-{
-    tmp<volScalarField> tdQ
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "dQ",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            mesh,
-            dimensionedScalar
-            (
-                "zero",
-                dimensionSet(1, -3, -1, 0, 0, 0, 0),
-                0.0
-            )
-        )
-    );
-
-    scalarField& dQ = tdQ();
-
-    scalarField cp(dQ.size(), 0.0);
-
-    forAll(Y, i)
-    {
-        volScalarField RRi = chemistry.RR(i);
-
-        forAll(h, celli)
-        {
-            scalar Ti = T[celli];
-            cp[celli] += Y[i][celli]*chemistry.specieThermo()[i].Cp(Ti);
-            scalar hi = chemistry.specieThermo()[i].h(Ti);
-            scalar RR = RRi[celli];
-            dQ[celli] -= hi*RR;
-        }
-    }
-
-    forAll(dQ, celli)
-    {
-        dQ[celli] /= cp[celli];
-    }
-
-    tdQ().write();
-}
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/chemistry.H
deleted file mode 100644
index 07b1e9953b0db867186f6c668d27a9415a26c265..0000000000000000000000000000000000000000
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/chemistry.H
+++ /dev/null
@@ -1,25 +0,0 @@
-{
-    Info << "Solving chemistry" << endl;
-
-    chemistry.solve
-    (
-        runTime.value() - runTime.deltaT().value(),
-        runTime.deltaT().value()
-    );
-
-    // turbulent time scale
-    if (turbulentReaction)
-    {
-        DimensionedField<scalar, volMesh> tk =
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
-        DimensionedField<scalar, volMesh> tc =
-            chemistry.tc()().dimensionedInternalField();
-
-        // Chalmers PaSR model
-        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
-    }
-    else
-    {
-        kappa = 1.0;
-    }
-}
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/createClouds.H
index aaf91c81cc6ce8c5ba6248863458b17d98b4fd79..3740171bebae348b385e6fcd08bb94482ed6e8e8 100644
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/trackedReactingParcelFoam/createClouds.H
@@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-BasicTrackedReactingCloud<specieReactingProperties> reactingParcels
+BasicTrackedReactingCloud<gasThermoPhysics> parcels
 (
     "reactingCloud1",
     rho,
     U,
     g,
-    thermo()
+    thermo
 );
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/createFields.H
deleted file mode 100644
index b9cea162f4a863babc5ef3536cf68e2451c84297..0000000000000000000000000000000000000000
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/createFields.H
+++ /dev/null
@@ -1,111 +0,0 @@
-    Info<< "Reading thermophysical properties" << nl << endl;
-
-    autoPtr<hCombustionThermo> thermo
-    (
-        hCombustionThermo::New(mesh)
-    );
-
-    combustionMixture& composition = thermo->composition();
-    PtrList<volScalarField>& Y = composition.Y();
-
-    word inertSpecie(thermo->lookup("inertSpecie"));
-
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
-    const volScalarField& T = thermo->T();
-    const volScalarField& psi = thermo->psi();
-
-    volScalarField rho
-    (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        thermo->rho()
-    );
-
-    Info<< "Reading field U" << nl << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "compressibleCreatePhi.H"
-
-    DimensionedField<scalar, volMesh> kappa
-    (
-        IOobject
-        (
-            "kappa",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimless, 0.0)
-    );
-
-    Info<< "Creating turbulence model" << nl << endl;
-    autoPtr<compressible::turbulenceModel> turbulence
-    (
-        compressible::turbulenceModel::New
-        (
-            rho,
-            U,
-            phi,
-            thermo()
-        )
-    );
-
-    Info<< "Creating field DpDt" << nl << endl;
-    volScalarField DpDt =
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
-
-    Info << "Constructing chemical mechanism" << nl << endl;
-    chemistryModel chemistry
-    (
-        thermo(),
-        rho
-    );
-
-    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-    forAll (Y, i)
-    {
-        fields.add(Y[i]);
-    }
-    fields.add(h);
-
-    Info<< "Creating porous zones" << nl << endl;
-    porousZones pZones(mesh);
-    Switch pressureImplicitPorosity(false);
-
-    label nUCorr = 0;
-    if (pZones.size())
-    {
-        // nUCorrectors for pressureImplicitPorosity
-        if (mesh.solutionDict().subDict("PISO").found("nUCorrectors"))
-        {
-            mesh.solutionDict().subDict("PISO").lookup("nUCorrectors")
-                >> nUCorr;
-        }
-
-        if (nUCorr > 0)
-        {
-            pressureImplicitPorosity = true;
-        }
-    }
-
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/createPorousZones.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/createPorousZones.H
new file mode 100644
index 0000000000000000000000000000000000000000..caed85ecba899c76d9f3dcce50d474ef1114e1bf
--- /dev/null
+++ b/applications/solvers/lagrangian/trackedReactingParcelFoam/createPorousZones.H
@@ -0,0 +1,21 @@
+    Info<< "Creating porous zones" << nl << endl;
+
+    porousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);
+
+    label nUCorr = 0;
+    if (pZones.size())
+    {
+        // nUCorrectors for pressureImplicitPorosity
+        if (mesh.solutionDict().subDict("PISO").found("nUCorrectors"))
+        {
+            mesh.solutionDict().subDict("PISO").lookup("nUCorrectors")
+                >> nUCorr;
+        }
+
+        if (nUCorr > 0)
+        {
+            pressureImplicitPorosity = true;
+        }
+    }
+
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/hEqn.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/hEqn.H
deleted file mode 100644
index 32e71555369290be7d019af0f4d6ad15473ee064..0000000000000000000000000000000000000000
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/hEqn.H
+++ /dev/null
@@ -1,20 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + fvm::div(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        DpDt
-     +  reactingParcels.Sh()
-     +  radiation->Sh(thermo())
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo->correct();
-
-    radiation->correct();
-}
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/pEqn.H
index b6207ce0c4df93ff1272e297b0eb5c7553976c1c..f8d6f4e6383ad19f6934d46490577cf657a4363b 100644
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/trackedReactingParcelFoam/pEqn.H
@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 if (pressureImplicitPorosity)
 {
@@ -14,7 +14,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
 //          + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -40,7 +40,7 @@ if (transonic)
           + fvm::div(phid, p)
           - lapTerm()
          ==
-            reactingParcels.Srho()
+            parcels.Srho()
           + pointMassSources.Su()
         );
 
@@ -80,7 +80,7 @@ else
           + fvc::div(phi)
           - lapTerm()
          ==
-            reactingParcels.Srho()
+            parcels.Srho()
           + pointMassSources.Su()
         );
 
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/readChemistryProperties.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/readChemistryProperties.H
deleted file mode 100644
index 1a60e6fb34645a004fd39321f7a54d3bd5b45381..0000000000000000000000000000000000000000
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/readChemistryProperties.H
+++ /dev/null
@@ -1,22 +0,0 @@
-Info<< "Reading chemistry properties\n" << endl;
-
-IOdictionary chemistryProperties
-(
-    IOobject
-    (
-        "chemistryProperties",
-        runTime.constant(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE
-    )
-);
-
-Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
-
-dimensionedScalar Cmix("Cmix", dimless, 1.0);
-
-if (turbulentReaction)
-{
-    chemistryProperties.lookup("Cmix") >> Cmix;
-}
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/trackedReactingParcelFoam/rhoEqn.H
index faadc1ecfbb1c8347bc8e6b1b6102f1a0ab911c4..76b105474ba491ac489e77f559ad9673b7941618 100644
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/trackedReactingParcelFoam/rhoEqn.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,7 +36,7 @@ Description
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
-        reactingParcels.Srho()
+        parcels.Srho()
       + pointMassSources.Su()
     );
 }
diff --git a/applications/solvers/lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C b/applications/solvers/lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C
index 62af964eba7dad29b1c2517c025d02af76df1c38..a8338a5e34901dfac9b07ec37d32b3438838d896 100644
--- a/applications/solvers/lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C
@@ -32,9 +32,9 @@ Description
 #include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "BasicTrackedReactingCloud.H"
-#include "chemistryModel.H"
+#include "psiChemistryModel.H"
 #include "chemistrySolver.H"
-#include "reactingThermoTypes.H"
+#include "thermoPhysicsTypes.H"
 #include "radiationModel.H"
 #include "porousZones.H"
 #include "timeActivatedExplicitMulticomponentPointSource.H"
@@ -53,6 +53,7 @@ int main(int argc, char *argv[])
     #include "createRadiationModel.H"
     #include "createClouds.H"
     #include "createMulticomponentPointSources.H"
+    #include "createPorousZones.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
@@ -74,9 +75,9 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        reactingParcels.evolve();
+        parcels.evolve();
 
-        reactingParcels.info();
+        parcels.info();
 
         #include "chemistry.H"
         #include "rhoEqn.H"
@@ -100,12 +101,9 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
-        if (runTime.write())
-        {
-            #include "additionalOutput.H"
-        }
+        runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
diff --git a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/createFields.H b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/createFields.H
index c76a81ca4666d882442b8c8f09663be378a134d0..78b0e8b552f8a9b63fce0a2bf52b545200d944f4 100644
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/createFields.H
@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -15,7 +16,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "\nReading field U\n" << endl;
@@ -42,7 +43,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
@@ -55,6 +56,6 @@
         kinematicCloudName,
         rho,
         U,
-        thermo().mu(),
+        thermo.mu(),
         g
     );
diff --git a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C
index e4c2b396652630423cc6ed925e55a61b0203ce03..2039346411b21b7ac84a633416040c4e1ecc71d0 100644
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C
@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "basicKinematicCloud.H"
 
diff --git a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C b/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
index 66df8cfb9f260d31d788d24f20819ee5ed208629..af6de6def0bb1d6561b24052aee46724b67a2f54 100644
--- a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
+++ b/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
@@ -26,7 +26,7 @@ Application
     bubbleFoam
 
 Description
-    Solver for a system of 2 incompressible fluid phases with one phase 
+    Solver for a system of 2 incompressible fluid phases with one phase
     dispersed, e.g. gas bubbles in a liquid.
 
 \*---------------------------------------------------------------------------*/
@@ -40,16 +40,15 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -57,30 +56,30 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readBubbleFoamControls.H"
-#       include "CourantNo.H"
+        #include "readBubbleFoamControls.H"
+        #include "CourantNo.H"
 
-#       include "alphaEqn.H"
-#       include "liftDragCoeffs.H"
-#       include "UEqns.H"
+        #include "alphaEqn.H"
+        #include "liftDragCoeffs.H"
+        #include "UEqns.H"
 
         // --- PISO loop
 
         for (int corr=0; corr<nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
 
             if (correctAlpha)
             {
-#               include "alphaEqn.H"
+                #include "alphaEqn.H"
             }
         }
 
-#       include "DDtU.H"
+        #include "DDtU.H"
 
-#       include "kEpsilon.H"
+        #include "kEpsilon.H"
 
-#       include "write.H"
+        #include "write.H"
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
diff --git a/applications/solvers/multiphase/bubbleFoam/createFields.H b/applications/solvers/multiphase/bubbleFoam/createFields.H
index 875a4f1be13550588c393b2d52e17a19e3990a4a..9aca6933fc300f754f064f3b2ad95c05168e7fdd 100644
--- a/applications/solvers/multiphase/bubbleFoam/createFields.H
+++ b/applications/solvers/multiphase/bubbleFoam/createFields.H
@@ -141,8 +141,8 @@
         transportProperties.lookup("Ct")
     );
 
-#   include "createPhia.H"
-#   include "createPhib.H"
+    #include "createPhia.H"
+    #include "createPhib.H"
 
     surfaceScalarField phi
     (
@@ -157,7 +157,7 @@
     );
 
 
-#   include "createRASTurbulence.H"
+    #include "createRASTurbulence.H"
 
     Info<< "Calculating field DDtUa and DDtUb\n" << endl;
 
diff --git a/applications/solvers/multiphase/bubbleFoam/kEpsilon.H b/applications/solvers/multiphase/bubbleFoam/kEpsilon.H
index 016d1f4488c43c53bb22e11802bff96667712220..b1c8951b04ab389de71cbdf1a5efbc9193c734fb 100644
--- a/applications/solvers/multiphase/bubbleFoam/kEpsilon.H
+++ b/applications/solvers/multiphase/bubbleFoam/kEpsilon.H
@@ -9,7 +9,7 @@ if(turbulence)
     volScalarField G = 2*nutb*(tgradUb() && dev(symm(tgradUb())));
     tgradUb.clear();
 
-#   include "wallFunctions.H"
+    #include "wallFunctions.H"
 
     // Dissipation equation
     fvScalarMatrix epsEqn
@@ -22,7 +22,7 @@ if(turbulence)
        - fvm::Sp(C2*beta*epsilon/k, epsilon)
     );
 
-#   include "wallDissipation.H"
+    #include "wallDissipation.H"
 
     epsEqn.relax();
     epsEqn.solve();
@@ -48,8 +48,7 @@ if(turbulence)
     //- Re-calculate turbulence viscosity
     nutb = Cmu*sqr(k)/epsilon;
 
-#   include "wallViscosity.H"
-
+    #include "wallViscosity.H"
 }
 
 nuEffa = sqr(Ct)*nutb + nua;
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
index 12a8f4d49991b3e8f280bcef7f3274f05498a921..62aa7a076fbff1fb215ce1cda8f04cfd88119c13 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C
@@ -41,40 +41,39 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readThermodynamicProperties.H"
+    #include "readControls.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readThermodynamicProperties.H"
-#   include "readControls.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readControls.H"
-#       include "CourantNo.H"
-#       include "setDeltaT.H"
+        #include "readControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         for (int outerCorr=0; outerCorr<nOuterCorr; outerCorr++)
         {
-#           include "rhoEqn.H"
-#           include "gammaPsi.H"
-#           include "UEqn.H"
+            #include "rhoEqn.H"
+            #include "gammaPsi.H"
+            #include "UEqn.H"
 
             for (int corr=0; corr<nCorr; corr++)
             {
-#               include "pEqn.H"
+                #include "pEqn.H"
             }
         }
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/createFields.H b/applications/solvers/multiphase/cavitatingFoam/createFields.H
index dc2f5e696129edaac8b20e59effedc957c975a24..49c7de1473355e69429cf309091dc81aaea5a563 100644
--- a/applications/solvers/multiphase/cavitatingFoam/createFields.H
+++ b/applications/solvers/multiphase/cavitatingFoam/createFields.H
@@ -71,8 +71,8 @@
         mesh
     );
 
-#   include "createPhiv.H"
-#   include "compressibleCreatePhi.H"
+    #include "createPhiv.H"
+    #include "compressibleCreatePhi.H"
 
     Info<< "Reading transportProperties\n" << endl;
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index c22aa76479a0d503d22bc60af591a82a7733d850..979834717be78414443c08afdf542d4355eb5b64 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -24,7 +24,7 @@
 
     phiv -= phiGradp/rhof;
 
-#   include "resetPhivPatches.H"
+    #include "resetPhivPatches.H"
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
@@ -82,6 +82,5 @@
     Info<< "max-min rho: " << max(rho).value()
         << " " << min(rho).value() << endl;
 
-#   include "gammaPsi.H"
-
+    #include "gammaPsi.H"
 }
diff --git a/applications/solvers/multiphase/settlingFoam/createFields.H b/applications/solvers/multiphase/settlingFoam/createFields.H
index 7c8445a71d4465bae5514c86afbf83f6611ce057..916acde39ed5264a9f7bc765d3ec445b83519637 100644
--- a/applications/solvers/multiphase/settlingFoam/createFields.H
+++ b/applications/solvers/multiphase/settlingFoam/createFields.H
@@ -128,7 +128,7 @@
         alpha.boundaryField().types()
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
 
     label pRefCell = 0;
diff --git a/applications/solvers/multiphase/settlingFoam/settlingFoam.C b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
index ad8dbe3646cf2deb96e6affbdf75d155f2c90d37..e68ddbb0f0710ab911cc6855763e624c818c9883 100644
--- a/applications/solvers/multiphase/settlingFoam/settlingFoam.C
+++ b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
@@ -43,17 +43,15 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -61,27 +59,27 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
 
-#       include "calcVdj.H"
+        #include "calcVdj.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
-#       include "alphaEqn.H"
+        #include "alphaEqn.H"
 
-#       include "correctViscosity.H"
+        #include "correctViscosity.H"
 
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
-#       include "kEpsilon.H"
+        #include "kEpsilon.H"
 
         runTime.write();
 
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
index b61566d924db04be243167e57d2250845b3c0b3c..b94bf6c23e65cb80f5557f2bccd1ecd8e4c23946 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
@@ -40,7 +40,7 @@
         mesh
     );
 
-#   include "createPhi.H"
+    #include "createPhi.H"
 
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index ce3fb2aa551eeb45011a373277d6d8e169c01b7e..a28d929578081eb03bee4d6599b6575dc40341b7 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -47,59 +47,58 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "readPPProperties.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "readPPProperties.H"
-#   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readTwoPhaseEulerFoamControls.H"
-#       include "CourantNos.H"
-#       include "setDeltaT.H"
+        #include "readTwoPhaseEulerFoamControls.H"
+        #include "CourantNos.H"
+        #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "alphaEqn.H"
+        #include "alphaEqn.H"
 
-#       include "liftDragCoeffs.H"
+        #include "liftDragCoeffs.H"
 
-#       include "UEqns.H"
+        #include "UEqns.H"
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
 
             if (correctAlpha && corr<nCorr-1)
             {
-#               include "alphaEqn.H"
+                #include "alphaEqn.H"
             }
         }
 
-#       include "DDtU.H"
+        #include "DDtU.H"
 
-#       include "kEpsilon.H"
+        #include "kEpsilon.H"
 
         if (kineticTheory.on())
         {
             kineticTheory.solve();
             nuEffa += kineticTheory.mua()/rhoa;
         }
-#       include "write.H"
+        #include "write.H"
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
diff --git a/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C b/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C
index 0843f3acf3b72ea161cfb5958493373792787fac..b9310ffc92c5a51be55475d778465399da550a32 100644
--- a/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C
+++ b/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C
@@ -42,7 +42,7 @@ Description
 #include "incompressible/RAS/RASModel/RASModel.H"
 #include "incompressible/LES/LESModel/LESModel.H"
 
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "compressible/RAS/RASModel/RASModel.H"
 #include "compressible/LES/LESModel/LESModel.H"
 
@@ -194,7 +194,7 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     }
     else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
     {
-        autoPtr<basicThermo> thermo(basicThermo::New(mesh));
+        autoPtr<basicPsiThermo> thermo(basicPsiThermo::New(mesh));
 
         volScalarField rho
         (
diff --git a/applications/utilities/postProcessing/velocityField/Mach/Mach.C b/applications/utilities/postProcessing/velocityField/Mach/Mach.C
index daa4aa2ccf3d7aef33eef506761f46cb7c168a1c..2e7c517e9d0358f7bbdb37ea0e121384b755f356 100644
--- a/applications/utilities/postProcessing/velocityField/Mach/Mach.C
+++ b/applications/utilities/postProcessing/velocityField/Mach/Mach.C
@@ -33,7 +33,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "calc.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
@@ -66,9 +66,9 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
         if (isFile(runTime.constantPath()/"thermophysicalProperties"))
         {
             // thermophysical Mach
-            autoPtr<basicThermo> thermo
+            autoPtr<basicPsiThermo> thermo
             (
-                basicThermo::New(mesh)
+                basicPsiThermo::New(mesh)
             );
 
             volScalarField Cp = thermo->Cp();
diff --git a/applications/utilities/postProcessing/velocityField/Mach/thermophysicalMach.H b/applications/utilities/postProcessing/velocityField/Mach/thermophysicalMach.H
index 4e4ea81b3d87463b4c8d91fb115e0ea465e4b41b..7283802f7790a6c0cf209195a9273e51697fc017 100644
--- a/applications/utilities/postProcessing/velocityField/Mach/thermophysicalMach.H
+++ b/applications/utilities/postProcessing/velocityField/Mach/thermophysicalMach.H
@@ -18,9 +18,9 @@
         {
             volVectorField U(Uheader, mesh);
 
-            autoPtr<basicThermo> thermo
+            autoPtr<basicPsiThermo> thermo
             (
-                basicThermo::New(mesh)
+                basicPsiThermo::New(mesh)
             );
 
             volScalarField Cp = thermo->Cp();
diff --git a/applications/utilities/postProcessing/velocityField/Pe/Pe.C b/applications/utilities/postProcessing/velocityField/Pe/Pe.C
index 905074131d9f32b2d3d199dc8f12aae410343cb1..05c9ca1a0d2488dcd717df604e96f230f9cb9c7b 100644
--- a/applications/utilities/postProcessing/velocityField/Pe/Pe.C
+++ b/applications/utilities/postProcessing/velocityField/Pe/Pe.C
@@ -39,7 +39,7 @@ Description
 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
 #include "incompressible/RAS/RASModel/RASModel.H"
 #include "incompressible/LES/LESModel/LESModel.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "compressible/RAS/RASModel/RASModel.H"
 #include "compressible/LES/LESModel/LESModel.H"
 
@@ -204,7 +204,7 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
             {
                 IOdictionary RASProperties(RASPropertiesHeader);
 
-                autoPtr<basicThermo> thermo(basicThermo::New(mesh));
+                autoPtr<basicPsiThermo> thermo(basicPsiThermo::New(mesh));
 
                 volScalarField rho
                 (
@@ -252,7 +252,7 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
             {
                 IOdictionary LESProperties(LESPropertiesHeader);
 
-                autoPtr<basicThermo> thermo(basicThermo::New(mesh));
+                autoPtr<basicPsiThermo> thermo(basicPsiThermo::New(mesh));
 
                 volScalarField rho
                 (
diff --git a/src/finiteVolume/cfdTools/compressible/compressibleContinuityErrs.H b/src/finiteVolume/cfdTools/compressible/compressibleContinuityErrs.H
index b846a36cde00196ab2210bc8d0ad705887ebc4f3..9500b6d6713d3f17b48f5a242680d34155f28b33 100644
--- a/src/finiteVolume/cfdTools/compressible/compressibleContinuityErrs.H
+++ b/src/finiteVolume/cfdTools/compressible/compressibleContinuityErrs.H
@@ -33,15 +33,11 @@ Description
 {
     dimensionedScalar totalMass = fvc::domainIntegrate(rho);
 
-    scalar sumLocalContErr = 
-    (
-        fvc::domainIntegrate(mag(rho - thermo->rho()))/totalMass
-    ).value();
-
-    scalar globalContErr = 
-    (
-        fvc::domainIntegrate(rho - thermo->rho())/totalMass
-    ).value();
+    scalar sumLocalContErr =
+        (fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass).value();
+
+    scalar globalContErr =
+        (fvc::domainIntegrate(rho - thermo.rho())/totalMass).value();
 
     cumulativeContErr += globalContErr;
 
diff --git a/src/lagrangian/dieselSpray/spray/spray.C b/src/lagrangian/dieselSpray/spray/spray.C
index 4845d8a2e546e6f572a7bfede6c5f5d6665e7483..07dbcdbe7ef3665a01201daf4eb7cf31bb0bb405 100644
--- a/src/lagrangian/dieselSpray/spray/spray.C
+++ b/src/lagrangian/dieselSpray/spray/spray.C
@@ -57,7 +57,7 @@ Foam::spray::spray
     const volScalarField& p,
     const volScalarField& T,
     const basicMultiComponentMixture& composition,
-    const PtrList<specieProperties>& gasProperties,
+    const PtrList<gasThermoPhysics>& gasProperties,
     const dictionary&,
     const dictionary& environmentalProperties
 )
@@ -264,7 +264,7 @@ Foam::spray::spray
                 "spray::spray(const volVectorField& U, "
                 "const volScalarField& rho, const volScalarField& p, "
                 "const volScalarField& T, const combustionMixture& composition,"
-                "const PtrList<specieProperties>& gaseousFuelProperties, "
+                "const PtrList<gasThermoPhsyics>& gaseousFuelProperties, "
                 "const dictionary& thermophysicalProperties, "
                 "const dictionary& environmentalProperties)"
             )   << "spray::(...) only one wedgePolyPatch found. "
diff --git a/src/lagrangian/dieselSpray/spray/spray.H b/src/lagrangian/dieselSpray/spray/spray.H
index cc1e0ac382386529f6c343ae128d50fa0da79d30..f416e3cf4ea2f175588fa9c95d3633494f996420 100644
--- a/src/lagrangian/dieselSpray/spray/spray.H
+++ b/src/lagrangian/dieselSpray/spray/spray.H
@@ -38,10 +38,10 @@ Description
 #include "IOPtrList.H"
 #include "interpolation.H"
 #include "liquid.H"
-#include "sprayThermoTypes.H"
 #include "autoPtr.H"
 #include "liquidMixture.H"
 #include "Random.H"
+#include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -127,7 +127,7 @@ class spray
 
         // Composition properties
 
-            const PtrList<specieProperties>& gasProperties_;
+            const PtrList<gasThermoPhysics>& gasProperties_;
             const basicMultiComponentMixture& composition_;
 
             List<label> liquidToGasIndex_;
@@ -193,7 +193,7 @@ public:
             const volScalarField& p,
             const volScalarField& T,
             const basicMultiComponentMixture& composition,
-            const PtrList<specieProperties>& gasProperties,
+            const PtrList<gasThermoPhysics>& gasProperties,
             const dictionary& thermophysicalProperties,
             const dictionary& environmentalProperties
         );
@@ -257,7 +257,7 @@ public:
             inline const vector& g() const;
 
             inline const liquidMixture& fuels() const;
-            inline const PtrList<specieProperties>& gasProperties() const;
+            inline const PtrList<gasThermoPhysics>& gasProperties() const;
             inline const basicMultiComponentMixture& composition() const;
 
             inline const List<label>& liquidToGasIndex() const;
diff --git a/src/lagrangian/dieselSpray/spray/sprayI.H b/src/lagrangian/dieselSpray/spray/sprayI.H
index 17672c2a1faf369721674e4958f72bafefc31be5..86c063ca6f7dc5bdae3823837b0ff472918efb10 100644
--- a/src/lagrangian/dieselSpray/spray/sprayI.H
+++ b/src/lagrangian/dieselSpray/spray/sprayI.H
@@ -246,7 +246,7 @@ inline const liquidMixture& spray::fuels() const
 }
 
 
-inline const PtrList<specieProperties>& spray::gasProperties() const
+inline const PtrList<gasThermoPhysics>& spray::gasProperties() const
 {
     return gasProperties_;
 }
diff --git a/src/lagrangian/dieselSpray/sprayThermoTypes/sprayThermoTypes.H b/src/lagrangian/dieselSpray/sprayThermoTypes/sprayThermoTypes.H
deleted file mode 100644
index b52f690b5f88384469aea4b0c402217db41993a7..0000000000000000000000000000000000000000
--- a/src/lagrangian/dieselSpray/sprayThermoTypes/sprayThermoTypes.H
+++ /dev/null
@@ -1,56 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Typedef
-    Foam::specieProperties
-
-Description
-    Foam::specieProperties
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef sprayThermoTypes_H
-#define sprayThermoTypes_H
-
-#include "sutherlandTransport.H"
-#include "specieThermo.H"
-#include "janafThermo.H"
-#include "perfectGas.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    typedef sutherlandTransport<specieThermo<janafThermo<perfectGas> > >
-        specieProperties;
-
-  // typedef sutherlandTransport<specieThermo<hConstThermo<perfectGas> > >
-  //    specieProperties;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/include/createRadiationModel.H b/src/thermophysicalModels/radiation/include/createRadiationModel.H
index babe3c4dbe103dfafed9d392eebf0ea11d51ea4d..d61082c36fdc563e400e4f47bb6ee0050597bc2c 100644
--- a/src/thermophysicalModels/radiation/include/createRadiationModel.H
+++ b/src/thermophysicalModels/radiation/include/createRadiationModel.H
@@ -1,5 +1,5 @@
     Info<< "Creating radiation model\n" << endl;
     autoPtr<radiation::radiationModel> radiation
     (
-        radiation::radiationModel::New(thermo->T())
+        radiation::radiationModel::New(thermo.T())
     );