diff --git a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
index 5caaedf2..722f6721 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
@@ -17,22 +17,38 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     composition.calculateDiffusivities(p, T);
 
+    EnTd = En.internalField();
+    EnTd *= EnToTableUnit;
+
+    Te.internalField() = TeOfEn.value(EnTd) * TeFac;
     forAll(rho, celli)
     {
-        const scalar Eni = En[celli] * EnFac;
-
-        De[celli] = DeN.value(Eni) * DeNFac;
-        mue[celli] = mueN.value(Eni) * mueNFac;
-        Te[celli] = TeOfEn.value(Eni) * TeFac;
+        Te[celli] = max(Te[celli], T[celli]);
     }
-    De.correctBoundaryConditions();
-    mue.correctBoundaryConditions();
     Te.correctBoundaryConditions();
 
+    if (mobility_f_of_Te)
+    {
+        EnTd = Te.internalField();
+        EnTd *= TeToTableUnit;
+    }
+
+    mue.internalField() = mueN.value(EnTd) * mueNFac;
+
+    if (calculateDe)
+    {
+        De = mue * Te * (kB / eCharge);
+    }
+    else
+    {
+        De.internalField() = DeN.value(EnTd) * DeNFac;
+    }
+
+    mue.correctBoundaryConditions();
+    De.correctBoundaryConditions();
 
     q = linearInterpolate(U) & mesh.Sf();
 
-
     forAll(Y, i)
     {
         volScalarField& Yi = Y[i];
@@ -66,7 +82,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
             }
 
             tmp<fvScalarMatrix> electronR(
-                    new fvScalarMatrix(ne, dimless/dimTime));
+                    new fvScalarMatrix(ne,
+                        ne.dimensions()*dimVol/dimTime));
             electronR->source() = reaction->R(Yi)->source();
 
             fvScalarMatrix neEqn
@@ -88,6 +105,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
             fvOptions.correct(ne);
 
+            ne.writeMinMax(Info);
+
             ne.max(0.0);
         }
         else if (Y[i].name() != inertSpecie)
diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
index 794066fb..0ace61f9 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/createFields.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
@@ -17,10 +17,20 @@ dimensionedScalar epsilon0
     physicalProperties.lookup("epsilon0")
 );
 
+Switch mobility_f_of_Te = physicalProperties.lookupOrDefault("mobility_f_of_Te", false);
+
+Switch calculateDe = physicalProperties.lookupOrDefault("calculateDe", false);
+
 // Convert E/n in SI unit to table unit.
 // Default V m^2 => Td
-scalar EnFac (
-        physicalProperties.lookupOrDefault("EnFac", 1.0e21)
+scalar TeToTableUnit (
+        physicalProperties.lookupOrDefault("TeToTableUnit", 1.0)
+        );
+
+// Convert E/n in SI unit to table unit.
+// Default V m^2 => Td
+scalar EnToTableUnit (
+        physicalProperties.lookupOrDefault("EnToTableUnit", 1.0e21)
         );
 
 
@@ -277,6 +287,7 @@ volScalarField ng("ng", p / R / T * NA);
 
 Info<< "Creating field reduced electric field\n" << endl;
 volScalarField En ("En", mag(E) / (ng));
+scalarField EnTd(En.internalField() * EnToTableUnit);
 
 Info<< "Creating field electron mobility\n" << endl;
 volScalarField mue