diff --git a/applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C b/applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C
index c4da39ff..dd281b5f 100644
--- a/applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C
+++ b/applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C
@@ -57,6 +57,8 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "postProcess.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
@@ -64,9 +66,10 @@ int main(int argc, char *argv[])
     #include "createTimeControls.H"
     #include "createRDeltaT.H"
 
-    #include "createFields.H"
-    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createFvOptions.H"
 
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
diff --git a/applications/solvers/cfdemSolverRhoPimpleChem/createFields.H b/applications/solvers/cfdemSolverRhoPimpleChem/createFields.H
index 4cda15ee..d5daf875 100644
--- a/applications/solvers/cfdemSolverRhoPimpleChem/createFields.H
+++ b/applications/solvers/cfdemSolverRhoPimpleChem/createFields.H
@@ -24,7 +24,7 @@
     volScalarField W(thermo.W());
 #endif
 
-    bool propagateInertSpecie = true;
+    Switch propagateInertSpecie(true);
 
     const word inertSpecie(thermo.lookup("inertSpecie"));
 
@@ -40,9 +40,9 @@
             << exit(FatalError);
     }
 
+    Info<< "inert will be bounded in [" << inertLowerBound << "," << inertUpperBound << "]" << endl;
+
     volScalarField& p = thermo.p();
-    const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo.psi();
 
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
diff --git a/applications/solvers/rcfdemSolverBase/createFields.H b/applications/solvers/rcfdemSolverBase/createFields.H
index 7088bab5..16b75ff3 100644
--- a/applications/solvers/rcfdemSolverBase/createFields.H
+++ b/applications/solvers/rcfdemSolverBase/createFields.H
@@ -44,7 +44,7 @@
         dimensionedVector("URec", dimensionSet(0, 1, -1, 0, 0), vector::zero)
     );
 
-    bool updateURec = false;
+    Switch updateURec(false);
     if (URec.headerOk())
     {
         updateURec = true;
@@ -65,7 +65,7 @@
         dimensionedScalar("voidfractionRec", dimensionSet(0, 0, 0, 0, 0), 1.0)
     );
 
-    bool updateVoidfractionRec = false;
+    Switch updateVoidfractionRec(false);
     if (voidfractionRec.headerOk())
     {
         updateVoidfractionRec = true;
@@ -86,7 +86,7 @@
         dimensionedVector("URec", dimensionSet(0, 1, -1, 0, 0), vector::zero)
     );
 
-    bool updateUsRec = false;
+    Switch updateUsRec(false);
     if (UsRec.headerOk())
     {
         updateUsRec = true;
@@ -140,7 +140,7 @@
         linearInterpolate(URec*voidfractionRec) & mesh.Sf()
     );
 
-    bool updatePhiRec = false;
+    Switch updatePhiRec(false);
     if (phiRec.headerOk())
     {
         updatePhiRec = true;
diff --git a/applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C b/applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C
index 0a4f5efb..b2d6993a 100644
--- a/applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C
+++ b/applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C
@@ -58,7 +58,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createFields.H"
-    #include "createFvOptions.H"
 
     #include "readGravitationalAcceleration.H"
 
diff --git a/applications/solvers/rcfdemSolverForcedTracers/rcfdemSolverForcedTracers.C b/applications/solvers/rcfdemSolverForcedTracers/rcfdemSolverForcedTracers.C
index 9039e038..80a55195 100644
--- a/applications/solvers/rcfdemSolverForcedTracers/rcfdemSolverForcedTracers.C
+++ b/applications/solvers/rcfdemSolverForcedTracers/rcfdemSolverForcedTracers.C
@@ -51,7 +51,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createFields.H"
-    #include "createFvOptions.H"
 
     cfdemCloud particleCloud(mesh);
     recBase recurrenceBase(mesh);
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H
index 164799e8..4a65a741 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H
@@ -168,7 +168,7 @@ Info<< "Reading thermophysical properties\n" << endl;
          linearInterpolate(rho*U*voidfraction) & mesh.Sf()
     );
 
-    bool transientEEqn(pimple.dict().lookupOrDefault<bool>("transientEEqn",false));
+    Switch transientEEqn(pimple.dict().lookupOrDefault<bool>("transientEEqn",false));
 
     dimensionedScalar rhoMax
     (
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C b/applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C
index 6add71a1..01098c87 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C
@@ -63,7 +63,6 @@ int main(int argc, char *argv[])
     #include "createControl.H"
     #include "createTimeControls.H"
     #include "createRDeltaT.H"
-    #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
     #include "createFvOptions.H"
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H
index 1556801e..e7640d48 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H
@@ -24,7 +24,7 @@ Info<< "Reading thermophysical properties\n" << endl;
     volScalarField W(thermo.W());
 #endif
 
-    bool propagateInertSpecie(thermo.lookupOrDefault<bool>("propagateInertSpecie",true));
+    Switch propagateInertSpecie(thermo.lookupOrDefault<bool>("propagateInertSpecie",true));
 
     const word inertSpecie(thermo.lookupOrDefault<word>("inertSpecie","none"));
 
@@ -40,9 +40,9 @@ Info<< "Reading thermophysical properties\n" << endl;
             << exit(FatalError);
     }
 
