diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 43622464153c435595c8cdaa8e133335204af336..3e229fa1e65a3cd8f321ed4a9f970fce25819ae5 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -212,7 +212,7 @@ void Foam::KinematicCloud<CloudType>::evolveCloud
         // before it is required.
         cloud.motion(cloud, td);
 
-        stochasticCollision().update(solution_.trackTime());
+        stochasticCollision().update(td, solution_.trackTime());
     }
     else
     {
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index 577c150706d87d3b2dc63c138b54d3bfa3075cd8..fc52c1abaf9957cbd00b0b5cd01786bad9b9372a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -221,7 +221,6 @@ void Foam::ReactingCloud<CloudType>::setParcelThermoProperties
 {
     CloudType::setParcelThermoProperties(parcel, lagrangianDt);
 
-    parcel.pc() = this->thermo().thermo().p()[parcel.cell()];
     parcel.Y() = composition().YMixture0();
 }
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 77dab8ad9f5914378a140326b2f0dfd7892ae849..bc5d13c69efa8ced5929ff027b2cb07a85be6a13 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -42,8 +42,7 @@ template<class TrackCloudType>
 void Foam::KinematicParcel<ParcelType>::setCellValues
 (
     TrackCloudType& cloud,
-    trackingData& td,
-    const scalar dt
+    trackingData& td
 )
 {
     tetIndices tetIs = this->currentTetIndices();
@@ -65,8 +64,18 @@ void Foam::KinematicParcel<ParcelType>::setCellValues
     td.Uc() = td.UInterp().interpolate(this->coordinates(), tetIs);
 
     td.muc() = td.muInterp().interpolate(this->coordinates(), tetIs);
+}
 
-    // Apply dispersion components to carrier phase velocity
+
+template<class ParcelType>
+template<class TrackCloudType>
+void Foam::KinematicParcel<ParcelType>::calcDispersion
+(
+    TrackCloudType& cloud,
+    trackingData& td,
+    const scalar dt
+)
+{
     td.Uc() = cloud.dispersion().update
     (
         dt,
@@ -307,7 +316,9 @@ bool Foam::KinematicParcel<ParcelType>::move
         if (dt > ROOTVSMALL)
         {
             // Update cell based properties
-            p.setCellValues(cloud, ttd, dt);
+            p.setCellValues(cloud, ttd);
+
+            p.calcDispersion(cloud, ttd, dt);
 
             if (solution.cellValueSourceCorrection())
             {
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index fe8874100919d1c691db0d57d61bd435f4aca02d..e32748ae92a07787c8f9c16cfa2de9dbc0fff67e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -593,7 +593,12 @@ public:
 
             //- Set cell values
             template<class TrackCloudType>
-            void setCellValues
+            void setCellValues(TrackCloudType& cloud, trackingData& td);
+
+            //- Apply dispersion to the carrier phase velocity and update
+            //  parcel turbulence parameters
+            template<class TrackCloudType>
+            void calcDispersion
             (
                 TrackCloudType& cloud,
                 trackingData& td,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 09f3e10277041b7d3fc59388ad0d01b605dea97f..ed8318074cc891f47c2623e1bd818f2fb14a5a90 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -143,11 +143,10 @@ template<class TrackCloudType>
 void Foam::ReactingMultiphaseParcel<ParcelType>::setCellValues
 (
     TrackCloudType& cloud,
-    trackingData& td,
-    const scalar dt
+    trackingData& td
 )
 {
-    ParcelType::setCellValues(cloud, td, dt);
+    ParcelType::setCellValues(cloud, td);
 }
 
 
@@ -187,7 +186,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     const scalar T0 = this->T_;
     const scalar mass0 = this->mass();
 
-    const scalar pc = this->pc_;
+    const scalar pc = td.pc();
 
     const scalarField& YMix = this->Y_;
     const label idG = composition.idGas();
@@ -571,8 +570,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     if (cloud.heatTransfer().BirdCorrection())
     {
         // Molar average molecular weight of carrier mix
-        const scalar Wc =
-            max(SMALL, td.rhoc()*RR*this->Tc_/this->pc_);
+        const scalar Wc = max(SMALL, td.rhoc()*RR*td.Tc()/td.pc());
 
         // Note: hardcoded gaseous diffusivities for now
         // TODO: add to carrier thermo
@@ -581,7 +579,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
         forAll(dMassDV, i)
         {
             const label id = composition.localToCarrierId(GAS, i);
-            const scalar Cp = composition.carrier().Cp(id, this->pc_, Ts);
+            const scalar Cp = composition.carrier().Cp(id, td.pc(), Ts);
             const scalar W = composition.carrier().W(id);
             const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 
@@ -589,7 +587,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
             const scalar Dab =
                 3.6059e-3*(pow(1.8*Ts, 1.75))
                *sqrt(1.0/W + 1.0/Wc)
-               /(this->pc_*beta);
+               /(td.pc()*beta);
 
             N += Ni;
             NCpW += Ni*Cp*W;
@@ -647,8 +645,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
         this->cell(),
         d,
         T,
-        this->Tc_,
-        this->pc_,
+        td.Tc(),
+        td.pc(),
         td.rhoc(),
         mass,
         YGas,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index d10af101005e73e7df5f6b5eadb9ac7fa856ffca..9fef053dd9469c3c10840f59897c3cc8d854f77f 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -421,12 +421,7 @@ public:
 
             //- Set cell values
             template<class TrackCloudType>
-            void setCellValues
-            (
-                TrackCloudType& cloud,
-                trackingData& td,
-                const scalar dt
-            );
+            void setCellValues(TrackCloudType& cloud, trackingData& td);
 
             //- Correct cell values using latest transfer information
             template<class TrackCloudType>
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index f66bfa9c89cb303aeb4d61a81dd45021b2d529cb..99a8679d73fb92d27030b04be707d1694db96ca7 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -74,7 +74,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
         return;
     }
 
-    const scalar TMax = phaseChange.TMax(pc_, X);
+    const scalar TMax = phaseChange.TMax(td.pc(), X);
     const scalar Tdash = min(T, TMax);
     const scalar Tsdash = min(Ts, TMax);
 
@@ -89,8 +89,8 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
         nus,
         Tdash,
         Tsdash,
-        pc_,
-        this->Tc_,
+        td.pc(),
+        td.Tc(),
         X,
         dMassPC
     );
@@ -107,7 +107,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     {
         const label cid = composition.localToCarrierId(idPhase, i);
 
-        const scalar dh = phaseChange.dh(cid, i, pc_, Tdash);
+        const scalar dh = phaseChange.dh(cid, i, td.pc(), Tdash);
         Sh -= dMassPC[i]*dh/dt;
     }
 
@@ -116,18 +116,18 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     if (cloud.heatTransfer().BirdCorrection())
     {
         // Average molecular weight of carrier mix - assumes perfect gas
-        const scalar Wc = td.rhoc()*RR*this->Tc_/this->pc_;
+        const scalar Wc = td.rhoc()*RR*td.Tc()/td.pc();
 
         forAll(dMassPC, i)
         {
             const label cid = composition.localToCarrierId(idPhase, i);
 
-            const scalar Cp = composition.carrier().Cp(cid, pc_, Tsdash);
+            const scalar Cp = composition.carrier().Cp(cid, td.pc(), Tsdash);
             const scalar W = composition.carrier().W(cid);
             const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
 
             const scalar Dab =
-                composition.liquids().properties()[i].D(pc_, Tsdash, Wc);
+                composition.liquids().properties()[i].D(td.pc(), Tsdash, Wc);
 
             // Molar flux of species coming from the particle (kmol/m^2/s)
             N += Ni;
@@ -175,8 +175,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
 :
     ParcelType(p),
     mass0_(p.mass0_),
-    Y_(p.Y_),
-    pc_(p.pc_)
+    Y_(p.Y_)
 {}
 
 
@@ -189,8 +188,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
 :
     ParcelType(p, mesh),
     mass0_(p.mass0_),
-    Y_(p.Y_),
-    pc_(p.pc_)
+    Y_(p.Y_)
 {}
 
 
@@ -201,19 +199,18 @@ template<class TrackCloudType>
 void Foam::ReactingParcel<ParcelType>::setCellValues
 (
     TrackCloudType& cloud,
-    trackingData& td,
-    const scalar dt
+    trackingData& td
 )
 {
-    ParcelType::setCellValues(cloud, td, dt);
+    ParcelType::setCellValues(cloud, td);
 
-    pc_ = td.pInterp().interpolate
+    td.pc() = td.pInterp().interpolate
     (
         this->coordinates(),
         this->currentTetIndices()
     );
 
-    if (pc_ < cloud.constProps().pMin())
+    if (td.pc() < cloud.constProps().pMin())
     {
         if (debug)
         {
@@ -222,7 +219,7 @@ void Foam::ReactingParcel<ParcelType>::setCellValues
                 << " to " << cloud.constProps().pMin() <<  nl << endl;
         }
 
-        pc_ = cloud.constProps().pMin();
+        td.pc() = cloud.constProps().pMin();
     }
 }
 
@@ -261,20 +258,15 @@ void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
     forAll(cloud.rhoTrans(), i)
     {
         scalar Y = cloud.rhoTrans(i)[this->cell()]/addedMass;
-        CpEff += Y*cloud.composition().carrier().Cp
-        (
-            i,
-            this->pc_,
-            this->Tc_
-        );
+        CpEff += Y*cloud.composition().carrier().Cp(i, td.pc(), td.Tc());
     }
 
     const scalar Cpc = td.CpInterp().psi()[this->cell()];
-    this->Cpc_ = (massCell*Cpc + addedMass*CpEff)/massCellNew;
+    td.Cpc() = (massCell*Cpc + addedMass*CpEff)/massCellNew;
 
-    this->Tc_ += cloud.hsTrans()[this->cell()]/(this->Cpc_*massCellNew);
+    td.Tc() += cloud.hsTrans()[this->cell()]/(td.Cpc()*massCellNew);
 
-    if (this->Tc_ < cloud.constProps().TMin())
+    if (td.Tc() < cloud.constProps().TMin())
     {
         if (debug)
         {
@@ -283,7 +275,7 @@ void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
                 << " to " << cloud.constProps().TMin() <<  nl << endl;
         }
 
-        this->Tc_ = cloud.constProps().TMin();
+        td.Tc() = cloud.constProps().TMin();
     }
 }
 
@@ -320,10 +312,10 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
     Xinf /= sum(Xinf);
 
     // Molar fraction of far field species at particle surface
-    const scalar Xsff = 1.0 - min(sum(Cs)*RR*this->T_/pc_, 1.0);
+    const scalar Xsff = 1.0 - min(sum(Cs)*RR*this->T_/td.pc(), 1.0);
 
     // Surface carrier total molar concentration
-    const scalar CsTot = pc_/(RR*this->T_);
+    const scalar CsTot = td.pc()/(RR*this->T_);
 
     // Surface carrier composition (molar fraction)
     scalarField Xs(Xinf.size());
@@ -357,9 +349,9 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
         const scalar cbrtW = cbrt(W);
 
         rhos += Xs[i]*W;
-        mus += Ys[i]*sqrtW*thermo.carrier().mu(i, pc_, T);
-        kappas += Ys[i]*cbrtW*thermo.carrier().kappa(i, pc_, T);
-        Cps += Xs[i]*thermo.carrier().Cp(i, pc_, T);
+        mus += Ys[i]*sqrtW*thermo.carrier().mu(i, td.pc(), T);
+        kappas += Ys[i]*cbrtW*thermo.carrier().kappa(i, td.pc(), T);
+        Cps += Xs[i]*thermo.carrier().Cp(i, td.pc(), T);
 
         sumYiSqrtW += Ys[i]*sqrtW;
         sumYiCbrtW += Ys[i]*cbrtW;
@@ -367,7 +359,7 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
 
     Cps = max(Cps, ROOTVSMALL);
 
-    rhos *= pc_/(RR*T);
+    rhos *= td.pc()/(RR*T);
     rhos = max(rhos, ROOTVSMALL);
 
     mus /= sumYiSqrtW;
@@ -478,7 +470,7 @@ void Foam::ReactingParcel<ParcelType>::calc
     scalarField dMass(dMassPC);
     scalar mass1 = updateMassFraction(mass0, dMass, Y_);
 
-    this->Cp_ = composition.Cp(0, Y_, pc_, T0);
+    this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);
 
     // Update particle density or diameter
     if (cloud.constProps().constantVolume())
@@ -504,7 +496,7 @@ void Foam::ReactingParcel<ParcelType>::calc
             {
                 scalar dmi = dm*Y_[i];
                 label gid = composition.localToCarrierId(0, i);
-                scalar hs = composition.carrier().Hs(gid, pc_, T0);
+                scalar hs = composition.carrier().Hs(gid, td.pc(), T0);
 
                 cloud.rhoTrans(gid)[this->cell()] += dmi;
                 cloud.hsTrans()[this->cell()] += dmi*hs;
@@ -544,7 +536,7 @@ void Foam::ReactingParcel<ParcelType>::calc
             Sph
         );
 
-    this->Cp_ = composition.Cp(0, Y_, pc_, T0);
+    this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);
 
 
     // Motion
@@ -565,7 +557,7 @@ void Foam::ReactingParcel<ParcelType>::calc
         {
             scalar dm = np0*dMass[i];
             label gid = composition.localToCarrierId(0, i);
-            scalar hs = composition.carrier().Hs(gid, pc_, T0);
+            scalar hs = composition.carrier().Hs(gid, td.pc(), T0);
 
             cloud.rhoTrans(gid)[this->cell()] += dm;
             cloud.UTrans()[this->cell()] += dm*U0;
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index 00ceaa00341d08e2ebd805540bdbdcf5ad7889ad..812b2c7ef55d49d58d7a5d5ef2ad338c5857098e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -128,6 +128,12 @@ public:
                 autoPtr<interpolation<scalar>> pInterp_;
 
 
+            // Cached continuous phase properties
+
+                //- Pressure [Pa]
+                scalar pc_;
+
+
     public:
 
         typedef typename ParcelType::trackingData::trackPart trackPart;
@@ -148,6 +154,12 @@ public:
             //- Return const access to the interpolator for continuous phase
             //  pressure field
             inline const interpolation<scalar>& pInterp() const;
+
+            //- Return the continuous phase pressure
+            inline scalar pc() const;
+
+            //- Access the continuous phase pressure
+            inline scalar& pc();
     };
 
 
@@ -164,12 +176,6 @@ protected:
             scalarField Y_;
 
 
-        // Cell-based quantities
-
-            //- Pressure [Pa]
-            scalar pc_;
-
-
     // Protected Member Functions
 
         //- Calculate Phase change
@@ -336,12 +342,6 @@ public:
             //- Return const access to mass fractions of mixture []
             inline const scalarField& Y() const;
 
-            //- Return the owner cell pressure [Pa]
-            inline scalar pc() const;
-
-            //- Return reference to the owner cell pressure [Pa]
-            inline scalar& pc();
-
 
         // Edit
 
@@ -356,12 +356,7 @@ public:
 
             //- Set cell values
             template<class TrackCloudType>
-            void setCellValues
-            (
-                TrackCloudType& cloud,
-                trackingData& td,
-                const scalar dt
-            );
+            void setCellValues(TrackCloudType& cloud, trackingData& td);
 
             //- Correct cell values using latest transfer information
             template<class TrackCloudType>
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
index 08699da9b876351f0dffd55f8f6b7de36e4fe9c3..9108f85e0934bf152dbafd1a7496534888f5e555 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
@@ -71,8 +71,7 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
 :
     ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
     mass0_(0.0),
-    Y_(0),
-    pc_(0.0)
+    Y_(0)
 {}
 
 
@@ -86,8 +85,7 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
 :
     ParcelType(mesh, position, celli),
     mass0_(0.0),
-    Y_(0),
-    pc_(0.0)
+    Y_(0)
 {}
 
 
@@ -129,8 +127,7 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
         constProps
     ),
     mass0_(0.0),
-    Y_(Y0),
-    pc_(0.0)
+    Y_(Y0)
 {
     // Set initial parcel mass
     mass0_ = this->mass();
@@ -171,20 +168,6 @@ inline const Foam::scalarField& Foam::ReactingParcel<ParcelType>::Y() const
 }
 
 
-template<class ParcelType>
-inline Foam::scalar Foam::ReactingParcel<ParcelType>::pc() const
-{
-    return pc_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar& Foam::ReactingParcel<ParcelType>::pc()
-{
-    return pc_;
-}
-
-
 template<class ParcelType>
 inline Foam::scalar& Foam::ReactingParcel<ParcelType>::mass0()
 {
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
index dccbd67800118545faece669c933025ca2f327a8..9541a85df8ad12d61a6b7b3d2af354889a6dd04b 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
@@ -55,8 +55,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
 :
     ParcelType(mesh, is, readFields),
     mass0_(0.0),
-    Y_(0),
-    pc_(0.0)
+    Y_(0)
 {
     if (readFields)
     {
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelTrackingDataI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelTrackingDataI.H
index cee1de2b91750a17ca6d40b97ef23ab3f6340140..80bccae5d3bc2e4d9aa6ea59b202c74e87d5f8d8 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelTrackingDataI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelTrackingDataI.H
@@ -39,7 +39,8 @@ inline Foam::ReactingParcel<ParcelType>::trackingData::trackingData
             cloud.solution().interpolationSchemes(),
             cloud.p()
         )
-    )
+    ),
+    pc_(Zero)
 {}
 
 
@@ -51,4 +52,18 @@ Foam::ReactingParcel<ParcelType>::trackingData::pInterp() const
 }
 
 
+template<class ParcelType>
+inline Foam::scalar Foam::ReactingParcel<ParcelType>::trackingData::pc() const
+{
+    return pc_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::ReactingParcel<ParcelType>::trackingData::pc()
+{
+    return pc_;
+}
+
+
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 006c45a9fdfcc3fe973de90d5159491e89ecf9db..f3da39feeccc533690cab22bf283aa7ffff5f2c5 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -35,19 +35,18 @@ template<class TrackCloudType>
 void Foam::ThermoParcel<ParcelType>::setCellValues
 (
     TrackCloudType& cloud,
-    trackingData& td,
-    const scalar dt
+    trackingData& td
 )
 {
-    ParcelType::setCellValues(cloud, td, dt);
+    ParcelType::setCellValues(cloud, td);
 
     tetIndices tetIs = this->currentTetIndices();
 
-    Cpc_ = td.CpInterp().interpolate(this->coordinates(), tetIs);
+    td.Cpc() = td.CpInterp().interpolate(this->coordinates(), tetIs);
 
-    Tc_ = td.TInterp().interpolate(this->coordinates(), tetIs);
+    td.Tc() = td.TInterp().interpolate(this->coordinates(), tetIs);
 
-    if (Tc_ < cloud.constProps().TMin())
+    if (td.Tc() < cloud.constProps().TMin())
     {
         if (debug)
         {
@@ -56,7 +55,7 @@ void Foam::ThermoParcel<ParcelType>::setCellValues
                 << " to " << cloud.constProps().TMin() <<  nl << endl;
         }
 
-        Tc_ = cloud.constProps().TMin();
+        td.Tc() = cloud.constProps().TMin();
     }
 }
 
@@ -71,7 +70,7 @@ void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
 )
 {
     const label celli = this->cell();
-    const scalar massCell = this->massCell(celli);
+    const scalar massCell = this->massCell(td);
 
     td.Uc() += cloud.UTrans()[celli]/massCell;
 
@@ -79,9 +78,9 @@ void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
 //    Tc_ = td.TInterp().interpolate(this->coordinates(), tetIs);
 
     const scalar CpMean = td.CpInterp().psi()[celli];
-    Tc_ += cloud.hsTrans()[celli]/(CpMean*massCell);
+    td.Tc() += cloud.hsTrans()[celli]/(CpMean*massCell);
 
-    if (Tc_ < cloud.constProps().TMin())
+    if (td.Tc() < cloud.constProps().TMin())
     {
         if (debug)
         {
@@ -90,7 +89,7 @@ void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
                 << " to " << cloud.constProps().TMin() <<  nl << endl;
         }
 
-        Tc_ = cloud.constProps().TMin();
+        td.Tc() = cloud.constProps().TMin();
     }
 }
 
@@ -110,7 +109,7 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
 ) const
 {
     // Surface temperature using two thirds rule
-    Ts = (2.0*T + Tc_)/3.0;
+    Ts = (2.0*T + td.Tc())/3.0;
 
     if (Ts < cloud.constProps().TMin())
     {
@@ -125,7 +124,7 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
     }
 
     // Assuming thermo props vary linearly with T for small d(T)
-    const scalar TRatio = Tc_/Ts;
+    const scalar TRatio = td.Tc()/Ts;
 
     rhos = td.rhoc()*TRatio;
 
@@ -133,7 +132,7 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
     mus = td.muInterp().interpolate(this->coordinates(), tetIs)/TRatio;
     kappas = td.kappaInterp().interpolate(this->coordinates(), tetIs)/TRatio;
 
-    Pr = Cpc_*mus/kappas;
+    Pr = td.Cpc()*mus/kappas;
     Pr = max(ROOTVSMALL, Pr);
 }
 
@@ -287,7 +286,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     htc = max(htc, ROOTVSMALL);
     const scalar As = this->areaS(d);
 
-    scalar ap = Tc_ + Sh/(As*htc);
+    scalar ap = td.Tc() + Sh/(As*htc);
     const scalar bp = 6.0*htc/max(rho*d*Cp_, ROOTVSMALL);
 
     if (cloud.radiation())
@@ -320,7 +319,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     Sph = dt*htc*As;
 
-    dhsTrans += Sph*(Tres.average() - Tc_);
+    dhsTrans += Sph*(Tres.average() - td.Tc());
 
     return Tnew;
 }
@@ -336,9 +335,7 @@ Foam::ThermoParcel<ParcelType>::ThermoParcel
 :
     ParcelType(p),
     T_(p.T_),
-    Cp_(p.Cp_),
-    Tc_(p.Tc_),
-    Cpc_(p.Cpc_)
+    Cp_(p.Cp_)
 {}
 
 
@@ -351,9 +348,7 @@ Foam::ThermoParcel<ParcelType>::ThermoParcel
 :
     ParcelType(p, mesh),
     T_(p.T_),
-    Cp_(p.Cp_),
-    Tc_(p.Tc_),
-    Cpc_(p.Cpc_)
+    Cp_(p.Cp_)
 {}
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index 8e97c35ab173fb5f2d6b13262c50fac491e71eb7..c45c3377e26831543549c6dc2397becc756e836f 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -179,6 +179,15 @@ public:
                 autoPtr<interpolation<scalar>> GInterp_;
 
 
+            // Cached continuous phase properties
+
+                //- Temperature [K]
+                scalar Tc_;
+
+                //- Specific heat capacity [J/(kg.K)]
+                scalar Cpc_;
+
+
     public:
 
         typedef typename ParcelType::trackingData::trackPart trackPart;
@@ -217,6 +226,18 @@ public:
             //- Return const access to the interpolator for continuous
             //  radiation field
             inline const interpolation<scalar>& GInterp() const;
+
+            //- Return the continuous phase temperature
+            inline scalar Tc() const;
+
+            //- Access the continuous phase temperature
+            inline scalar& Tc();
+
+            //- Return the continuous phase specific heat capacity
+            inline scalar Cpc() const;
+
+            //- Access the continuous phase specific heat capacity
+            inline scalar& Cpc();
     };
 
 
@@ -233,15 +254,6 @@ protected:
             scalar Cp_;
 
 
-        // Cell-based quantities
-
-            //- Temperature [K]
-            scalar Tc_;
-
-            //- Specific heat capacity [J/(kg.K)]
-            scalar Cpc_;
-
-
     // Protected Member Functions
 
         //- Calculate new particle temperature
@@ -387,12 +399,6 @@ public:
             //- Return the parcel sensible enthalpy
             inline scalar hs() const;
 
-            //- Return const access to carrier temperature
-            inline scalar Tc() const;
-
-            //- Return const access to carrier specific heat capacity
-            inline scalar Cpc() const;
-
 
         // Edit
 
@@ -407,12 +413,7 @@ public:
 
             //- Set cell values
             template<class TrackCloudType>
-            void setCellValues
-            (
-                TrackCloudType& cloud,
-                trackingData& td,
-                const scalar dt
-            );
+            void setCellValues(TrackCloudType& cloud, trackingData& td);
 
             //- Correct cell values using latest transfer information
             template<class TrackCloudType>
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
index 9c35a55a6f1edb22da92dc753d696f5ef9275eab..8513b33ea52a14bb39e07aa4a190ad82478d22b3 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
@@ -82,9 +82,7 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
 :
     ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
     T_(0.0),
-    Cp_(0.0),
-    Tc_(0.0),
-    Cpc_(0.0)
+    Cp_(0.0)
 {}
 
 
@@ -98,9 +96,7 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
 :
     ParcelType(mesh, position, celli),
     T_(0.0),
-    Cp_(0.0),
-    Tc_(0.0),
-    Cpc_(0.0)
+    Cp_(0.0)
 {}
 
 
@@ -141,9 +137,7 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
         constProps
     ),
     T_(constProps.T0()),
