diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options
index c14a8799f..043be5d13 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options
@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I. \
     -I../buoyantBoussinesqSimpleFoam \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
index d522323fd..e082648b7 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
@@ -56,19 +56,18 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "postProcess.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-
-    pimpleControl pimple(mesh);
-
+    #include "createControl.H"
     #include "createFields.H"
-    #include "createIncompressibleRadiationModel.H"
     #include "createFvOptions.H"
-    #include "initContinuityErrs.H"
     #include "createTimeControls.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
+    #include "initContinuityErrs.H"
 
     turbulence->validate();
 
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
index cc2c5f0b4..def60953c 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
@@ -122,3 +122,4 @@ if (p_rgh.needReference())
 mesh.setFluxRequired(p_rgh.name());
 
 #include "createMRF.H"
+#include "createIncompressibleRadiationModel.H"
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
index 6c4955dfb..de844ad96 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
@@ -56,14 +56,13 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "postProcess.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-
-    simpleControl simple(mesh);
-
+    #include "createControl.H"
     #include "createFields.H"
-    #include "createIncompressibleRadiationModel.H"
     #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
index 583cbaf2d..043440763 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
@@ -122,3 +122,4 @@ if (p_rgh.needReference())
 mesh.setFluxRequired(p_rgh.name());
 
 #include "createMRF.H"
+#include "createIncompressibleRadiationModel.H"
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index 56ebdc5c8..d960e053d 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -44,15 +44,15 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "postProcess.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-
-    pimpleControl pimple(mesh);
-
+    #include "createControl.H"
     #include "createFields.H"
+    #include "createFieldRefs.H"
     #include "createFvOptions.H"
-    #include "createRadiationModel.H"
     #include "initContinuityErrs.H"
     #include "createTimeControls.H"
     #include "compressibleCourantNo.H"
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFieldRefs.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFieldRefs.H
new file mode 100644
index 000000000..502b3b423
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFieldRefs.H
@@ -0,0 +1 @@
+const volScalarField& psi = thermo.psi();
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index b3c45fcaa..ff835c8a7 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@@ -18,8 +18,6 @@ volScalarField rho
 );
 
 volScalarField& p = thermo.p();
-const volScalarField& psi = thermo.psi();
-
 
 Info<< "Reading field U\n" << endl;
 volVectorField U
@@ -92,3 +90,4 @@ Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
 #include "createMRF.H"
+#include "createRadiationModel.H"
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
index cb4d67a0c..e9add9cd0 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
@@ -41,15 +41,15 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "postProcess.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-
-    simpleControl simple(mesh);
-
+    #include "createControl.H"
     #include "createFields.H"
+    #include "createFieldRefs.H"
     #include "createFvOptions.H"
-    #include "createRadiationModel.H"
     #include "initContinuityErrs.H"
 
     turbulence->validate();
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFieldRefs.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFieldRefs.H
new file mode 100644
index 000000000..502b3b423
--- /dev/null
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFieldRefs.H
@@ -0,0 +1 @@
+const volScalarField& psi = thermo.psi();
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index ca107c092..b9f275407 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -18,7 +18,6 @@ volScalarField rho
 );
 
 volScalarField& p = thermo.p();
-const volScalarField& psi = thermo.psi();
 
 Info<< "Reading field U\n" << endl;
 volVectorField U
@@ -88,3 +87,4 @@ dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 dimensionedScalar totalVolume = sum(mesh.V());
 
 #include "createMRF.H"
+#include "createRadiationModel.H"
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
index 2b6f5e7e7..fd53d162e 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
@@ -1,5 +1,6 @@
 
 EXE_INC = \
+    -I. \
     -I./fluid \
     -I./solid \
     -I./porousFluid \
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
index 1ad39c9d0..7a396ee38 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -50,22 +50,17 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #define NO_CONTROL
+    #define CREATE_MESH createMeshesPostProcess.H
+    #include "postProcess.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
-
-    regionProperties rp(runTime);
-
-    #include "createFluidMeshes.H"
-    #include "createSolidMeshes.H"
-
-    #include "createFluidFields.H"
-    #include "createSolidFields.H"
-
+    #include "createMeshes.H"
+    #include "createFields.H"
     #include "initContinuityErrs.H"
     #include "createTimeControls.H"
     #include "readSolidTimeControls.H"
-
-
     #include "compressibleMultiRegionCourantNo.H"
     #include "solidRegionDiffusionNo.H"
     #include "setInitialMultiRegionDeltaT.H"
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options
index 252ecae02..78fc5be02 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options
@@ -1,4 +1,6 @@
 EXE_INC = \
