diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index c51587c174..77dab8ad9f 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -48,9 +48,9 @@ void Foam::KinematicParcel<ParcelType>::setCellValues
 {
     tetIndices tetIs = this->currentTetIndices();
 
-    rhoc_ = td.rhoInterp().interpolate(this->coordinates(), tetIs);
+    td.rhoc() = td.rhoInterp().interpolate(this->coordinates(), tetIs);
 
-    if (rhoc_ < cloud.constProps().rhoMin())
+    if (td.rhoc() < cloud.constProps().rhoMin())
     {
         if (debug)
         {
@@ -59,20 +59,20 @@ void Foam::KinematicParcel<ParcelType>::setCellValues
                 << " to " << cloud.constProps().rhoMin() <<  nl << endl;
         }
 
-        rhoc_ = cloud.constProps().rhoMin();
+        td.rhoc() = cloud.constProps().rhoMin();
     }
 
-    Uc_ = td.UInterp().interpolate(this->coordinates(), tetIs);
+    td.Uc() = td.UInterp().interpolate(this->coordinates(), tetIs);
 
-    muc_ = td.muInterp().interpolate(this->coordinates(), tetIs);
+    td.muc() = td.muInterp().interpolate(this->coordinates(), tetIs);
 
     // Apply dispersion components to carrier phase velocity
-    Uc_ = cloud.dispersion().update
+    td.Uc() = cloud.dispersion().update
     (
         dt,
         this->cell(),
         U_,
-        Uc_,
+        td.Uc(),
         UTurb_,
         tTurb_
     );
@@ -88,7 +88,7 @@ void Foam::KinematicParcel<ParcelType>::cellValueSourceCorrection
     const scalar dt
 )
 {
-    Uc_ += cloud.UTrans()[this->cell()]/massCell(this->cell());
+    td.Uc() += cloud.UTrans()[this->cell()]/massCell(td);
 }
 
 
@@ -107,7 +107,7 @@ void Foam::KinematicParcel<ParcelType>::calc
     const scalar mass0 = mass();
 
     // Reynolds number
-    const scalar Re = this->Re(U_, d_, rhoc_, muc_);
+    const scalar Re = this->Re(td);
 
 
     // Sources
@@ -127,7 +127,8 @@ void Foam::KinematicParcel<ParcelType>::calc
     // ~~~~~~
 
     // Calculate new particle velocity
-    this->U_ = calcVelocity(cloud, td, dt, Re, muc_, mass0, Su, dUTrans, Spu);
+    this->U_ =
+        calcVelocity(cloud, td, dt, Re, td.muc(), mass0, Su, dUTrans, Spu);
 
 
     // Accumulate carrier phase source terms