+    Info<< "inert will be bounded in [" << inertLowerBound << "," << inertUpperBound << "]" << endl;
+
     volScalarField& p = thermo.p();
-    const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo.psi();
 
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
@@ -212,7 +212,7 @@ Info<< "Reading thermophysical properties\n" << endl;
          linearInterpolate(rho*U*voidfraction) & mesh.Sf()
     );
 
-    bool transientEEqn(pimple.dict().lookupOrDefault<bool>("transientEEqn",false));
+    Switch transientEEqn(pimple.dict().lookupOrDefault<bool>("transientEEqn",false));
 
     dimensionedScalar rhoMax
     (
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C
index 682fc5b0..ef6420eb 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C
@@ -73,7 +73,7 @@ int main(int argc, char *argv[])
 
     #include "initContinuityErrs.H"
     #include "createFields.H"
-    //#include "createFieldRefs.H"
+    #include "createFieldRefs.H"
     #include "createFvOptions.H"
 
     // create cfdemCloud
diff --git a/applications/solvers/rctfSpeciesTransport/createFields.H b/applications/solvers/rctfSpeciesTransport/createFields.H
index 3b6ef0e5..71de9463 100755
--- a/applications/solvers/rctfSpeciesTransport/createFields.H
+++ b/applications/solvers/rctfSpeciesTransport/createFields.H
@@ -255,7 +255,6 @@ volScalarField alphat
     
     fvScalarMatrix CEqn(C, dimless*dimVolume/(dimTime));
 
-    scalar relaxCoeff(0.0);
 
  Info<< "reading clockProperties\n" << endl;
 
diff --git a/applications/solvers/rctfSpeciesTransport/rctfSpeciesTransport.C b/applications/solvers/rctfSpeciesTransport/rctfSpeciesTransport.C
index 061a0fb2..3cd080ad 100755
--- a/applications/solvers/rctfSpeciesTransport/rctfSpeciesTransport.C
+++ b/applications/solvers/rctfSpeciesTransport/rctfSpeciesTransport.C
@@ -58,10 +58,11 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createControl.H"
+    #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFvOptions.H"
-    scalar cumulativeContErr = 0;
-  
+    scalar relaxCoeff(0.0);
+
     //create recBases according to a list of recProperties
     #include "createRecBase.H"
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/recSolverTurbTransport/createFields.H b/applications/solvers/recSolverTurbTransport/createFields.H
index 3b9b3af1..8d520e5c 100755
--- a/applications/solvers/recSolverTurbTransport/createFields.H
+++ b/applications/solvers/recSolverTurbTransport/createFields.H
@@ -147,7 +147,6 @@ surfaceScalarField phiRec
     
     fvScalarMatrix TEqn(T, dimless*dimVolume/(dimTime));
 
-    scalar relaxCoeff(0.0);
 
  Info<< "reading clockProperties\n" << endl;
 
diff --git a/applications/solvers/recSolverTurbTransport/recSolverTurbTransport.C b/applications/solvers/recSolverTurbTransport/recSolverTurbTransport.C
index a93997a4..92c2e9c6 100755
--- a/applications/solvers/recSolverTurbTransport/recSolverTurbTransport.C
+++ b/applications/solvers/recSolverTurbTransport/recSolverTurbTransport.C
@@ -55,6 +55,7 @@ int main(int argc, char *argv[])
     #include "createControl.H"
     #include "createFields.H"
     #include "createFvOptions.H"
+    scalar relaxCoeff(0.0);
 
     recBase recurrenceBase(mesh);
 
diff --git a/applications/solvers/rtfmSolverSpecies/rtfmSolverSpecies.C b/applications/solvers/rtfmSolverSpecies/rtfmSolverSpecies.C
index 0096009a..12f11b00 100644
--- a/applications/solvers/rtfmSolverSpecies/rtfmSolverSpecies.C
+++ b/applications/solvers/rtfmSolverSpecies/rtfmSolverSpecies.C
@@ -58,6 +58,7 @@ int main(int argc, char *argv[])
     #include "createControl.H"
     #include "createFields.H"
     #include "createFvOptions.H"
+    scalar relaxCoeff(0.0);
 
     cfdemCloudRec<cfdemCloud> particleCloud(mesh);
     recBase recurrenceBase(mesh);