diff --git a/applications/solvers/combustion/reactingFoam/chemistry.H b/applications/solvers/combustion/reactingFoam/chemistry.H
index 3a36aaac520a9bf63fd20c6b0373ae7e4751a61f..99f418af6f10ac69b44d08fc2271e2dd79c73fc4 100644
--- a/applications/solvers/combustion/reactingFoam/chemistry.H
+++ b/applications/solvers/combustion/reactingFoam/chemistry.H
@@ -1,3 +1,4 @@
+if (chemistry.chemistry())
 {
     Info<< "Solving chemistry" << endl;
 
@@ -10,17 +11,29 @@
     // turbulent time scale
     if (turbulentReaction)
     {
-        volScalarField tk
-        (
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
-        );
-        volScalarField tc
-        (
-            chemistry.tc()
-        );
+        tmp<volScalarField> tepsilon(turbulence->epsilon());
+        const volScalarField& epsilon = tepsilon();
+        tmp<volScalarField> tmuEff(turbulence->muEff());
+        const volScalarField& muEff = tmuEff();
+        tmp<volScalarField> ttc(chemistry.tc());
+        const volScalarField& tc = ttc();
 
-        // Chalmers PaSR model
-        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
+        forAll(epsilon, i)
+        {
+            if (epsilon[i] > 0)
+            {
+                // Chalmers PaSR model
+                scalar tk = Cmix.value()*Foam::sqrt(muEff[i]/rho[i]/epsilon[i]);
+                kappa[i] =
+                    (runTime.deltaTValue() + tc[i])
+                   /(runTime.deltaTValue() + tc[i] + tk);
+            }
+            else
+            {
+                // Return to laminar combustion
+                kappa[i] = 1.0;
+            }
+        }
     }
     else
     {
diff --git a/applications/solvers/combustion/rhoReactingFoam/Make/options b/applications/solvers/combustion/rhoReactingFoam/Make/options
index 6386af8c557fcea5ace76bfe7c203c36ab85cd6e..d6306816fd8174d34268711983ddb53cdb4eb040 100644
--- a/applications/solvers/combustion/rhoReactingFoam/Make/options
+++ b/applications/solvers/combustion/rhoReactingFoam/Make/options
@@ -5,7 +5,9 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
 
 EXE_LIBS = \
     -lcompressibleTurbulenceModel \
diff --git a/applications/solvers/combustion/rhoReactingFoam/chemistry.H b/applications/solvers/combustion/rhoReactingFoam/chemistry.H
deleted file mode 100644
index 841b434a7176d259a77f9d6f756bd1834af773a5..0000000000000000000000000000000000000000
--- a/applications/solvers/combustion/rhoReactingFoam/chemistry.H
+++ /dev/null
@@ -1,28 +0,0 @@
-{
-    Info<< "Solving chemistry" << endl;
-
-    chemistry.solve
-    (
-        runTime.value() - runTime.deltaTValue(),
-        runTime.deltaTValue()
-    );
-
-    // turbulent time scale
-    if (turbulentReaction)
-    {
-        volScalarField tk
-        (
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
-        );
-        volScalarField tc(chemistry.tc());
-
-        // Chalmers PaSR model
-        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
-    }
-    else
-    {
-        kappa = 1.0;
-    }
-
-    chemistrySh = kappa*chemistry.Sh()();
-}
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/Make/options b/applications/solvers/lagrangian/coalChemistryFoam/Make/options
index ab0bd2c08a582fa8af1eb5acd781c2773d696e5f..75bbcab490ee4ed3e85e4c2fbcc864e300f03455 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/Make/options
+++ b/applications/solvers/lagrangian/coalChemistryFoam/Make/options
@@ -19,7 +19,9 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/ODE/lnInclude
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
 
 EXE_LIBS = \
     -lfiniteVolume \
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/chemistry.H b/applications/solvers/lagrangian/coalChemistryFoam/chemistry.H
deleted file mode 100644
index a1a5cb39ba27fde8c7492d6b43bc9a28a2141188..0000000000000000000000000000000000000000
--- a/applications/solvers/lagrangian/coalChemistryFoam/chemistry.H
+++ /dev/null
@@ -1,32 +0,0 @@
-if (chemistry.chemistry())
-{
-    Info<< "Solving chemistry" << endl;
-
-    chemistry.solve
-    (
-        runTime.value() - runTime.deltaTValue(),
-        runTime.deltaTValue()
-    );
-
-    // 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;
-    }
-
-    chemistrySh = kappa*chemistry.Sh()();
-}
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options
index 1ff1aaa6129259368b371693f66f75fa3b1bac4f..386e25dd477398c5d64f2f6324c866241e3cc13c 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options
@@ -19,7 +19,9 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
 
 EXE_LIBS = \
     -lfiniteVolume \
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H
deleted file mode 100644
index 1690487b45db2e49964b43a9219cdbfa1d08853f..0000000000000000000000000000000000000000
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H
+++ /dev/null
@@ -1,31 +0,0 @@
-{
-    Info<< "Solving chemistry" << endl;
-
-    chemistry.solve
-    (
-        runTime.value() - runTime.deltaTValue(),
-        runTime.deltaTValue()
-    );
-
-    // 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;
-    }
-
-    chemistrySh = kappa*chemistry.Sh()();
-}
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options
index 1efc598b5b1dedd44f98783de44fab8e39f650a7..cd1c25ffea95df023e9392b53723c0df5173e3f6 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options
@@ -18,7 +18,9 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/ODE/lnInclude
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
 
 EXE_LIBS = \
     -lfiniteVolume \
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H
deleted file mode 100644
index ce36fa0d1b5b6b472958ed9c38a0010c4f710888..0000000000000000000000000000000000000000
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H
+++ /dev/null
@@ -1,32 +0,0 @@
-if (chemistry.chemistry())
-{
-    Info << "Solving chemistry" << endl;
-
-    chemistry.solve
-    (
-        runTime.value() - runTime.deltaTValue(),
-        runTime.deltaTValue()
-    );
-
-    // 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;
-    }
-
-    chemistrySh = kappa*chemistry.Sh()();
-}
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFoam/Make/options
index fb54b88e26984daa386646fdb38b5cc2b968e25a..1380c457e269c3ac598407c7b200c72032e9170c 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/Make/options
@@ -18,7 +18,9 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
 
 EXE_LIBS = \
     -lfiniteVolume \
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/reactingParcelFoam/chemistry.H
deleted file mode 100644
index 1690487b45db2e49964b43a9219cdbfa1d08853f..0000000000000000000000000000000000000000
--- a/applications/solvers/lagrangian/reactingParcelFoam/chemistry.H
+++ /dev/null
@@ -1,31 +0,0 @@
-{
-    Info<< "Solving chemistry" << endl;
-
-    chemistry.solve
-    (
-        runTime.value() - runTime.deltaTValue(),
-        runTime.deltaTValue()
-    );
-
-    // 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;
-    }
-
-    chemistrySh = kappa*chemistry.Sh()();
-}