@@ -158,24 +159,25 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     scalar& Spu
 ) const
 {
-    typedef typename TrackCloudType::parcelType parcelType;
-    typedef typename TrackCloudType::forceType forceType;
+    const typename TrackCloudType::parcelType& p =
+        static_cast<const typename TrackCloudType::parcelType&>(*this);
+    typename TrackCloudType::parcelType::trackingData& ttd =
+        static_cast<typename TrackCloudType::parcelType::trackingData&>(td);
 
-    const forceType& forces = cloud.forces();
+    const typename TrackCloudType::forceType& forces = cloud.forces();
 
     // Momentum source due to particle forces
-    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 Fcp = forces.calcCoupled(p, ttd, dt, mass, Re, mu);
+    const forceSuSp Fncp = forces.calcNonCoupled(p, ttd, dt, mass, Re, mu);
     const forceSuSp Feff = Fcp + Fncp;
-    const scalar massEff = forces.massEff(p, mass);
+    const scalar massEff = forces.massEff(p, ttd, mass);
 
 
     // New particle velocity
     //~~~~~~~~~~~~~~~~~~~~~~
 
     // Update velocity - treat as 3-D
-    const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/massEff;
+    const vector abp = (Feff.Sp()*td.Uc() + (Feff.Su() + Su))/massEff;
     const scalar bp = Feff.Sp()/massEff;
 
     Spu = dt*Feff.Sp();
@@ -186,7 +188,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     vector Unew = Ures.value();
 
     // note: Feff.Sp() and Fc.Sp() must be the same
-    dUTrans += dt*(Feff.Sp()*(Ures.average() - Uc_) - Fcp.Su());
+    dUTrans += dt*(Feff.Sp()*(Ures.average() - td.Uc()) - Fcp.Su());
 
     // Apply correction to velocity and dUTrans for reduced-D cases
     const polyMesh& mesh = cloud.pMesh();
@@ -215,10 +217,7 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(p.rho_),
     age_(p.age_),
     tTurb_(p.tTurb_),
-    UTurb_(p.UTurb_),
-    rhoc_(p.rhoc_),
-    Uc_(p.Uc_),
-    muc_(p.muc_)
+    UTurb_(p.UTurb_)
 {}
 
 
@@ -239,10 +238,7 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(p.rho_),
     age_(p.age_),
     tTurb_(p.tTurb_),
-    UTurb_(p.UTurb_),
-    rhoc_(p.rhoc_),
-    Uc_(p.Uc_),
-    muc_(p.muc_)
+    UTurb_(p.UTurb_)
 {}
 
 
@@ -259,14 +255,12 @@ bool Foam::KinematicParcel<ParcelType>::move
 {
     typename TrackCloudType::parcelType& p =
         static_cast<typename TrackCloudType::parcelType&>(*this);
-    typename TrackCloudType::particleType::trackingData& ttd =
-        static_cast<typename TrackCloudType::particleType::trackingData&>(td);
+    typename TrackCloudType::parcelType::trackingData& ttd =
+        static_cast<typename TrackCloudType::parcelType::trackingData&>(td);
 
     ttd.switchProcessor = false;
     ttd.keepParticle = true;
 
-    const polyMesh& mesh = cloud.pMesh();
-    const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
     const cloudSolution& solution = cloud.solution();
     const scalarField& cellLengthScale = cloud.cellLengthScale();
     const scalar maxCo = solution.maxCo();
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index 44d0514dfb..fe88741009 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -184,6 +184,18 @@ public:
                 autoPtr<interpolation<scalar>> muInterp_;
 
 
+            // Cached continuous phase properties
+
+                //- Density [kg/m3]
+                scalar rhoc_;
+
+                //- Velocity [m/s]
+                vector Uc_;
+
+                //- Viscosity [Pa.s]
+                scalar muc_;
+
+
             //- Local gravitational or other body-force acceleration
             const vector& g_;
 
@@ -219,6 +231,24 @@ public:
             //  phase dynamic viscosity field
             inline const interpolation<scalar>& muInterp() const;
 
+            //- Return the continuous phase density
+            inline scalar rhoc() const;
+
+            //- Access the continuous phase density
+            inline scalar& rhoc();
+
+            //- Return the continuous phase velocity
+            inline const vector& Uc() const;
+
+            //- Access the continuous phase velocity
+            inline vector& Uc();
+
+            //- Return the continuous phase viscosity
+            inline scalar muc() const;
+
+            //- Access the continuous phase viscosity
+            inline scalar& muc();
+
             // Return const access to the gravitational acceleration vector
             inline const vector& g() const;
 
@@ -268,18 +298,6 @@ protected:
             vector UTurb_;
 
 
-        // Cell-based quantities
-
-            //- Density [kg/m3]
-            scalar rhoc_;
-
-            //- Velocity [m/s]
-            vector Uc_;
-
-            //- Viscosity [Pa.s]
-            scalar muc_;
-
-
     // Protected Member Functions
 
         //- Calculate new particle velocity
@@ -459,15 +477,6 @@ public:
             //- Return const access to turbulent velocity fluctuation
             inline const vector& UTurb() const;
 
-            //- Return const access to carrier density [kg/m3]
-            inline scalar rhoc() const;
-
-            //- Return const access to carrier velocity [m/s]
-            inline const vector& Uc() const;
-
-            //- Return const access to carrier viscosity [Pa.s]
-            inline scalar muc() const;
-
 
         // Edit
 
@@ -505,7 +514,7 @@ public:
         // Helper functions
 
             //- Cell owner mass
-            inline scalar massCell(const label celli) const;
+            inline scalar massCell(const trackingData& td) const;
 
             //- Particle mass
             inline scalar mass() const;
@@ -532,31 +541,53 @@ public:
             inline static scalar areaS(const scalar d);
 
             //- Reynolds number
-            inline scalar Re
+            inline scalar Re(const trackingData& td) const;
+
+            //- Reynolds number for given conditions
+            inline static scalar Re
             (
-                const vector& U,        // particle velocity
-                const scalar d,         // particle diameter
-                const scalar rhoc,      // carrier density
-                const scalar muc        // carrier dynamic viscosity
-            ) const;
+                const scalar rhoc,
+                const vector& U,
+                const vector& Uc,
+                const scalar d,
+                const scalar muc
+            );
 
             //- Weber number
             inline scalar We
             (
-                const vector& U,        // particle velocity
-                const scalar d,         // particle diameter
-                const scalar rhoc,      // carrier density
-                const scalar sigma      // particle surface tension
+                const trackingData& td,
+                const scalar sigma
             ) const;
 
+            //- Weber number for given conditions
+            inline static scalar We
+            (
+                const scalar rhoc,
+                const vector& U,
+                const vector& Uc,
+                const scalar d,
+                const scalar sigma
+            );
+
             //- Eotvos number
             inline scalar Eo
             (
-                const vector& a,        // acceleration
-                const scalar d,         // particle diameter
-                const scalar sigma      // particle surface tension
+                const trackingData& td,
+                const scalar sigma
             ) const;
 
+            //- Eotvos number for given conditions
+            inline static scalar Eo
+            (
+                const vector& g,
+                const scalar rho,
+                const scalar rhoc,
+                const vector& U,
+                const scalar d,
+                const scalar sigma
+            );
+
 
         // Main calculation loop
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index 749c136883..d2fdfbf67a 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -89,10 +89,7 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(0.0),
     age_(0.0),
     tTurb_(0.0),
-    UTurb_(Zero),
-    rhoc_(0.0),
-    Uc_(Zero),
-    muc_(0.0)
+    UTurb_(Zero)
 {}
 
 
@@ -114,10 +111,7 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(0.0),
     age_(0.0),
     tTurb_(0.0),
-    UTurb_(Zero),
-    rhoc_(0.0),
-    Uc_(Zero),
-    muc_(0.0)
+    UTurb_(Zero)
 {}
 
 
@@ -147,10 +141,7 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(constProps.rho0()),
     age_(0.0),
     tTurb_(0.0),
-    UTurb_(Zero),
-    rhoc_(0.0),
-    Uc_(Zero),
-    muc_(0.0)
+    UTurb_(Zero)
 {}
 
 
@@ -269,28 +260,7 @@ inline const Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb() const
 
 
 template<class ParcelType>
-inline Foam::scalar Foam::KinematicParcel<ParcelType>::rhoc() const
-{
-    return rhoc_;
-}
-
-
-template<class ParcelType>
-inline const Foam::vector& Foam::KinematicParcel<ParcelType>::Uc() const
-{
-    return Uc_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar Foam::KinematicParcel<ParcelType>::muc() const
-{
-    return muc_;
-}
-
-
-template<class ParcelType>
-inline void Foam::KinematicParcel<ParcelType>::active(const bool state)
+inline bool& Foam::KinematicParcel<ParcelType>::active()
 {
     active_ = state;
 }
@@ -362,10 +332,10 @@ inline Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb()
 template<class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::massCell
 (
-    const label celli
+    const trackingData& td
 ) const
 {
-    return rhoc_*this->mesh().cellVolumes()[celli];
+    return td.rhoc()*this->mesh().cellVolumes()[this->cell()];
 }
 
 
@@ -428,39 +398,76 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::areaS(const scalar d)
 template<class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::Re
 (
-    const vector& U,
-    const scalar d,
-    const scalar rhoc,
-    const scalar muc
+    const trackingData& td
 ) const
 {
-    return rhoc*mag(U - Uc_)*d/(muc + ROOTVSMALL);
+    return Re(td.rhoc(), U_, td.Uc(), d_, td.muc());
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::KinematicParcel<ParcelType>::Re
+(
+    const scalar rhoc,
+    const vector& U,
+    const vector& Uc,
+    const scalar d,
+    const scalar muc
+)
+{
+    return rhoc*mag(U - Uc)*d/max(muc, ROOTVSMALL);
 }
 
 
 template<class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::We
 (
-    const vector& U,
-    const scalar d,
-    const scalar rhoc,
+    const trackingData& td,
     const scalar sigma
 ) const
 {
-    return rhoc*magSqr(U - Uc_)*d/(sigma + ROOTVSMALL);
+    return We(td.rhoc(), U_, td.Uc(), d_, sigma);
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::KinematicParcel<ParcelType>::We
+(
+    const scalar rhoc,
+    const vector& U,
+    const vector& Uc,
+    const scalar d,
+    const scalar sigma
+)
+{
+    return rhoc*magSqr(U - Uc)*d/max(sigma, ROOTVSMALL);
 }
 
 
 template<class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::Eo
 (
-    const vector& a,
-    const scalar d,
+    const trackingData& td,
     const scalar sigma
 ) const
 {
-    vector dir = U_/(mag(U_) + ROOTVSMALL);
-    return mag(a & dir)*(rho_ - rhoc_)*sqr(d)/(sigma + ROOTVSMALL);
+    return Eo(td.g(), rho_, td.rhoc(), U_, d_, sigma);
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::KinematicParcel<ParcelType>::Eo
+(
+    const vector& g,
+    const scalar rho,
+    const scalar rhoc,
+    const vector& U,
+    const scalar d,
+    const scalar sigma
+)
+{
+    const vector dir = U/max(mag(U), ROOTVSMALL);
+    return mag(g & dir)*(rho - rhoc)*sqr(d)/max(sigma, ROOTVSMALL);
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
index 7f4800fdb4..37fdfa8c00 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
@@ -41,7 +41,7 @@ Foam::string Foam::KinematicParcel<ParcelType>::propertyTypes_ =
 template<class ParcelType>
 const std::size_t Foam::KinematicParcel<ParcelType>::sizeofFields
 (
-    offsetof(KinematicParcel<ParcelType>, rhoc_)
+    sizeof(KinematicParcel<ParcelType>)
   - offsetof(KinematicParcel<ParcelType>, active_)
 );
 
@@ -66,10 +66,7 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(0.0),
     age_(0.0),
     tTurb_(0.0),
-    UTurb_(Zero),
-    rhoc_(0.0),
-    Uc_(Zero),
-    muc_(0.0)
+    UTurb_(Zero)
 {
     if (readFields)
     {
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelTrackingDataI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelTrackingDataI.H
index bf558611ae..72806a9fe3 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelTrackingDataI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelTrackingDataI.H
@@ -56,6 +56,9 @@ inline Foam::KinematicParcel<ParcelType>::trackingData::trackingData
             cloud.mu()
         )
     ),
+    rhoc_(Zero),
+    Uc_(Zero),
+    muc_(Zero),
     g_(cloud.g().value()),
     part_(part)
 {}
@@ -93,6 +96,50 @@ Foam::KinematicParcel<ParcelType>::trackingData::g() const
 }
 
 
+template<class ParcelType>
+inline Foam::scalar
+Foam::KinematicParcel<ParcelType>::trackingData::rhoc() const
+{
+    return rhoc_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::KinematicParcel<ParcelType>::trackingData::rhoc()
+{
+    return rhoc_;
+}
+
+
+template<class ParcelType>
+inline const Foam::vector&
+Foam::KinematicParcel<ParcelType>::trackingData::Uc() const
+{
+    return Uc_;
+}
+
+
+template<class ParcelType>
+inline Foam::vector& Foam::KinematicParcel<ParcelType>::trackingData::Uc()
+{
+    return Uc_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::KinematicParcel<ParcelType>::trackingData::muc() const
+{
+    return muc_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::KinematicParcel<ParcelType>::trackingData::muc()
+{
+    return muc_;
+}
+
+
 template<class ParcelType>
 inline typename Foam::KinematicParcel<ParcelType>::trackingData::trackPart
 Foam::KinematicParcel<ParcelType>::trackingData::part() const
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index aae4db0f73..09f3e10277 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -198,7 +198,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // Calc surface values
     scalar Ts, rhos, mus, Prs, kappas;
     this->calcSurfaceValues(cloud, td, T0, Ts, rhos, mus, Prs, kappas);
-    scalar Res = this->Re(U0, d0, rhos, mus);
+    scalar Res = this->Re(rhos, U0, td.Uc(), d0, mus);
 
 
     // Sources
@@ -390,7 +390,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 
     // Correct surface values due to emitted species
     this->correctSurfaceValues(cloud, td, Ts, Cs, rhos, mus, Prs, kappas);
-    Res = this->Re(U0, this->d_, rhos, mus);
+    Res = this->Re(rhos, U0, td.Uc(), this->d_, mus);
 
 
     // 3. Compute heat- and momentum transfers
@@ -572,7 +572,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     {
         // Molar average molecular weight of carrier mix
         const scalar Wc =
-            max(SMALL, this->rhoc_*RR*this->Tc_/this->pc_);
+            max(SMALL, td.rhoc()*RR*this->Tc_/this->pc_);
 
         // Note: hardcoded gaseous diffusivities for now
         // TODO: add to carrier thermo
@@ -649,7 +649,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
         T,
         this->Tc_,
         this->pc_,
-        this->rhoc_,
+        td.rhoc(),
         mass,
         YGas,
         YLiquid,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 8211e3eb6d..f66bfa9c89 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -116,7 +116,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     if (cloud.heatTransfer().BirdCorrection())
     {
         // Average molecular weight of carrier mix - assumes perfect gas
-        const scalar Wc = this->rhoc_*RR*this->Tc_/this->pc_;
+        const scalar Wc = td.rhoc()*RR*this->Tc_/this->pc_;
 
         forAll(dMassPC, i)
         {
@@ -250,12 +250,12 @@ void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
         return;
     }
 
-    const scalar massCell = this->massCell(this->cell());
+    const scalar massCell = this->massCell(td);
 
-    this->rhoc_ += addedMass/cloud.pMesh().cellVolumes()[this->cell()];
+    td.rhoc() += addedMass/cloud.pMesh().cellVolumes()[this->cell()];
 
     const scalar massCellNew = massCell + addedMass;
-    this->Uc_ = (this->Uc_*massCell + cloud.UTrans()[this->cell()])/massCellNew;
+    td.Uc() = (td.Uc()*massCell + cloud.UTrans()[this->cell()])/massCellNew;
 
     scalar CpEff = 0.0;
     forAll(cloud.rhoTrans(), i)
@@ -405,7 +405,7 @@ void Foam::ReactingParcel<ParcelType>::calc
     // Calc surface values
     scalar Ts, rhos, mus, Prs, kappas;
     this->calcSurfaceValues(cloud, td, T0, Ts, rhos, mus, Prs, kappas);
-    scalar Res = this->Re(U0, d0, rhos, mus);
+    scalar Res = this->Re(rhos, U0, td.Uc(), d0, mus);
 
 
     // Sources
@@ -519,7 +519,7 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     // Correct surface values due to emitted species
     correctSurfaceValues(cloud, td, Ts, Cs, rhos, mus, Prs, kappas);
-    Res = this->Re(U0, this->d_, rhos, mus);
+    Res = this->Re(rhos, U0, td.Uc(), this->d(), mus);
 
 
     // 3. Compute heat- and momentum transfers
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index d20900cbc3..006c45a9fd 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -73,7 +73,7 @@ void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
     const label celli = this->cell();
     const scalar massCell = this->massCell(celli);
 
-    this->Uc_ += cloud.UTrans()[celli]/massCell;
+    td.Uc() += cloud.UTrans()[celli]/massCell;
 
 //    tetIndices tetIs = this->currentTetIndices();
 //    Tc_ = td.TInterp().interpolate(this->coordinates(), tetIs);
@@ -127,7 +127,7 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
     // Assuming thermo props vary linearly with T for small d(T)
     const scalar TRatio = Tc_/Ts;
 
-    rhos = this->rhoc_*TRatio;
+    rhos = td.rhoc()*TRatio;
 
     tetIndices tetIs = this->currentTetIndices();
     mus = td.muInterp().interpolate(this->coordinates(), tetIs)/TRatio;
@@ -162,7 +162,7 @@ void Foam::ThermoParcel<ParcelType>::calc
     calcSurfaceValues(cloud, td, this->T_, Ts, rhos, mus, Pr, kappas);
 
     // Reynolds number
-    scalar Re = this->Re(this->U_, this->d_, rhos, mus);
+    scalar Re = this->Re(rhos, this->U_, td.Uc(), this->d_, mus);
 
 
     // Sources
diff --git a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C
index 9ecc679c24..e74695516d 100644
--- a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C
+++ b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,6 +148,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::ParticleForceList<CloudType>::calcCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
@@ -160,7 +161,7 @@ Foam::forceSuSp Foam::ParticleForceList<CloudType>::calcCoupled
     {
         forAll(*this, i)
         {
-            value += this->operator[](i).calcCoupled(p, dt, mass, Re, muc);
+            value += this->operator[](i).calcCoupled(p, td, dt, mass, Re, muc);
         }
     }
 
@@ -172,6 +173,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::ParticleForceList<CloudType>::calcNonCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
@@ -184,7 +186,8 @@ Foam::forceSuSp Foam::ParticleForceList<CloudType>::calcNonCoupled
     {
         forAll(*this, i)
         {
-            value += this->operator[](i).calcNonCoupled(p, dt, mass, Re, muc);
+            value +=
+                this->operator[](i).calcNonCoupled(p, td, dt, mass, Re, muc);
         }
     }
 
@@ -196,13 +199,14 @@ template<class CloudType>
 Foam::scalar Foam::ParticleForceList<CloudType>::massEff
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar mass
 ) const
 {
     scalar massEff = mass;
     forAll(*this, i)
     {
-        massEff += this->operator[](i).massAdd(p, mass);
+        massEff += this->operator[](i).massAdd(p, td, mass);
     }
 
     return massEff;
diff --git a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H
index 9b115cbf88..e068b8199c 100644
--- a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H
+++ b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,6 +127,7 @@ public:
             virtual forceSuSp calcCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
@@ -137,6 +138,7 @@ public:
             virtual forceSuSp calcNonCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
@@ -147,6 +149,7 @@ public:
             virtual scalar massEff
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar mass
             ) const;
 };
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C
index 7d66ae2b37..fc0e07454b 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,6 +81,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::DistortedSphereDragForce<CloudType>::calcCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H
index 2852b631b3..4a0d0bb163 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,6 +105,7 @@ public:
             virtual forceSuSp calcCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.C
index fde0072ae2..de17e60966 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,6 +96,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::ErgunWenYuDragForce<CloudType>::calcCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.H
index 4dcf5fb9e4..5e4476cff2 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,6 +105,7 @@ public:
             virtual forceSuSp calcCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.C
index 674d4b1824..9bfa59b048 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -89,6 +89,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::NonSphereDragForce<CloudType>::calcCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.H
index b10a28748f..fcb1a4ddb4 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -145,6 +145,7 @@ public:
             virtual forceSuSp calcCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.C
index 2afe4f620f..632d9b4262 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,6 +96,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::PlessisMasliyahDragForce<CloudType>::calcCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H
index 8ea4f2d663..c46dd94662 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/PlessisMasliyahDrag/PlessisMasliyahDragForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,6 +105,7 @@ public:
             virtual forceSuSp calcCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.C
index b93c2650aa..00f0b621b1 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,6 +78,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.H
index 552db6e686..d16e19c005 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -97,6 +97,7 @@ public:
             virtual forceSuSp calcCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.C
index c6747ce768..7fe9d28e2b 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -93,6 +93,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::WenYuDragForce<CloudType>::calcCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H
index f61cc6d86f..7a03000974 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/WenYuDrag/WenYuDragForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,6 +105,7 @@ public:
             virtual forceSuSp calcCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Gravity/GravityForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Gravity/GravityForce.C
index 54f7dc9c63..d73a78b249 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Gravity/GravityForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Gravity/GravityForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,6 +61,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::GravityForce<CloudType>::calcNonCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
@@ -69,7 +70,7 @@ Foam::forceSuSp Foam::GravityForce<CloudType>::calcNonCoupled
 {
     forceSuSp value(Zero, 0.0);
 
-    value.Su() = mass*g_*(1.0 - p.rhoc()/p.rho());
+    value.Su() = mass*g_*(1.0 - td.rhoc()/p.rho());
 
     return value;
 }
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Gravity/GravityForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Gravity/GravityForce.H
index 9ba6251751..43e55523bd 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Gravity/GravityForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Gravity/GravityForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -110,6 +110,7 @@ public:
             virtual forceSuSp calcNonCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C
index 1268179b73..cc45ca417d 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C
@@ -32,6 +32,7 @@ template<class CloudType>
 Foam::scalar Foam::LiftForce<CloudType>::LiftForce::Cl
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const vector& curlUc,
     const scalar Re,
     const scalar muc
@@ -128,6 +129,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::LiftForce<CloudType>::calcCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
@@ -139,9 +141,9 @@ Foam::forceSuSp Foam::LiftForce<CloudType>::calcCoupled
     vector curlUc =
         curlUcInterp().interpolate(p.coordinates(), p.currentTetIndices());
 
-    scalar Cl = this->Cl(p, curlUc, Re, muc);
+    scalar Cl = this->Cl(p, td, curlUc, Re, muc);
 
-    value.Su() = mass/p.rho()*p.rhoc()*Cl*((p.Uc() - p.U())^curlUc);
+    value.Su() = mass/p.rho()*td.rhoc()*Cl*((td.Uc() - p.U())^curlUc);
 
     return value;
 }
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.H
index d03eaf4b77..79063a04c5 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,6 +74,7 @@ protected:
         virtual scalar Cl
         (
             const typename CloudType::parcelType& p,
+            const typename CloudType::parcelType::trackingData& td,
             const vector& curlUc,
             const scalar Re,
             const scalar muc
@@ -127,6 +128,7 @@ public:
             virtual forceSuSp calcCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/SaffmanMeiLift/SaffmanMeiLiftForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/SaffmanMeiLift/SaffmanMeiLiftForce.C
index 26da78568f..04a3615899 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/SaffmanMeiLift/SaffmanMeiLiftForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/SaffmanMeiLift/SaffmanMeiLiftForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,12 +34,13 @@ template<class CloudType>
 Foam::scalar Foam::SaffmanMeiLiftForce<CloudType>::SaffmanMeiLiftForce::Cl
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const vector& curlUc,
     const scalar Re,
     const scalar muc
 ) const
 {
-    scalar Rew = p.rhoc()*mag(curlUc)*sqr(p.d())/(muc + ROOTVSMALL);
+    scalar Rew = td.rhoc()*mag(curlUc)*sqr(p.d())/(muc + ROOTVSMALL);
     scalar beta = 0.5*(Rew/(Re + ROOTVSMALL));
     scalar alpha = 0.3314*sqrt(beta);
     scalar f = (1.0 - alpha)*exp(-0.1*Re) + alpha;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/SaffmanMeiLift/SaffmanMeiLiftForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/SaffmanMeiLift/SaffmanMeiLiftForce.H
index b5d12560e1..1563b19d8e 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/SaffmanMeiLift/SaffmanMeiLiftForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/SaffmanMeiLift/SaffmanMeiLiftForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -62,6 +62,7 @@ protected:
         virtual scalar Cl
         (
             const typename CloudType::parcelType& p,
+            const typename CloudType::parcelType::trackingData& td,
             const vector& curlUc,
             const scalar Re,
             const scalar muc
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/TomiyamaLift/TomiyamaLiftForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/TomiyamaLift/TomiyamaLiftForce.C
index 87c482bec5..5bb1d84065 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/TomiyamaLift/TomiyamaLiftForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/TomiyamaLift/TomiyamaLiftForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,7 @@ template<class CloudType>
 Foam::scalar Foam::TomiyamaLiftForce<CloudType>::TomiyamaLiftForce::Cl
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const vector& curlUc,
     const scalar Re,
     const scalar muc
@@ -38,9 +39,9 @@ Foam::scalar Foam::TomiyamaLiftForce<CloudType>::TomiyamaLiftForce::Cl
 {
     const vector& g = this->owner().g().value();
 
-    scalar Eo = p.Eo(g, p.d(), sigma_);
+    scalar Eo = p.Eo(td, sigma_);
     scalar dH = p.d()*cbrt(1.0 + 0.163*pow(Eo, 0.757));
-    scalar Eod = p.Eo(g, dH, sigma_);
+    scalar Eod = p.Eo(g, p.rho(), td.rhoc(), p.U(), dH, sigma_);
     scalar f = 0.00105*pow3(Eod) - 0.0159*sqr(Eod) - 0.0204*Eod + 0.474;
 
     if (Eod <= 4)
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/TomiyamaLift/TomiyamaLiftForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/TomiyamaLift/TomiyamaLiftForce.H
index 4ce9f7874a..01c8cffd90 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/TomiyamaLift/TomiyamaLiftForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/TomiyamaLift/TomiyamaLiftForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,7 @@ protected:
         virtual scalar Cl
         (
             const typename CloudType::parcelType& p,
+            const typename CloudType::parcelType::trackingData& td,
             const vector& curlUc,
             const scalar Re,
             const scalar muc
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/NonInertialFrame/NonInertialFrameForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/NonInertialFrame/NonInertialFrameForce.C
index 23b5045a7e..dec9a71663 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/NonInertialFrame/NonInertialFrameForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/NonInertialFrame/NonInertialFrameForce.C
@@ -181,6 +181,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::NonInertialFrameForce<CloudType>::calcNonCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/NonInertialFrame/NonInertialFrameForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/NonInertialFrame/NonInertialFrameForce.H
index 01b3592532..a1afeb4ceb 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/NonInertialFrame/NonInertialFrameForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/NonInertialFrame/NonInertialFrameForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -143,6 +143,7 @@ public:
             virtual forceSuSp calcNonCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.C
index a268a07b09..b0fc66c7dc 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.C
@@ -97,6 +97,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::ParamagneticForce<CloudType>::calcNonCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.H
index be7f0d409e..410117729c 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Paramagnetic/ParamagneticForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -123,6 +123,7 @@ public:
             virtual forceSuSp calcNonCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.C
index b337b4222f..e0b8d6fb3e 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,6 +79,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::ParticleForce<CloudType>::calcCoupled
 (
     const typename CloudType::parcelType&,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
@@ -97,6 +98,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::ParticleForce<CloudType>::calcNonCoupled
 (
     const typename CloudType::parcelType&,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
@@ -115,6 +117,7 @@ template<class CloudType>
 Foam::scalar Foam::ParticleForce<CloudType>::massAdd
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar mass
 ) const
 {
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.H
index 872d7f87f6..23ab368c53 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/ParticleForce/ParticleForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -158,6 +158,7 @@ public:
             virtual forceSuSp calcCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
@@ -168,6 +169,7 @@ public:
             virtual forceSuSp calcNonCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
@@ -178,6 +180,7 @@ public:
             virtual scalar massAdd
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar mass
             ) const;
 };
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C
index 0059a6627f..2d1b322950 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.C
@@ -119,6 +119,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::PressureGradientForce<CloudType>::calcCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
@@ -130,7 +131,7 @@ Foam::forceSuSp Foam::PressureGradientForce<CloudType>::calcCoupled
     vector DUcDt =
         DUcDtInterp().interpolate(p.coordinates(), p.currentTetIndices());
 
-    value.Su() = mass*p.rhoc()/p.rho()*DUcDt;
+    value.Su() = mass*td.rhoc()/p.rho()*DUcDt;
 
     return value;
 }
@@ -140,6 +141,7 @@ template<class CloudType>
 Foam::scalar Foam::PressureGradientForce<CloudType>::massAdd
 (
     const typename CloudType::parcelType&,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar
 ) const
 {
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H
index fb91c626b0..737b2c1797 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/PressureGradient/PressureGradientForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -119,6 +119,7 @@ public:
             virtual forceSuSp calcCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
@@ -129,6 +130,7 @@ public:
             virtual scalar massAdd
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar mass
             ) const;
 };
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/SRF/SRFForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/SRF/SRFForce.C
index 08f065539e..83443fda39 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/SRF/SRFForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/SRF/SRFForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,6 +80,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::SRFForce<CloudType>::calcNonCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
@@ -96,7 +97,7 @@ Foam::forceSuSp Foam::SRFForce<CloudType>::calcNonCoupled
 
     // Coriolis and centrifugal acceleration terms
     value.Su() =
-        mass*(1.0 - p.rhoc()/p.rho())
+        mass*(1.0 - td.rhoc()/p.rho())
        *(2.0*(p.U() ^ omega) + (omega ^ (r ^ omega)));
 
     return value;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/SRF/SRFForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/SRF/SRFForce.H
index 1cb72cf4d4..d0472f3dc7 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/SRF/SRFForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/SRF/SRFForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,6 +107,7 @@ public:
             virtual forceSuSp calcNonCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C
index c701e07a00..54408f9b24 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,6 +72,7 @@ template<class CloudType>
 Foam::forceSuSp Foam::VirtualMassForce<CloudType>::calcCoupled
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar dt,
     const scalar mass,
     const scalar Re,
@@ -79,7 +80,7 @@ Foam::forceSuSp Foam::VirtualMassForce<CloudType>::calcCoupled
 ) const
 {
     forceSuSp value =
-        PressureGradientForce<CloudType>::calcCoupled(p, dt, mass, Re, muc);
+        PressureGradientForce<CloudType>::calcCoupled(p, td, dt, mass, Re, muc);
 
     value.Su() *= Cvm_;
 
@@ -91,10 +92,11 @@ template<class CloudType>
 Foam::scalar Foam::VirtualMassForce<CloudType>::massAdd
 (
     const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
     const scalar mass
 ) const
 {
-    return mass*p.rhoc()/p.rho()*Cvm_;
+    return mass*td.rhoc()/p.rho()*Cvm_;
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H
index b1bb602bca..f9a4d540dc 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/VirtualMass/VirtualMassForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,6 +105,7 @@ public:
             virtual forceSuSp calcCoupled
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar dt,
                 const scalar mass,
                 const scalar Re,
@@ -115,6 +116,7 @@ public:
             virtual scalar massAdd
             (
                 const typename CloudType::parcelType& p,
+                const typename CloudType::parcelType::trackingData& td,
                 const scalar mass
             ) const;
 };
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C
index 9667fcb03b..cd27aab7ba 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C
@@ -164,11 +164,11 @@ Foam::vector Foam::PackingModels::Explicit<CloudType>::velocityCorrection
     vector dU = Zero;
 
     //// existing forces
-    //const scalar Re = p.Re(p.U(), p.d(), p.rhoc(), p.muc());
+    //const scalar Re = p.Re(td);
     //const typename CloudType::forceType& forces = this->owner().forces();
     //const forceSuSp F =
-    //    forces.calcCoupled(p, deltaT, p.mass(), Re, p.muc())
-    //  + forces.calcNonCoupled(p, deltaT, p.mass(), Re, p.muc());
+    //    forces.calcCoupled(p, td, deltaT, p.mass(), Re, td.muc())
+    //  + forces.calcNonCoupled(p, td, deltaT, p.mass(), Re, td.muc());
 
     // correction velocity
     if ((uRelative & alphaGrad) > 0)
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
index 3b5d63db13..229f47b698 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
@@ -164,7 +164,7 @@ void Foam::SprayParcel<ParcelType>::calcAtomization
     const auto& liquids = composition.liquids();
 
     // Average molecular weight of carrier mix - assumes perfect gas
-    scalar Wc = this->rhoc_*RR*this->Tc()/this->pc();
+    scalar Wc = td.rhoc()*RR*this->Tc()/this->pc();
     scalar R = RR/Wc;
     scalar Tav = atomization.Taverage(this->T(), this->Tc());
 
@@ -229,32 +229,27 @@ void Foam::SprayParcel<ParcelType>::calcBreakup
         return;
     }
 
-    typedef typename TrackCloudType::parcelType parcelType;
-    typedef typename TrackCloudType::forceType forceType;
-
-    const parcelType& p = static_cast<const parcelType&>(*this);
-    const forceType& forces = cloud.forces();
-
     if (breakup.solveOscillationEq())
     {
         solveTABEq(cloud, td, dt);
     }
 
     // Average molecular weight of carrier mix - assumes perfect gas
-    scalar Wc = this->rhoc()*RR*this->Tc()/this->pc();
+    scalar Wc = td.rhoc()*RR*this->Tc()/this->pc();
     scalar R = RR/Wc;
     scalar Tav = cloud.atomization().Taverage(this->T(), this->Tc());
 
     // Calculate average gas density based on average temperature
     scalar rhoAv = this->pc()/(R*Tav);
-    scalar muAv = this->muc();
-    vector Urel = this->U() - this->Uc();
+    scalar muAv = td.muc();
+    vector Urel = this->U() - td.Uc();
     scalar Urmag = mag(Urel);
-    scalar Re = this->Re(this->U(), this->d(), rhoAv, muAv);
+    scalar Re = this->Re(rhoAv, this->U(), td.Uc(), this->d(), muAv);
 
     const scalar mass = p.mass();
-    const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, muAv);
-    const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, muAv);
+    const typename TrackCloudType::forceType& forces = cloud.forces();
+    const forceSuSp Fcp = forces.calcCoupled(p, ttd, dt, mass, Re, muAv);
+    const forceSuSp Fncp = forces.calcNonCoupled(p, ttd, dt, mass, Re, muAv);
     this->tMom() = mass/(Fcp.Sp() + Fncp.Sp() + ROOTVSMALL);
 
     const vector g = cloud.g().value();
@@ -301,9 +296,9 @@ void Foam::SprayParcel<ParcelType>::calcBreakup
         child->nParticle() = parcelMassChild/massChild;
 
         const forceSuSp Fcp =
-            forces.calcCoupled(*child, dt, massChild, Re, muAv);
+            forces.calcCoupled(*child, ttd, dt, massChild, Re, muAv);
         const forceSuSp Fncp =
-            forces.calcNonCoupled(*child, dt, massChild, Re, muAv);
+            forces.calcNonCoupled(*child, ttd, dt, massChild, Re, muAv);
 
         child->age() = 0.0;
         child->liquidCore() = 0.0;
@@ -392,8 +387,8 @@ void Foam::SprayParcel<ParcelType>::solveTABEq
     if (omega2 > 0)
     {
         scalar omega = sqrt(omega2);
-        scalar rhoc = this->rhoc();
-        scalar We = this->We(this->U(), r, rhoc, sigma_)/TABtwoWeCrit;
+        scalar We =
+            this->We(td.rhoc(), this->U(), td.Uc(), r, sigma_)/TABtwoWeCrit;
 
         // Initial values for y and yDot
         scalar y0 = this->y() - We;