-    Cp_(constProps.Cp0()),
-    Tc_(0.0),
-    Cpc_(0.0)
+    Cp_(constProps.Cp0())
 {}
 
 
@@ -228,20 +222,6 @@ inline Foam::scalar Foam::ThermoParcel<ParcelType>::hs() const
 }
 
 
-template<class ParcelType>
-inline Foam::scalar Foam::ThermoParcel<ParcelType>::Tc() const
-{
-    return Tc_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar Foam::ThermoParcel<ParcelType>::Cpc() const
-{
-    return Cpc_;
-}
-
-
 template<class ParcelType>
 inline Foam::scalar& Foam::ThermoParcel<ParcelType>::T()
 {
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
index e9b1b64194925eec36283fc53687304168a98f64..0b658e033a0d1ffd815b6b2310c4fecbefaecb15 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
@@ -39,7 +39,7 @@ Foam::string Foam::ThermoParcel<ParcelType>::propertyTypes_ =
 template<class ParcelType>
 const std::size_t Foam::ThermoParcel<ParcelType>::sizeofFields
 (
-    offsetof(ThermoParcel<ParcelType>, Tc_)
+    sizeof(ThermoParcel<ParcelType>)
   - offsetof(ThermoParcel<ParcelType>, T_)
 );
 
@@ -56,9 +56,7 @@ Foam::ThermoParcel<ParcelType>::ThermoParcel
 :
     ParcelType(mesh, is, readFields),
     T_(0.0),
-    Cp_(0.0),
-    Tc_(0.0),
-    Cpc_(0.0)
+    Cp_(0.0)
 {
     if (readFields)
     {
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelTrackingDataI.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelTrackingDataI.H
index 73adcd177f59e5b27277b080cf5dfd2626a2d745..5ea8189bce95f0bde9fb9df3b4311142da617422 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelTrackingDataI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelTrackingDataI.H
@@ -58,7 +58,9 @@ inline Foam::ThermoParcel<ParcelType>::trackingData::trackingData
             kappa_
         )
     ),
-    GInterp_(nullptr)
+    GInterp_(nullptr),
+    Tc_(Zero),
+    Cpc_(Zero)
 {
     if (cloud.radiation())
     {
@@ -130,4 +132,32 @@ Foam::ThermoParcel<ParcelType>::trackingData::GInterp() const
 }
 
 
+template<class ParcelType>
+inline Foam::scalar Foam::ThermoParcel<ParcelType>::trackingData::Tc() const
+{
+    return Tc_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::ThermoParcel<ParcelType>::trackingData::Tc()
+{
+    return Tc_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::ThermoParcel<ParcelType>::trackingData::Cpc() const
+{
+    return Cpc_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::ThermoParcel<ParcelType>::trackingData::Cpc()
+{
+    return Cpc_;
+}
+
+
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/NoStochasticCollision/NoStochasticCollision.C b/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/NoStochasticCollision/NoStochasticCollision.C
index 9a62cf70e4b6f4ba96b80c11a22f8ffa8fbd530b..65eb4d2979f527edbe83836b1ca1e5febb8a8053 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/NoStochasticCollision/NoStochasticCollision.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/NoStochasticCollision/NoStochasticCollision.C
@@ -28,7 +28,11 @@ License
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
 template<class CloudType>
-void Foam::NoStochasticCollision<CloudType>::collide(const scalar)
+void Foam::NoStochasticCollision<CloudType>::collide
+(
+    typename CloudType::parcelType::trackingData&,
+    const scalar
+)
 {}
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/NoStochasticCollision/NoStochasticCollision.H b/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/NoStochasticCollision/NoStochasticCollision.H
index a2ad6ca6a1edf3873bb691aba018a105a37e338e..7e785fc7555e417a0745db667c00a981790b8f01 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/NoStochasticCollision/NoStochasticCollision.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/NoStochasticCollision/NoStochasticCollision.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
@@ -55,7 +55,11 @@ protected:
     // Protected Member Functions
 
         //- Update the model
-        virtual void collide(const scalar dt);
+        virtual void collide
+        (
+            typename CloudType::parcelType::trackingData& td,
+            const scalar dt
+        );
 
 
 public:
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/StochasticCollisionModel/StochasticCollisionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/StochasticCollisionModel/StochasticCollisionModel.C
index 0ed36fac1bb34672774e8722cac056e7e46a6ea1..a53ded90052ce3ae2414d36f1751fdd943be48da 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/StochasticCollisionModel/StochasticCollisionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/StochasticCollisionModel/StochasticCollisionModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,11 +69,15 @@ Foam::StochasticCollisionModel<CloudType>::~StochasticCollisionModel()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
-void Foam::StochasticCollisionModel<CloudType>::update(const scalar dt)
+void Foam::StochasticCollisionModel<CloudType>::update
+(
+    typename CloudType::parcelType::trackingData& td,
+    const scalar dt
+)
 {
     if (this->active())
     {
-        this->collide(dt);
+        this->collide(td, dt);
     }
 }
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/StochasticCollisionModel/StochasticCollisionModel.H b/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/StochasticCollisionModel/StochasticCollisionModel.H
index 9be2fc6bb5808d687dba800637c2b3c2317fca1b..2210b0b8c345c127f3b7b588f3ee1995f1af5358 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/StochasticCollisionModel/StochasticCollisionModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/StochasticCollision/StochasticCollisionModel/StochasticCollisionModel.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
@@ -61,7 +61,11 @@ class StochasticCollisionModel
 protected:
 
     //- Main collision routine
-    virtual void collide(const scalar dt) = 0;
+    virtual void collide
+    (
+        typename CloudType::parcelType::trackingData& td,
+        const scalar dt
+    ) = 0;
 
 
 public:
@@ -118,7 +122,11 @@ public:
     // Member Functions
 
         //- Update the model
-        void update(const scalar dt);
+        void update
+        (
+            typename CloudType::parcelType::trackingData& td,
+            const scalar dt
+        );
 };
 
 
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/StochasticCollision/SuppressionCollision/SuppressionCollision.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/StochasticCollision/SuppressionCollision/SuppressionCollision.C
index c7b49097b90467d722a6e5f101c255b93b6b0c8f..1bd646338c4cc89860c61f39a67e58fdc6471ed3 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/StochasticCollision/SuppressionCollision/SuppressionCollision.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/StochasticCollision/SuppressionCollision/SuppressionCollision.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
@@ -29,7 +29,11 @@ License
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
 template<class CloudType>
-void Foam::SuppressionCollision<CloudType>::collide(const scalar dt)
+void Foam::SuppressionCollision<CloudType>::collide
+(
+    typename CloudType::parcelType::trackingData& td,
+    const scalar dt
+)
 {
     const kinematicCloud& sc =
         this->owner().mesh().template
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/StochasticCollision/SuppressionCollision/SuppressionCollision.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/StochasticCollision/SuppressionCollision/SuppressionCollision.H
index 806b46af8c67811b6fc1bf69e192663363df22f5..81135fbb5972d5ea7688bf852c9904ceffa2245b 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/StochasticCollision/SuppressionCollision/SuppressionCollision.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/StochasticCollision/SuppressionCollision/SuppressionCollision.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
@@ -65,7 +65,11 @@ protected:
    // Protected Member Functions
 
         //- Update the model
-        virtual void collide(const scalar dt);
+        virtual void collide
+        (
+            typename CloudType::parcelType::trackingData& td,
+            const scalar dt
+        );
 
 
 public:
diff --git a/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.C b/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.C
index cfb76f21248639dabb7a52c9c8bd65fea019d6b3..7144ad4e510a0023363028095de28d0b911750ec 100644
--- a/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.C
+++ b/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.C
@@ -162,12 +162,13 @@ void Foam::SprayCloud<CloudType>::setParcelThermoProperties
 
     const scalarField& Y(parcel.Y());
     scalarField X(liqMix.X(Y));
+    const scalar pc = this->p()[parcel.cell()];
 
     // override rho and Cp from constantProperties
-    parcel.Cp() = liqMix.Cp(parcel.pc(), parcel.T(), X);
-    parcel.rho() = liqMix.rho(parcel.pc(), parcel.T(), X);
-    parcel.sigma() = liqMix.sigma(parcel.pc(), parcel.T(), X);
-    parcel.mu() = liqMix.mu(parcel.pc(), parcel.T(), X);
+    parcel.Cp() = liqMix.Cp(pc, parcel.T(), X);
+    parcel.rho() = liqMix.rho(pc, parcel.T(), X);
+    parcel.sigma() = liqMix.sigma(pc, parcel.T(), X);
+    parcel.mu() = liqMix.mu(pc, parcel.T(), X);
 }
 
 
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
index abdcac0d03102f5db5e0979fb3dfcd73b819ceae..87ca3fbff107dd9b5bbb202519b1b9f0c7bc1f3c 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
@@ -35,11 +35,10 @@ template<class TrackCloudType>
 void Foam::SprayParcel<ParcelType>::setCellValues
 (
     TrackCloudType& cloud,
-    trackingData& td,
-    const scalar dt
+    trackingData& td
 )
 {
-    ParcelType::setCellValues(cloud, td, dt);
+    ParcelType::setCellValues(cloud, td);
 }
 
 
@@ -81,7 +80,7 @@ void Foam::SprayParcel<ParcelType>::calc
     // Check if we have critical or boiling conditions
     scalar TMax = liquids.Tc(X0);
     const scalar T0 = this->T();
-    const scalar pc0 = this->pc_;
+    const scalar pc0 = td.pc();
     if (liquids.pv(pc0, T0, X0) >= pc0*0.999)
     {
         // Set TMax to boiling temperature
@@ -112,14 +111,14 @@ void Foam::SprayParcel<ParcelType>::calc
         scalar T1 = this->T();
         scalarField X1(liquids.X(this->Y()));
 
-        this->Cp() = liquids.Cp(this->pc_, T1, X1);
+        this->Cp() = liquids.Cp(td.pc(), T1, X1);
 
-        sigma_ = liquids.sigma(this->pc_, T1, X1);
+        sigma_ = liquids.sigma(td.pc(), T1, X1);
 
-        scalar rho1 = liquids.rho(this->pc_, T1, X1);
+        scalar rho1 = liquids.rho(td.pc(), T1, X1);
         this->rho() = rho1;
 
-        mu_ = liquids.mu(this->pc_, T1, X1);
+        mu_ = liquids.mu(td.pc(), T1, X1);
 
         scalar d1 = this->d()*cbrt(rho0/rho1);
         this->d() = d1;
@@ -164,12 +163,12 @@ void Foam::SprayParcel<ParcelType>::calcAtomization
     const auto& liquids = composition.liquids();
 
     // Average molecular weight of carrier mix - assumes perfect gas
-    scalar Wc = td.rhoc()*RR*this->Tc()/this->pc();
+    scalar Wc = td.rhoc()*RR*td.Tc()/td.pc();
     scalar R = RR/Wc;
-    scalar Tav = atomization.Taverage(this->T(), this->Tc());
+    scalar Tav = atomization.Taverage(this->T(), td.Tc());
 
     // Calculate average gas density based on average temperature
-    scalar rhoAv = this->pc()/(R*Tav);
+    scalar rhoAv = td.pc()/(R*Tav);
 
     scalar soi = cloud.injectors().timeStart();
     scalar currentTime = cloud.db().time().value();
@@ -235,12 +234,12 @@ void Foam::SprayParcel<ParcelType>::calcBreakup
     }
 
     // Average molecular weight of carrier mix - assumes perfect gas
-    scalar Wc = td.rhoc()*RR*this->Tc()/this->pc();
+    scalar Wc = td.rhoc()*RR*td.Tc()/td.pc();
     scalar R = RR/Wc;
-    scalar Tav = cloud.atomization().Taverage(this->T(), this->Tc());
+    scalar Tav = cloud.atomization().Taverage(this->T(), td.Tc());
 
     // Calculate average gas density based on average temperature
-    scalar rhoAv = this->pc()/(R*Tav);
+    scalar rhoAv = td.pc()/(R*Tav);
     scalar muAv = td.muc();
     vector Urel = this->U() - td.Uc();
     scalar Urmag = mag(Urel);
@@ -314,7 +313,7 @@ void Foam::SprayParcel<ParcelType>::calcBreakup
         child->injector() = this->injector();
         child->tMom() = massChild/(Fcp.Sp() + Fncp.Sp());
         child->user() = 0.0;
-        child->setCellValues(cloud, td, dt);
+        child->calcDispersion(cloud, td, dt);
 
         cloud.addParticle(child);
     }
@@ -337,7 +336,7 @@ Foam::scalar Foam::SprayParcel<ParcelType>::chi
 
     scalar chi = 0.0;
     scalar T0 = this->T();
-    scalar p0 = this->pc();
+    scalar p0 = td.pc();
     scalar pAmb = cloud.pAmbient();
 
     scalar pv = liquids.pv(p0, T0, X);
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
index 11b617af229688b2186b446d4a0c0ab39385fc27..0fb14b2cf4e867315942aedb952eb641f45f5d29 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
@@ -418,12 +418,7 @@ public:
 
             //- Set cell values
             template<class TrackCloudType>
-            void setCellValues
-            (
-                TrackCloudType& cloud,
-                trackingData& td,
-                const scalar dt
-            );
+            void setCellValues(TrackCloudType& cloud, trackingData& td);
 
             //- Correct parcel properties according to atomization model
             template<class TrackCloudType>
diff --git a/src/lagrangian/spray/submodels/StochasticCollision/ORourkeCollision/ORourkeCollision.C b/src/lagrangian/spray/submodels/StochasticCollision/ORourkeCollision/ORourkeCollision.C
index 686e39093e02ed8f5e6d6db4df708a6df65005e0..376953ba90f496a6b0415b321e8639709f76772e 100644
--- a/src/lagrangian/spray/submodels/StochasticCollision/ORourkeCollision/ORourkeCollision.C
+++ b/src/lagrangian/spray/submodels/StochasticCollision/ORourkeCollision/ORourkeCollision.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
@@ -33,7 +33,11 @@ using namespace Foam::constant::mathematical;
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
 template<class CloudType>
-void Foam::ORourkeCollision<CloudType>::collide(const scalar dt)
+void Foam::ORourkeCollision<CloudType>::collide
+(
+    typename CloudType::parcelType::trackingData& td,
+    const scalar dt
+)
 {
     // Create the occupancy list for the cells
     labelList occupancy(this->owner().mesh().nCells(), 0);
@@ -77,20 +81,22 @@ void Foam::ORourkeCollision<CloudType>::collide(const scalar dt)
                         if (m1 > ROOTVSMALL)
                         {
                             const scalarField X(liquids_.X(p1.Y()));
-                            p1.rho() = liquids_.rho(p1.pc(), p1.T(), X);
-                            p1.Cp() = liquids_.Cp(p1.pc(), p1.T(), X);
-                            p1.sigma() = liquids_.sigma(p1.pc(), p1.T(), X);
-                            p1.mu() = liquids_.mu(p1.pc(), p1.T(), X);
+                            p1.setCellValues(this->owner(), td);
+                            p1.rho() = liquids_.rho(td.pc(), p1.T(), X);
+                            p1.Cp() = liquids_.Cp(td.pc(), p1.T(), X);
+                            p1.sigma() = liquids_.sigma(td.pc(), p1.T(), X);
+                            p1.mu() = liquids_.mu(td.pc(), p1.T(), X);
                             p1.d() = cbrt(6.0*m1/(p1.nParticle()*p1.rho()*pi));
                         }
 
                         if (m2 > ROOTVSMALL)
                         {
                             const scalarField X(liquids_.X(p2.Y()));
-                            p2.rho() = liquids_.rho(p2.pc(), p2.T(), X);
-                            p2.Cp() = liquids_.Cp(p2.pc(), p2.T(), X);
-                            p2.sigma() = liquids_.sigma(p2.pc(), p2.T(), X);
-                            p2.mu() = liquids_.mu(p2.pc(), p2.T(), X);
+                            p2.setCellValues(this->owner(), td);
+                            p2.rho() = liquids_.rho(td.pc(), p2.T(), X);
+                            p2.Cp() = liquids_.Cp(td.pc(), p2.T(), X);
+                            p2.sigma() = liquids_.sigma(td.pc(), p2.T(), X);
+                            p2.mu() = liquids_.mu(td.pc(), p2.T(), X);
                             p2.d() = cbrt(6.0*m2/(p2.nParticle()*p2.rho()*pi));
                         }
                     }
diff --git a/src/lagrangian/spray/submodels/StochasticCollision/ORourkeCollision/ORourkeCollision.H b/src/lagrangian/spray/submodels/StochasticCollision/ORourkeCollision/ORourkeCollision.H
index aa26c45341512fb3db9c18ef83ba26a11d2063a2..7972ac83e66957bb1a4c6fa37226166a14b3da27 100644
--- a/src/lagrangian/spray/submodels/StochasticCollision/ORourkeCollision/ORourkeCollision.H
+++ b/src/lagrangian/spray/submodels/StochasticCollision/ORourkeCollision/ORourkeCollision.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
@@ -65,7 +65,11 @@ protected:
     // Protected Member Functions
 
         //- Main collision routine
-        virtual void collide(const scalar dt);
+        virtual void collide
+        (
+            typename CloudType::parcelType::trackingData& td,
+            const scalar dt
+        );
 
         //- Collide parcels and return true if mass has changed
         virtual bool collideParcels
diff --git a/src/lagrangian/spray/submodels/StochasticCollision/TrajectoryCollision/TrajectoryCollision.C b/src/lagrangian/spray/submodels/StochasticCollision/TrajectoryCollision/TrajectoryCollision.C
index 4abc032432b2683f26988af4f2f3ce40f85e84c3..2d608941415234a7b50df0a8ab6cc473e28100ad 100644
--- a/src/lagrangian/spray/submodels/StochasticCollision/TrajectoryCollision/TrajectoryCollision.C
+++ b/src/lagrangian/spray/submodels/StochasticCollision/TrajectoryCollision/TrajectoryCollision.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,9 +28,13 @@ License
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
 template<class CloudType>
-void Foam::TrajectoryCollision<CloudType>::collide(const scalar dt)
+void Foam::TrajectoryCollision<CloudType>::collide
+(
+    typename CloudType::parcelType::trackingData& td,
+    const scalar dt
+)
 {
-    ORourkeCollision<CloudType>::collide(dt);
+    ORourkeCollision<CloudType>::collide(td, dt);
 }
 
 
diff --git a/src/lagrangian/spray/submodels/StochasticCollision/TrajectoryCollision/TrajectoryCollision.H b/src/lagrangian/spray/submodels/StochasticCollision/TrajectoryCollision/TrajectoryCollision.H
index 722ebd1db0f6a3a8b16bd61c818ae8ac43a55bb0..91c53dfca0fa2c8d4f90857552acdb948f8ea1f0 100644
--- a/src/lagrangian/spray/submodels/StochasticCollision/TrajectoryCollision/TrajectoryCollision.H
+++ b/src/lagrangian/spray/submodels/StochasticCollision/TrajectoryCollision/TrajectoryCollision.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
@@ -65,7 +65,11 @@ protected:
     // Protected Member Functions
 
         //- Main collision routine
-        virtual void collide(const scalar dt);
+        virtual void collide
+        (
+            typename CloudType::parcelType::trackingData& td,
+            const scalar dt
+        );
 
         //- Collide parcels and return true if mass has changed
         virtual bool collideParcels
diff --git a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
index 695427b9d4209ae9abe6e1023346b22bdd229789..a771a0837d7202bb2cdeaabe36261817c97d1ff2 100644
--- a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
+++ b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.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
@@ -170,7 +170,7 @@ Foam::forceSuSp Foam::BrownianMotionForce<CloudType>::calcCoupled
     forceSuSp value(Zero, 0.0);
 
     const scalar dp = p.d();
-    const scalar Tc = p.Tc();
+    const scalar Tc = td.Tc();
 
     const scalar alpha = 2.0*lambda_/dp;
     const scalar cc = 1.0 + alpha*(1.257 + 0.4*exp(-1.1/alpha));
diff --git a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
index 611b64407fe627226935734608ddaaf157c6d87a..8ae7134aea443e726ebb86684aeafd51f4b270e9 100644
--- a/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
+++ b/src/lagrangian/turbulence/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.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