diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H
index 50fc7f575b5e721711a250fa2b7d9d8ddf184fda..ec92091de378b6480529702d917b682c804254ca 100644
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H
@@ -1,22 +1,23 @@
+if (chemistry.chemistry())
 {
     Info<< "Solving chemistry" << endl;
 
-    chemistry.solve
-    (
-        runTime.value() - runTime.deltaTValue(),
-        runTime.deltaTValue()
-    );
+    // update reaction rates
+    chemistry.calculate();
 
     // turbulent time scale
     if (turbulentReaction)
     {
-        DimensionedField<scalar, volMesh> tk =
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
-        DimensionedField<scalar, volMesh> tc =
-            chemistry.tc()().dimensionedInternalField();
+        typedef DimensionedField<scalar, volMesh> dsfType;
 
-        // Chalmers PaSR model
-        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
+        const dimensionedScalar e0("e0", sqr(dimLength)/pow3(dimTime), SMALL);
+
+        const dsfType tk =
+            Cmix*sqrt(turbulence->muEff()/rho/(turbulence->epsilon() + e0));
+
+        const dsfType tc = chemistry.tc()().dimensionedInternalField();
+
+        kappa = tc/(tc + tk);
     }
     else
     {
diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H
index 6ae2712ea42a9490051f2397a0096f77bc6e896c..e7d7c55dbbee62df067024b98aff1d5d23ddcba1 100644
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H
@@ -43,6 +43,11 @@
         }
     }
 
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
+
     // Second part of thermodynamic density update
     thermo.rho() += psi*p;
 
diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H
index f0bcf7597fcf71f1e9b8ee2dbc879200a85fa2cc..e742e9fea78c2d196b07e96196b20ea8a2cf53f1 100644
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H
@@ -1,4 +1,4 @@
-Info<< "Reading chemistry properties\n" << endl;
+// Info<< "Reading chemistry properties\n" << endl;
 
 IOdictionary chemistryProperties
 (
diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H
index 33b7943cadd9e1ecb3433927e4e1a5db4fb543df..e42db723993327c323945353becd4cae0317d6c9 100644
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H
@@ -40,6 +40,9 @@ Description
     );
 
     rhoEqn.solve();
+
+    Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
+        << endl;
 }
 
 // ************************************************************************* //