diff --git a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C
index 89c72685a2..d7995dd405 100644
--- a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C
+++ b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C
@@ -61,7 +61,7 @@ Foam::Analytical<Type>::integrate
 (
     const Type& phi,
     const scalar dt,
-    const Type& alpha,
+    const Type& alphaBeta,
     const scalar beta
 ) const
 {
@@ -69,9 +69,18 @@ Foam::Analytical<Type>::integrate
 
     const scalar expTerm = exp(min(50, -beta*dt));
 
-    retValue.average() =
-        alpha + (phi - alpha)*(1 - expTerm)/(beta*dt + ROOTVSMALL);
-    retValue.value() =  alpha + (phi - alpha)*expTerm;
+    if (beta > ROOTVSMALL)
+    {
+        const Type alpha = alphaBeta/beta;
+        retValue.average() = alpha + (phi - alpha)*(1 - expTerm)/(beta*dt);
+        retValue.value() =  alpha + (phi - alpha)*expTerm;
+    }
+    else
+    {
+        retValue.average() = phi;
+        retValue.value() =  phi;
+    }
+
 
     return retValue;
 }
diff --git a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H
index 91a188dc3a..5eca4f17a8 100644
--- a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H
+++ b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H
@@ -83,7 +83,7 @@ public:
         (
             const Type& phi,
             const scalar dt,
-            const Type& alpha,
+            const Type& alphaBeta,
             const scalar beta
         ) const;
 };
diff --git a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C
index 1ebaa9dc29..ba9cc4c072 100644
--- a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C
+++ b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C
@@ -61,12 +61,12 @@ Foam::Euler<Type>::integrate
 (
     const Type& phi,
     const scalar dt,
-    const Type& alpha,
+    const Type& alphaBeta,
     const scalar beta
 ) const
 {
     typename IntegrationScheme<Type>::integrationResult retValue;
-    retValue.value() = (phi + beta*dt*alpha)/(1.0 + beta*dt);
+    retValue.value() = (phi + alphaBeta*dt)/(1.0 + beta*dt);
     retValue.average() = 0.5*(phi + retValue.value());
 
     return retValue;
diff --git a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H
index de58579a9d..3737d6e902 100644
--- a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H
+++ b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H
@@ -80,7 +80,7 @@ public:
         (
             const Type& phi,
             const scalar dt,
-            const Type& alpha,
+            const Type& alphaBeta,
             const scalar beta
         ) const;
 };
diff --git a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C
index d3cac7f4e0..a8ea549fce 100644
--- a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C
+++ b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C
@@ -62,7 +62,7 @@ Foam::IntegrationScheme<Type>::integrate
 (
     const Type& phi,
     const scalar dt,
-    const Type& alpha,
+    const Type& alphaBeta,
     const scalar beta
 ) const
 {
diff --git a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H
index d7f0d0ef07..941185dcd9 100644
--- a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H
+++ b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H
@@ -185,7 +185,7 @@ public:
         (
             const Type& phi,
             const scalar dt,
-            const Type& alpha,
+            const Type& alphaBeta,
             const scalar beta
         ) const;
 };
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 58390134b3..9a83e69908 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -192,21 +192,20 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     const parcelType& p = static_cast<const parcelType&>(*this);
     const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
     const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
-    forceSuSp Feff = Fcp + Fncp;
-    Feff.Sp() += ROOTVSMALL;
+    const forceSuSp Feff = Fcp + Fncp;
 
 
     // New particle velocity
     //~~~~~~~~~~~~~~~~~~~~~~
 
     // Update velocity - treat as 3-D
-    const vector ap = Uc_ + (Feff.Su() + Su)/Feff.Sp();
+    const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/mass;
     const scalar bp = Feff.Sp()/mass;
 
     Spu = Feff.Sp();
 
     IntegrationScheme<vector>::integrationResult Ures =
-        td.cloud().UIntegrator().integrate(U, dt, ap, bp);
+        td.cloud().UIntegrator().integrate(U, dt, abp, bp);
 
     vector Unew = Ures.value();
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 622ef85fcf..064f945ad1 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -333,7 +333,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     // Integrate to find the new parcel temperature
     IntegrationScheme<scalar>::integrationResult Tres =
-        td.cloud().TIntegrator().integrate(T, dt, ap, bp);
+        td.cloud().TIntegrator().integrate(T, dt, ap*bp, bp);
 
     scalar Tnew = max(Tres.value(), td.cloud().constProps().TMin());