+    -I. \
+    -I.. \
     -Ifluid \
     -Isolid \
     -I../solid \
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
index 6b4811a02..9992f8f70 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
@@ -43,17 +43,14 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #define NO_CONTROL
+    #define CREATE_MESH createMeshesPostProcess.H
+    #include "postProcess.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
-
-    regionProperties rp(runTime);
-
-    #include "createFluidMeshes.H"
-    #include "createSolidMeshes.H"
-
-    #include "createFluidFields.H"
-    #include "createSolidFields.H"
-
+    #include "createMeshes.H"
+    #include "createFields.H"
     #include "initContinuityErrs.H"
 
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/createFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/createFields.H
new file mode 100644
index 000000000..30f445230
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/createFields.H
@@ -0,0 +1,2 @@
+#include "createFluidFields.H"
+#include "createSolidFields.H"
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/createMeshes.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/createMeshes.H
new file mode 100644
index 000000000..97aed2177
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/createMeshes.H
@@ -0,0 +1,4 @@
+regionProperties rp(runTime);
+
+#include "createFluidMeshes.H"
+#include "createSolidMeshes.H"
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/createMeshesPostProcess.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/createMeshesPostProcess.H
new file mode 100644
index 000000000..61d94eb5c
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/createMeshesPostProcess.H
@@ -0,0 +1,2 @@
+#include "createMeshes.H"
+fvMesh& mesh = fluidRegions[0];
diff --git a/applications/solvers/heatTransfer/thermoFoam/createFields.H b/applications/solvers/heatTransfer/thermoFoam/createFields.H
index ffa6c3f3b..3a217895f 100644
--- a/applications/solvers/heatTransfer/thermoFoam/createFields.H
+++ b/applications/solvers/heatTransfer/thermoFoam/createFields.H
@@ -1,54 +1,56 @@
-    Info<< "Reading thermophysical properties\n" << endl;
+Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<rhoThermo> pThermo(rhoThermo::New(mesh));
-    rhoThermo& thermo = pThermo();
-    thermo.validate(args.executable(), "h", "e");
+autoPtr<rhoThermo> pThermo(rhoThermo::New(mesh));
+rhoThermo& thermo = pThermo();
+thermo.validate(args.executable(), "h", "e");
 
-    volScalarField rho
+volScalarField rho
+(
+    IOobject
     (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        thermo.rho()
-    );
-
-    volScalarField& p = thermo.p();
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    #include "compressibleCreatePhi.H"
-
-    #include "setAlphaEff.H"
-
-    Info<< "Creating field dpdt\n" << endl;
-    volScalarField dpdt
-    (
-        IOobject
-        (
-            "dpdt",
-            runTime.timeName(),
-            mesh
-        ),
+        "rho",
+        runTime.timeName(),
         mesh,
-        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
-    );
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    thermo.rho()
+);
 
-    Info<< "Creating field kinetic energy K\n" << endl;
-    volScalarField K("K", 0.5*magSqr(U));
+volScalarField& p = thermo.p();
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+#include "setAlphaEff.H"
+
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+#include "createRadiationModel.H"
diff --git a/applications/solvers/heatTransfer/thermoFoam/setAlphaEff.H b/applications/solvers/heatTransfer/thermoFoam/setAlphaEff.H
index fff960526..73a157fc1 100644
--- a/applications/solvers/heatTransfer/thermoFoam/setAlphaEff.H
+++ b/applications/solvers/heatTransfer/thermoFoam/setAlphaEff.H
@@ -1,49 +1,47 @@
-    Info<< "Creating turbulence model\n" << endl;
-    tmp<volScalarField> talphaEff;
+Info<< "Creating turbulence model\n" << endl;
+tmp<volScalarField> talphaEff;
 
-    IOobject turbulencePropertiesHeader
+IOobject turbulencePropertiesHeader
+(
+    "turbulenceProperties",
+    runTime.constant(),
+    mesh,
+    IOobject::MUST_READ_IF_MODIFIED,
+    IOobject::NO_WRITE,
+    false
+);
+
+if (turbulencePropertiesHeader.headerOk())
+{
+    autoPtr<compressible::turbulenceModel> turbulence
     (
-        "turbulenceProperties",
-        runTime.constant(),
-        mesh,
-        IOobject::MUST_READ_IF_MODIFIED,
-        IOobject::NO_WRITE,
-        false
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
     );
 
-    if (turbulencePropertiesHeader.headerOk())
-    {
-        autoPtr<compressible::turbulenceModel> turbulence
+    talphaEff = turbulence->alphaEff();
+}
+else
+{
+    talphaEff = tmp<volScalarField>
+    (
+        new volScalarField
         (
-            compressible::turbulenceModel::New
+            IOobject
             (
-                rho,
-                U,
-                phi,
-                thermo
-            )
-        );
-
-        talphaEff = turbulence->alphaEff();
-    }
-    else
-    {
-        talphaEff = tmp<volScalarField>
-        (
-            new volScalarField
-            (
-                IOobject
-                (
-                    "alphaEff",
-                    runTime.timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
+                "alphaEff",
+                runTime.timeName(),
                 mesh,
-                dimensionedScalar("0", dimMass/dimLength/dimTime, 0.0)
-            )
-        );
-    }
-
-    const volScalarField& alphaEff = talphaEff();
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("0", dimMass/dimLength/dimTime, 0.0)
+        )
+    );
+}
diff --git a/applications/solvers/heatTransfer/thermoFoam/thermoFoam.C b/applications/solvers/heatTransfer/thermoFoam/thermoFoam.C
index a3687213d..577091232 100644
--- a/applications/solvers/heatTransfer/thermoFoam/thermoFoam.C
+++ b/applications/solvers/heatTransfer/thermoFoam/thermoFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,13 +43,16 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #include "setRootCase.H"
+    #define NO_CONTROL
+    #include "postProcess.H"
 
+    #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
     #include "createFvOptions.H"
-    #include "createRadiationModel.H"
+
+    const volScalarField& alphaEff = talphaEff();
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //