diff --git a/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.C b/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.C
index eac5aac69a..925468ad61 100644
--- a/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.C
@@ -44,13 +44,17 @@ void Foam::CollidingCloud<CloudType>::setModels()
 
 template<class CloudType>
 template<class TrackData>
-void  Foam::CollidingCloud<CloudType>::moveCollide(TrackData& td)
+void  Foam::CollidingCloud<CloudType>::moveCollide
+(
+    TrackData& td,
+    const scalar deltaT
+)
 {
     td.part() = TrackData::tpVelocityHalfStep;
-    CloudType::move(td,  this->solution().deltaT());
+    CloudType::move(td,  deltaT);
 
     td.part() = TrackData::tpLinearTrack;
-    CloudType::move(td,  this->solution().deltaT());
+    CloudType::move(td,  deltaT);
 
     // td.part() = TrackData::tpRotationalTrack;
     // CloudType::move(td);
@@ -60,7 +64,7 @@ void  Foam::CollidingCloud<CloudType>::moveCollide(TrackData& td)
     this->collision().collide();
 
     td.part() = TrackData::tpVelocityHalfStep;
-    CloudType::move(td,  this->solution().deltaT());
+    CloudType::move(td,  deltaT);
 }
 
 
@@ -211,14 +215,14 @@ void  Foam::CollidingCloud<CloudType>::motion(TrackData& td)
 
         while(!(++moveCollideSubCycle).end())
         {
-            moveCollide(td);
+            moveCollide(td, this->db().time().deltaTValue());
         }
 
         moveCollideSubCycle.endSubCycle();
     }
     else
     {
-        moveCollide(td);
+        moveCollide(td, this->db().time().deltaTValue());
     }
 }
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.H b/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.H
index 9ca29cfd56..68fbc0000e 100644
--- a/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/CollidingCloud/CollidingCloud.H
@@ -108,7 +108,7 @@ protected:
 
             //- Move-collide particles
             template<class TrackData>
-            void moveCollide(TrackData& td);
+            void moveCollide(TrackData& td, const scalar deltaT);
 
             //- Reset state of cloud
             void cloudReset(CollidingCloud<CloudType>& c);
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index d2aa164294..58390134b3 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -192,14 +192,15 @@ 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);
-    const forceSuSp Feff = Fcp + Fncp;
+    forceSuSp Feff = Fcp + Fncp;
+    Feff.Sp() += ROOTVSMALL;
 
 
     // New particle velocity
     //~~~~~~~~~~~~~~~~~~~~~~
 
     // Update velocity - treat as 3-D
-    const vector ap = Uc_ + (Feff.Su() + Su)/(Feff.Sp() + ROOTVSMALL);
+    const vector ap = Uc_ + (Feff.Su() + Su)/Feff.Sp();
     const scalar bp = Feff.Sp()/mass;
 
     Spu = Feff.Sp();
@@ -295,6 +296,7 @@ bool Foam::KinematicParcel<ParcelType>::move
     const polyMesh& mesh = td.cloud().pMesh();
     const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
     const scalarField& V = mesh.cellVolumes();
+    const scalar maxCo = td.cloud().solution().maxCo();
 
     scalar tEnd = (1.0 - p.stepFraction())*trackTime;
     const scalar dtMax = tEnd;
@@ -315,7 +317,6 @@ bool Foam::KinematicParcel<ParcelType>::move
         if (p.active() && magU > ROOTVSMALL)
         {
             const scalar d = dt*magU;
-            const scalar maxCo = td.cloud().solution().maxCo();
             const scalar dCorr = min(d, maxCo*cbrt(V[cellI]));
             dt *=
                 dCorr/d