diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
index 321c8235..109d2212 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
@@ -792,6 +792,10 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
     scalarField c(nSpecie_);
     scalarField c0(nSpecie_);
 
+    const bool neFound = this->mesh().template foundObject<volScalarField>("ne");
+
+    if (!neFound)
+    {
     forAll(rho, celli)
     {
         const scalar rhoi = rho[celli];
@@ -822,6 +826,57 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
             RR_[i][celli] = (c[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
         }
     }
+    }
+    else
+    {
+    const scalarField& Te = this->mesh().template lookupObject<volScalarField>("Te").internalField();
+    const scalarField& ne = this->mesh().template lookupObject<volScalarField>("ne").internalField();
+
+    forAll(rho, celli)
+    {
+        const scalar rhoi = rho[celli];
+        scalar pi = p[celli];
+        scalar Ti = T[celli];
+        scalar nei = ne[celli];
+        scalar Tei = Te[celli];
+
+        nei = Tei;
+        Tei = nei;
+
+        for (label i=0; i<nReaction_; i++)
+        {
+            reactions_[i].species().contains("E-");
+        }
+
+        for (label i=0; i<nSpecie_; i++)
+        {
+            {
+                c[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
+                c0[i] = c[i];
+            }
+        }
+
+        // Initialise time progress
+        scalar timeLeft = deltaT[celli];
+
+        // Calculate the chemical source terms
+        while (timeLeft > SMALL)
+        {
+            scalar dt = timeLeft;
+            this->solve(c, Ti, pi, dt, this->deltaTChem_[celli]);
+            timeLeft -= dt;
+        }
+
+        deltaTMin = min(this->deltaTChem_[celli], deltaTMin);
+
+        for (label i=0; i<nSpecie_; i++)
+        {
+            {
+                RR_[i][celli] = (c[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
+            }
+        }
+    }
+    }
 
     return deltaTMin;
 }