From a0d36e3f90a88a1f30b0eea86594a883c34ce89b Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Fri, 27 Mar 2009 13:15:06 +0000
Subject: [PATCH] further dev on coupling, bits of restructuring

---
 .../Templates/KinematicCloud/KinematicCloud.C |   4 +-
 .../Templates/ReactingCloud/ReactingCloud.C   |   4 +-
 .../ReactingMultiphaseCloud.C                 |   6 +-
 .../Templates/ThermoCloud/ThermoCloud.C       |   4 +-
 .../KinematicParcel/KinematicParcel.C         |  30 ++-
 .../KinematicParcel/KinematicParcel.H         |  77 ++++---
 .../KinematicParcel/KinematicParcelI.H        | 195 ++++++++++-------
 .../KinematicParcel/KinematicParcelIO.C       |  27 +--
 .../ReactingMultiphaseParcel.C                | 207 ++++++++++++------
 .../ReactingMultiphaseParcel.H                |  56 +++--
 .../ReactingMultiphaseParcelI.H               |  51 ++---
 .../ReactingMultiphaseParcelIO.C              |   8 +-
 .../Templates/ReactingParcel/ReactingParcel.C |  62 ++++--
 .../Templates/ReactingParcel/ReactingParcel.H |  28 +--
 .../ReactingParcel/ReactingParcelI.H          |  12 +-
 .../ReactingParcel/ReactingParcelIO.C         |   6 +-
 .../Templates/ThermoParcel/ThermoParcel.C     |  58 +++--
 .../Templates/ThermoParcel/ThermoParcel.H     |  19 +-
 .../Templates/ThermoParcel/ThermoParcelI.H    |  16 +-
 .../Templates/ThermoParcel/ThermoParcelIO.C   |   4 +-
 .../basicKinematicParcel.C                    |   8 +-
 .../basicKinematicParcel.H                    |   9 +-
 .../basicReactingMultiphaseParcel.C           |   8 +-
 .../basicReactingMultiphaseParcel.H           |   6 +-
 .../basicReactingParcel/basicReactingParcel.C |   8 +-
 .../basicReactingParcel/basicReactingParcel.H |   4 +-
 .../basicThermoParcel/basicThermoParcel.C     |   8 +-
 .../basicThermoParcel/basicThermoParcel.H     |  11 +-
 .../IO/DataEntry/polynomial/polynomial.H      |   2 +-
 .../CompositionModel/CompositionModel.C       |   2 +-
 .../CompositionModel/CompositionModel.H       |   2 +-
 .../LiquidEvaporation/LiquidEvaporation.C     |   6 +-
 .../LiquidEvaporation/LiquidEvaporation.H     |   6 +-
 .../NoPhaseChange/NoPhaseChange.C             |   2 +-
 .../NoPhaseChange/NoPhaseChange.H             |   4 +-
 .../PhaseChangeModel/PhaseChangeModel.H       |   4 +-
 .../NoSurfaceReaction/NoSurfaceReaction.C     |   6 +-
 .../NoSurfaceReaction/NoSurfaceReaction.H     |   6 +-
 .../SurfaceReactionModel.H                    |   7 +-
 39 files changed, 592 insertions(+), 391 deletions(-)

diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index f60f69261c..eeee670009 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -50,12 +50,12 @@ void Foam::KinematicCloud<ParcelType>::addNewParcel
     ParcelType* pPtr = new ParcelType
     (
         *this,
-        parcelTypeId_,
         position,
         cellId,
+        parcelTypeId_,
+        nParticles,
         d,
         U,
-        nParticles,
         constProps_
     );
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index b14fcc17cd..4fec0e0666 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -45,12 +45,12 @@ void Foam::ReactingCloud<ParcelType>::addNewParcel
     ParcelType* pPtr = new ParcelType
     (
         *this,
-        this->parcelTypeId(),
         position,
         cellId,
+        this->parcelTypeId(),
+        nParticles,
         d,
         U,
-        nParticles,
         composition().YMixture0(),
         constProps_
     );
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
index 29936f328c..7936ddea40 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
@@ -49,16 +49,16 @@ void Foam::ReactingMultiphaseCloud<ParcelType>::addNewParcel
     ParcelType* pPtr = new ParcelType
     (
         *this,
-        this->parcelTypeId(),
         position,
         cellId,
+        this->parcelTypeId(),
+        nParticles,
         d,
         U,
-        nParticles,
+        this->composition().YMixture0(),
         this->composition().Y0(idGas),
         this->composition().Y0(idLiquid),
         this->composition().Y0(idSolid),
-        this->composition().YMixture0(),
         constProps_
     );
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index 40731ec254..6872f2428e 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -46,12 +46,12 @@ void Foam::ThermoCloud<ParcelType>::addNewParcel
     ParcelType* pPtr = new ParcelType
     (
         *this,
-        this->parcelTypeId(),
         position,
         cellId,
+        this->parcelTypeId(),
+        nParticles,
         d,
         U,
-        nParticles,
         constProps_
     );
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 70adc771b1..007e5aac47 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -79,6 +79,10 @@ void Foam::KinematicParcel<ParcelType>::calc
 {
     // Define local properties at beginning of time step
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    const scalar np0 = nParticle_;
+    const scalar d0 = d_;
+    const scalar U0 = U_;
+    const scalar rho0 = rho_;
     const scalar mass0 = mass();
 
 
@@ -97,7 +101,8 @@ void Foam::KinematicParcel<ParcelType>::calc
 
     // Calculate new particle velocity
     scalar Cud = 0.0;
-    vector U1 = calcVelocity(td, dt, cellI, Fx, mass0, Cud, dUTrans);
+    vector U1 =
+        calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx, Cud, dUTrans);
 
 
     // Accumulate carrier phase source terms
@@ -105,10 +110,10 @@ void Foam::KinematicParcel<ParcelType>::calc
     if (td.cloud().coupled())
     {
         // Update momentum transfer
-        td.cloud().UTrans()[cellI] += nParticle_*dUTrans;
+        td.cloud().UTrans()[cellI] += np0*dUTrans;
 
         // Coefficient to be applied in carrier phase momentum coupling
-        td.cloud().UCoeff()[cellI] += nParticle_*mass0*Cud;
+        td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
     }
 
 
@@ -125,14 +130,17 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     TrackData& td,
     const scalar dt,
     const label cellI,
-    const vector& Fx,
+    const scalar d,
+    const vector& U,
+    const scalar rho,
     const scalar mass,
+    const vector& Fx,
     scalar& Cud,
     vector& dUTrans
 ) const
 {
     // Return linearised term from drag model
-    Cud = td.cloud().drag().Cu(U_ - Uc_, d_, rhoc_, rho_, muc_);
+    Cud = td.cloud().drag().Cu(U - Uc_, d, rhoc_, rho, muc_);
 
     // Initialise total force (per unit mass)
     vector Ftot = vector::zero;
@@ -140,20 +148,20 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     // Gravity force
     if (td.cloud().forceGravity())
     {
-        Ftot += td.g()*(1 - rhoc_/rho_);
+        Ftot += td.g()*(1 - rhoc_/rho);
     }
 
     // Virtual mass force
     if (td.cloud().forceVirtualMass())
     {
-//        Ftot += td.constProps().Cvm()*rhoc_/rho_*d(Uc - U_)/dt;
+//        Ftot += td.constProps().Cvm()*rhoc_/rho*d(Uc - U)/dt;
     }
 
     // Pressure gradient force
     if (td.cloud().forcePressureGradient())
     {
-        const vector& d = this->mesh().deltaCoeffs()[cellI];
-        Ftot += rhoc_/rho_*(U_ & (d^Uc_));
+//        const vector& delta = td.cloud().mesh().deltaCoeffs()[cellI];
+//        Ftot += rhoc_/rho*(U & (delta^Uc_));
     }
 
 
@@ -164,11 +172,11 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     const vector ap = Uc_ + (Ftot + Fx)/(Cud + VSMALL);
     const scalar bp = Cud;
 
-    vector Unew = td.cloud().UIntegrator().integrate(U_, dt, ap, bp).value();
+    vector Unew = td.cloud().UIntegrator().integrate(U, dt, ap, bp).value();
 
     // Calculate the momentum transfer to the continuous phase
     // - do not include gravity impulse
-    dUTrans = -mass*(Unew - U_ - dt*td.g());
+    dUTrans = -mass*(Unew - U - dt*td.g());
 
     return Unew;
 }
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index 35b46329dc..1717ee1feb 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -31,7 +31,7 @@ Description
 
     Sub-models include:
     - drag
-    - break-up
+    - turbulent dispersion
     - wall interactions
 
 SourceFiles
@@ -193,15 +193,15 @@ protected:
             //- Parcel type id
             label typeId_;
 
+            //- Number of particles in Parcel
+            scalar nParticle_;
+
             //- Diameter [m]
             scalar d_;
 
             //- Velocity of Parcel [m/s]
             vector U_;
 
-            //- Number of particles in Parcel
-            scalar nParticle_;
-
             //- Density [kg/m3]
             scalar rho_;
 
@@ -214,13 +214,13 @@ protected:
 
         // Cell-based quantities
 
-            // - Density [kg/m3]
+            //- Density [kg/m3]
             scalar rhoc_;
 
-            // - Velocity [m/s]
+            //- Velocity [m/s]
             vector Uc_;
 
-            // - Viscosity [Pa.s]
+            //- Viscosity [Pa.s]
             scalar muc_;
 
 
@@ -231,12 +231,15 @@ protected:
         const vector calcVelocity
         (
             TrackData& td,
-            const scalar dt,
-            const label cellI,
-            const vector& Fx,
-            const scalar mass,
-            scalar& Cud,
-            vector& dUTrans
+            const scalar dt,           // timestep
+            const label cellI,         // owner cell
+            const scalar d,            // diameter
+            const vector& U,           // velocity
+            const scalar rho,          // density
+            const scalar mass,         // mass
+            const vector& Fx,          // additional forces
+            scalar& Cud,               // linearised drag term coeff
+            vector& dUTrans            // momentum transfer to carrier phase
         ) const;
 
 
@@ -254,12 +257,12 @@ public:
         inline KinematicParcel
         (
             KinematicCloud<ParcelType>& owner,
-            const label typeId,
             const vector& position,
             const label cellI,
+            const label typeId,
+            const scalar nParticle0,
             const scalar d0,
             const vector& U0,
-            const scalar nParticle0,
             const constantProperties& constProps
         );
 
@@ -282,21 +285,18 @@ public:
 
         // Access
 
-            //- Return type id
+            //- Return const access to type id
             inline label typeId() const;
 
+            //- Return const access to number of particles
+            inline scalar nParticle() const;
+
             //- Return const access to diameter
             inline scalar d() const;
 
             //- Return const access to velocity
             inline const vector& U() const;
 
-            //- Return const access to relative velocity
-            inline const vector& Ur() const;
-
-            //- Return const access to number of particles
-            inline scalar nParticle() const;
-
             //- Return const access to density
             inline scalar rho() const;
 
@@ -306,25 +306,21 @@ public:
             //- Return const access to turbulent velocity fluctuation
             inline const vector& UTurb() const;
 
-            //- The nearest distance to a wall that
-            //  the particle can be in the n direction
-            inline scalar wallImpactDistance(const vector& n) const;
-
 
         // Edit
 
+            //- Return access to type id
+            inline label typeId();
+
+            //- Return access to number of particles
+            inline scalar& nParticle();
+
             //- Return access to diameter
             inline scalar& d();
 
             //- Return access to velocity
             inline vector& U();
 
-            //- Return access to relative velocity
-            inline vector& Ur();
-
-            //- Return access to number of particles
-            inline scalar& nParticle();
-
             //- Return access to density
             inline scalar& rho();
 
@@ -337,22 +333,35 @@ public:
 
         // Helper functions
 
+            //- The nearest distance to a wall that
+            //  the particle can be in the n direction
+            inline scalar wallImpactDistance(const vector& n) const;
+
             //- Return the index of the face to be used in the interpolation
             //  routine
             inline label faceInterpolation() const;
 
+            //- Particle mass
+            inline scalar mass() const;
+
             //- Particle volume
             inline scalar volume() const;
 
-            //- Particle mass
-            inline scalar mass() const;
+            //- Particle volume for a given diameter
+            inline scalar volume(const scalar d) const;
 
             //- Particle projected area
             inline scalar areaP() const;
 
+            //- Projected area for given diameter
+            inline scalar areaP(const scalar d) const;
+
             //- Particle surface area
             inline scalar areaS() const;
 
+            //- Surface area for given diameter
+            inline scalar areaS(const scalar d) const;
+
 
         // Main calculation loop
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index 4a6be04a97..8582792460 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -63,20 +63,20 @@ template <class ParcelType>
 inline Foam::KinematicParcel<ParcelType>::KinematicParcel
 (
     KinematicCloud<ParcelType>& owner,
-    const label typeId,
     const vector& position,
     const label cellI,
+    const label typeId,
+    const scalar nParticle0,
     const scalar d0,
     const vector& U0,
-    const scalar nParticle0,
     const constantProperties& constProps
 )
 :
     Particle<ParcelType>(owner, position, cellI),
     typeId_(typeId),
+    nParticle_(nParticle0),
     d_(d0),
     U_(U0),
-    nParticle_(nParticle0),
     rho_(constProps.rho0()),
     tTurb_(0.0),
     UTurb_(vector::zero),
@@ -171,6 +171,13 @@ inline Foam::label Foam::KinematicParcel<ParcelType>::typeId() const
 }
 
 
+template <class ParcelType>
+inline Foam::scalar Foam::KinematicParcel<ParcelType>::nParticle() const
+{
+    return nParticle_;
+}
+
+
 template <class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::d() const
 {
@@ -178,6 +185,48 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::d() const
 }
 
 
+template <class ParcelType>
+inline const Foam::vector& Foam::KinematicParcel<ParcelType>::U() const
+{
+    return U_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar Foam::KinematicParcel<ParcelType>::rho() const
+{
+    return rho_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar Foam::KinematicParcel<ParcelType>::tTurb() const
+{
+    return tTurb_;
+}
+
+
+template <class ParcelType>
+inline const Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb() const
+{
+    return UTurb_;
+}
+
+
+template <class ParcelType>
+inline Foam::label Foam::KinematicParcel<ParcelType>::typeId()
+{
+    return typeId_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar& Foam::KinematicParcel<ParcelType>::nParticle()
+{
+    return nParticle_;
+}
+
+
 template <class ParcelType>
 inline Foam::scalar& Foam::KinematicParcel<ParcelType>::d()
 {
@@ -185,6 +234,34 @@ inline Foam::scalar& Foam::KinematicParcel<ParcelType>::d()
 }
 
 
+template <class ParcelType>
+inline Foam::vector& Foam::KinematicParcel<ParcelType>::U()
+{
+    return U_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar& Foam::KinematicParcel<ParcelType>::rho()
+{
+    return rho_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar& Foam::KinematicParcel<ParcelType>::tTurb()
+{
+    return tTurb_;
+}
+
+
+template <class ParcelType>
+inline Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb()
+{
+    return UTurb_;
+}
+
+
 template <class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::wallImpactDistance
 (
@@ -210,83 +287,6 @@ inline Foam::label Foam::KinematicParcel<ParcelType>::faceInterpolation() const
 }
 
 
-template <class ParcelType>
-inline const Foam::vector& Foam::KinematicParcel<ParcelType>::U() const
-{
-    return U_;
-}
-
-
-template <class ParcelType>
-inline Foam::vector& Foam::KinematicParcel<ParcelType>::U()
-{
-    return U_;
-}
-
-
-template <class ParcelType>
-inline Foam::scalar Foam::KinematicParcel<ParcelType>::nParticle() const
-{
-    return nParticle_;
-}
-
-
-template <class ParcelType>
-inline Foam::scalar& Foam::KinematicParcel<ParcelType>::nParticle()
-{
-    return nParticle_;
-}
-
-
-template <class ParcelType>
-inline Foam::scalar Foam::KinematicParcel<ParcelType>::rho() const
-{
-    return rho_;
-}
-
-
-template <class ParcelType>
-inline Foam::scalar& Foam::KinematicParcel<ParcelType>::rho()
-{
-    return rho_;
-}
-
-
-template <class ParcelType>
-inline Foam::scalar Foam::KinematicParcel<ParcelType>::tTurb() const
-{
-    return tTurb_;
-}
-
-
-template <class ParcelType>
-inline Foam::scalar& Foam::KinematicParcel<ParcelType>::tTurb()
-{
-    return tTurb_;
-}
-
-
-template <class ParcelType>
-inline const Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb() const
-{
-    return UTurb_;
-}
-
-
-template <class ParcelType>
-inline Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb()
-{
-    return UTurb_;
-}
-
-
-template <class ParcelType>
-inline Foam::scalar Foam::KinematicParcel<ParcelType>::volume() const
-{
-    return mathematicalConstant::pi/6.0*pow3(d_);
-}
-
-
 template <class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::mass() const
 {
@@ -294,17 +294,48 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::mass() const
 }
 
 
+template <class ParcelType>
+inline Foam::scalar Foam::KinematicParcel<ParcelType>::volume() const
+{
+    return volume(d_);
+}
+
+
+template <class ParcelType>
+inline Foam::scalar
+Foam::KinematicParcel<ParcelType>::volume(const scalar d) const
+{
+    return mathematicalConstant::pi/6.0*pow3(d);
+}
+
+
 template <class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::areaP() const
 {
-    return 0.25*areaS();
+    return areaP(d_);
+}
+
+
+template <class ParcelType>
+inline Foam::scalar
+Foam::KinematicParcel<ParcelType>::areaP(const scalar d) const
+{
+    return 0.25*areaS(d);
 }
 
 
 template <class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::areaS() const
 {
-    return mathematicalConstant::pi*d_*d_;
+    return areaS(d_);
+}
+
+
+template <class ParcelType>
+inline Foam::scalar
+Foam::KinematicParcel<ParcelType>::areaS(const scalar d) const
+{
+    return mathematicalConstant::pi*d*d;
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
index 36535ab6a0..ae96d25a14 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
@@ -41,9 +41,9 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
 :
     Particle<ParcelType>(cloud, is, readFields),
     typeId_(0),
+    nParticle_(0.0),
     d_(0.0),
     U_(vector::zero),
-    nParticle_(0.0),
     rho_(0.0),
     tTurb_(0.0),
     UTurb_(vector::zero),
@@ -56,9 +56,9 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
         if (is.format() == IOstream::ASCII)
         {
             typeId_ = readLabel(is);
+            nParticle_ = readScalar(is);
             d_ = readScalar(is);
             is >> U_;
-            nParticle_ = readScalar(is);
             rho_ = readScalar(is);
             tTurb_ = readScalar(is);
             is >> UTurb_;
@@ -69,9 +69,9 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
             (
                 reinterpret_cast<char*>(&typeId_),
                 sizeof(typeId_)
+              + sizeof(nParticle_)
               + sizeof(d_)
               + sizeof(U_)
-              + sizeof(nParticle_)
               + sizeof(rho_)
               + sizeof(tTurb_)
               + sizeof(UTurb_)
@@ -108,7 +108,8 @@ void Foam::KinematicParcel<ParcelType>::readFields
     IOField<vector> U(c.fieldIOobject("U", IOobject::MUST_READ));
     c.checkFieldIOobject(c, U);
 
-    IOField<scalar> nParticle(c.fieldIOobject("nParticle", IOobject::MUST_READ));
+    IOField<scalar>
+        nParticle(c.fieldIOobject("nParticle", IOobject::MUST_READ));
     c.checkFieldIOobject(c, nParticle);
 
     IOField<scalar> rho(c.fieldIOobject("rho", IOobject::MUST_READ));
@@ -126,9 +127,9 @@ void Foam::KinematicParcel<ParcelType>::readFields
         ParcelType& p = iter();
 
         p.typeId_ = typeId[i];
+        p.nParticle_ = nParticle[i];
         p.d_ = d[i];
         p.U_ = U[i];
-        p.nParticle_ = nParticle[i];
         p.rho_ = rho[i];
         p.tTurb_ = tTurb[i];
         p.UTurb_ = UTurb[i];
@@ -148,13 +149,13 @@ void Foam::KinematicParcel<ParcelType>::writeFields
     label np =  c.size();
 
     IOField<label> typeId(c.fieldIOobject("typeId", IOobject::NO_READ), np);
-    IOField<scalar> d(c.fieldIOobject("d", IOobject::NO_READ), np);
-    IOField<vector> U(c.fieldIOobject("U", IOobject::NO_READ), np);
     IOField<scalar> nParticle
     (
         c.fieldIOobject("nParticle", IOobject::NO_READ),
         np
     );
+    IOField<scalar> d(c.fieldIOobject("d", IOobject::NO_READ), np);
+    IOField<vector> U(c.fieldIOobject("U", IOobject::NO_READ), np);
     IOField<scalar> rho(c.fieldIOobject("rho", IOobject::NO_READ), np);
     IOField<scalar> tTurb(c.fieldIOobject("tTurb", IOobject::NO_READ), np);
     IOField<vector> UTurb(c.fieldIOobject("UTurb", IOobject::NO_READ), np);
@@ -165,9 +166,9 @@ void Foam::KinematicParcel<ParcelType>::writeFields
         const KinematicParcel<ParcelType>& p = iter();
 
         typeId[i] = p.typeId();
+        nParticle[i] = p.nParticle();
         d[i] = p.d();
         U[i] = p.U();
-        nParticle[i] = p.nParticle();
         rho[i] = p.rho();
         tTurb[i] = p.tTurb();
         UTurb[i] = p.UTurb();
@@ -175,9 +176,9 @@ void Foam::KinematicParcel<ParcelType>::writeFields
     }
 
     typeId.write();
+    nParticle.write();
     d.write();
     U.write();
-    nParticle.write();
     rho.write();
     tTurb.write();
     UTurb.write();
@@ -195,25 +196,25 @@ Foam::Ostream& Foam::operator<<
 {
     if (os.format() == IOstream::ASCII)
     {
-        os  << static_cast<const Particle<ParcelType>& >(p)
+        os  << static_cast<const Particle<ParcelType>&>(p)
             << token::SPACE << p.typeId()
+            << token::SPACE << p.nParticle()
             << token::SPACE << p.d()
             << token::SPACE << p.U()
-            << token::SPACE << p.nParticle()
             << token::SPACE << p.rho()
             << token::SPACE << p.tTurb()
             << token::SPACE << p.UTurb();
     }
     else
     {
-        os  << static_cast<const Particle<ParcelType>& >(p);
+        os  << static_cast<const Particle<ParcelType>&>(p);
         os.write
         (
             reinterpret_cast<const char*>(&p.typeId_),
             sizeof(p.typeId())
+          + sizeof(p.nParticle())
           + sizeof(p.d())
           + sizeof(p.U())
-          + sizeof(p.nParticle())
           + sizeof(p.rho())
           + sizeof(p.tTurb())
           + sizeof(p.UTurb())
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index aacd51f04a..70d04a910e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -87,7 +87,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
     const scalarField& dMassSolid
 )
 {
-    scalarList& YMix = this->Y_;
+    scalarField& YMix = this->Y_;
 
     scalar dMassGasTot = sum(dMassGas);
     scalar dMassLiquidTot = sum(dMassLiquid);
@@ -99,8 +99,8 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
 
     scalar massNew = mass0 - (dMassGasTot + dMassLiquidTot + dMassSolidTot);
 
-    YMix[GAS] = (mass0*YMix[GAS] - dMassGas)/massNew;
-    YMix[LIQUID] = (mass0*YMix[LIQUID] - dMassLiquid)/massNew;
+    YMix[GAS] = (mass0*YMix[GAS] - dMassGasTot)/massNew;
+    YMix[LIQUID] = (mass0*YMix[LIQUID] - dMassLiquidTot)/massNew;
     YMix[SOLID] = 1.0 - YMix[GAS] - YMix[LIQUID];
 
     return massNew;
@@ -133,11 +133,17 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 {
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const scalar mass0 = this->mass();
     const scalar np0 = this->nParticle_;
+    const scalar d0 = this->d_;
+    const vector& U0 = this->U_;
+    const scalar rho0 = this->rho_;
     const scalar T0 = this->T_;
+    const scalar cp0 = this->cp_;
+    const scalar mass0 = this->mass();
+
     const scalar pc = this->pc_;
-    scalarField& YMix = this->Y_;
+
+    const scalarField& YMix = this->Y_;
     const label idG = td.cloud().composition().idGas();
     const label idL = td.cloud().composition().idLiquid();
     const label idS = td.cloud().composition().idSolid();
@@ -170,14 +176,39 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 
     // Return enthalpy source and calc mass transfer due to phase change
     scalar ShPC =
-        calcPhaseChange(td, dt, cellI, T0, idL, YMix[idL], YLiquid_, dMassPC);
+        calcPhaseChange
+        (
+            td,
+            dt,
+            cellI,
+            d0,
+            T0,
+            U0,
+            idL,
+            YMix[LIQUID],
+            YLiquid_,
+            dMassPC
+        );
 
 
     // Devolatilisation
     // ~~~~~~~~~~~~~~~~
 
     // Return enthalpy source and calc mass transfer due to devolatilisation
-    scalar ShDV = calcDevolatilisation(td, dt, T0, mass0, idG, YMix, dMassDV);
+    scalar ShDV =
+        calcDevolatilisation
+        (
+            td,
+            dt,
+            T0,
+            mass0,
+            this->mass0_,
+            idG,
+            YMix[GAS],
+            YGas_,
+            canCombust_,
+            dMassDV
+        );
 
 
     // Surface reactions
@@ -190,8 +221,15 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             td,
             dt,
             cellI,
+            d0,
             T0,
-            dMassPC + dMassDV, // total volatiles evolved
+            mass0,
+            canCombust_,
+            dMassPC + dMassDV, // total mass of volatiles evolved
+            YMix,
+            YGas_,
+            YLiquid_,
+            YSolid_,
             dMassSRGas,
             dMassSRLiquid,
             dMassSRSolid,
@@ -208,7 +246,21 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 
     // Calculate new particle temperature
     scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, dt, cellI, Sh, htc, dhTrans);
+    scalar T1 =
+        calcHeatTransfer
+        (
+            td,
+            dt,
+            cellI,
+            d0,
+            U0,
+            rho0,
+            T0,
+            cp0,
+            Sh,
+            htc,
+            dhTrans
+        );
 
 
     // Update component mass fractions
@@ -231,7 +283,19 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // Calculate new particle velocity
     scalar Cud = 0.0;
     vector U1 =
-        calcVelocity(td, dt, cellI, Fx, 0.5*(mass0 + mass1), Cud, dUTrans);
+        calcVelocity
+        (
+            td,
+            dt,
+            cellI,
+            d0,
+            U0,
+            rho0,
+            0.5*(mass0 + mass1),
+            Fx,
+            Cud,
+            dUTrans
+        );
 
 
     // Accumulate carrier phase source terms
@@ -242,18 +306,18 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         // Transfer mass lost from particle to carrier mass source
         forAll(YGas_, i)
         {
-            label id = td.composition().localToGlobalGasId(GAS, i);
+            label id = td.cloud().composition().localToGlobalGasId(GAS, i);
             td.cloud().rhoTrans(id)[cellI] += np0*dMassGas[i];
         }
         forAll(YLiquid_, i)
         {
-            label id = td.composition().localToGlobalGasId(LIQUID, i);
+            label id = td.cloud().composition().localToGlobalGasId(LIQUID, i);
             td.cloud().rhoTrans(id)[cellI] += np0*dMassLiquid[i];
         }
-        // No mapping between solid components and carrier phase
+//        // No mapping between solid components and carrier phase
 //        forAll(YSolid_, i)
 //        {
-//            label id = td.composition().localToGlobalGasId(SOLID, i);
+//            label id = td.cloud().composition().localToGlobalGasId(SOLID, i);
 //            td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
 //        }
         forAll(dMassSRCarrier, i)
@@ -268,8 +332,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
 
         // Update enthalpy transfer
-        // - enthalpy of lost solids already accounted for
-        td.cloud().hTrans()[cellI] += np0*(dhTrans + Sh);
+        td.cloud().hTrans()[cellI] += np0*dhTrans;
 
         // Coefficient to be applied in carrier phase enthalpy coupling
         td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
@@ -288,22 +351,24 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             // Absorb parcel into carrier phase
             forAll(YGas_, i)
             {
-                label id = td.composition().localToGlobalGasId(GAS, i);
+                label id = td.cloud().composition().localToGlobalGasId(GAS, i);
                 td.cloud().rhoTrans(id)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
             }
             forAll(YLiquid_, i)
             {
-                label id = td.composition().localToGlobalGasId(LIQUID, i);
+                label id =
+                    td.cloud().composition().localToGlobalGasId(LIQUID, i);
                 td.cloud().rhoTrans(id)[cellI] +=
                     np0
                    *mass1
                    *YMix[LIQUID]
                    *YLiquid_[i];
             }
-            // No mapping between solid components and carrier phase
+//            // No mapping between solid components and carrier phase
 //            forAll(YSolid_, i)
 //            {
-//                label id = td.composition().localToGlobalGasId(SOLID, i);
+//                label id =
+//                    td.cloud().composition().localToGlobalGasId(SOLID, i);
 //                td.cloud().rhoTrans(id)[cellI] +=
 //                    np0
 //                   *mass1
@@ -350,9 +415,12 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     const scalar dt,
     const scalar T,
     const scalar mass,
+    const scalar mass0,
     const label idVolatile,
-    const scalarField& YMixture,
-    scalarList& dMass
+    const scalar YVolatileTot,
+    const scalarField& YVolatile,
+    bool& canCombust,
+    scalarField& dMassDV
 )
 {
     // Check that model is active, and that the parcel temperature is
@@ -367,31 +435,27 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
         return 0.0;
     }
 
-    // Determine mass to add to carrier phase
+    // Total mass of volatiles evolved
     const scalar dMassTot = td.cloud().devolatilisation().calculate
     (
         dt,
-        this->mass0_,
+        mass0,
         mass,
         T,
         td.cloud().composition().YMixture0()[idVolatile],
-        YMixture[GAS],
-        canCombust_
+        YVolatileTot,
+        canCombust
     );
 
-    // Add to cummulative mass transfer
-    forAll (YGas_, i)
+    // Volatile mass transfer - equal components of each volatile specie
+    forAll(YVolatile, i)
     {
-        label id = td.cloud().composition().globalIds(idVolatile)[i];
-
-        // Volatiles mass transfer
-        scalar volatileMass = YGas_[i]*dMassTot;
-        dMass[i] += volatileMass;
+        dMassDV[i] = YVolatile[i]*dMassTot;
     }
 
     td.cloud().addToMassDevolatilisation(dMassTot);
 
-    return td.constProps().Ldevol()*dMassTot;
+    return -dMassTot*td.constProps().LDevol();
 }
 
 
@@ -402,36 +466,46 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
     TrackData& td,
     const scalar dt,
     const label cellI,
+    const scalar d,
     const scalar T,
-    const scalarList& dMassVolatile,
+    const scalar mass,
+    const bool canCombust,
+    const scalarField& dMassVolatile,
+    const scalarField& YMix,
+    const scalarField& YGas,
+    const scalarField& YLiquid,
+    const scalarField& YSolid,
     scalarField& dMassSRGas,
     scalarField& dMassSRLiquid,
     scalarField& dMassSRSolid,
-    scalarList& dMassSRCarrier,
+    scalarField& dMassSRCarrier,
     scalar& dhTrans
 )
 {
     // Check that model is active
-    if (!td.cloud().surfaceReaction().active() || !canCombust_)
+    if (!td.cloud().surfaceReaction().active() || !canCombust)
     {
-        return;
+        return 0.0;
     }
 
     // Update surface reactions
-    scalar HReaction = td.cloud().surfaceReaction().calculate
+    const scalar pc = this->pc_;
+    const scalar Tc = this->Tc_;
+    const scalar rhoc = this->rhoc_;
+    const scalar HReaction = td.cloud().surfaceReaction().calculate
     (
         dt,
         cellI,
-        this->d_,
+        d,
         T,
-        this->Tc_,
-        this->pc_,
-        this->rhoc_,
-        this->mass(),
-        YGas_,
-        YLiquid_,
-        YSolid_,
-        this->Y_,
+        Tc,
+        pc,
+        rhoc,
+        mass,
+        YGas,
+        YLiquid,
+        YSolid,
+        YMix,
         dMassVolatile,
         dMassSRGas,
         dMassSRLiquid,
@@ -439,23 +513,32 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
         dMassSRCarrier
     );
 
+    td.cloud().addToMassSurfaceReaction
+    (
+        sum(dMassSRGas) + sum(dMassSRLiquid) + sum(dMassSRSolid)
+    );
+
+    // Add enthalpy of consumed components to the carrier phase enthalpy
+    // transfer
+    const label idG = td.cloud().composition().idGas();
+    const label idL = td.cloud().composition().idLiquid();
+    const label idS = td.cloud().composition().idSolid();
+    scalar dhGas =
+        sum(dMassSRGas)*td.cloud().composition().H(idG, YGas, pc, T);
+
+    scalar dhLiquid =
+        sum(dMassSRLiquid)*td.cloud().composition().H(idL, YLiquid, pc, T);
+
+    scalar dhSolid =
+        sum(dMassSRSolid)*td.cloud().composition().H(idS, YSolid, pc, T);
+
+    dhTrans += dhGas + dhLiquid + dhSolid;
+
     // Heat of reaction divided between particle and carrier phase by the
     // fraction fh and (1 - fh)
-    dhTrans -= (1.0 - td.constProps().fh())*HReaction;
+    dhTrans += (1.0 - td.constProps().hRetentionCoeff())*HReaction;
 
-/*
-    // Correct dhTrans to account for enthalpy of consumed solids
-    dhTrans +=
-        sum(dMassSR)*td.cloud().composition().H(idSolid, YSolid_, pc, T0);
-
-    // Correct dhTrans to account for enthalpy of evolved volatiles
-    dhTrans +=
-        sum(dMassMT)*td.cloud().composition().H(idGas, YGas_, pc, T0);
-
-    // TODO: td.cloud().addToMassSurfaceReaction(sum(dMassSR));
-*/
-
-    return td.constProps().fh()*HReaction;
+    return td.constProps().hRetentionCoeff()*HReaction;
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index 58e2ef60dc..5053c106b4 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -83,11 +83,11 @@ public:
         // Private data
 
             //- Latent heat of devolatilisation [J/kg]
-            const scalar Ldevol_;
+            const scalar LDevol_;
 
             //- Fraction of enthalpy retained by parcel due to surface
             //  reactions
-            const scalar fh_;
+            const scalar hRetentionCoeff_;
 
 
     public:
@@ -98,11 +98,11 @@ public:
         // Access
 
             //- Return const access to the latent heat of devolatilisation
-            inline scalar Ldevol() const;
+            inline scalar LDevol() const;
 
             //- Return const access to the fraction of enthalpy retained by
-            // parcel due to surface reactions
-            inline scalar fh() const;
+            //  parcel due to surface reactions
+            inline scalar hRetentionCoeff() const;
     };
 
 
@@ -214,12 +214,15 @@ protected:
         scalar calcDevolatilisation
         (
             TrackData& td,
-            const scalar dt,
-            const scalar T,
-            const scalar mass,
-            const label idVolatile,
-            const scalarField& YMixture,
-            scalarList& dMass
+            const scalar dt,           // timestep
+            const scalar T,            // temperature
+            const scalar mass,         // mass
+            const scalar mass0,        // mass (initial on injection)
+            const label idVolatile,    // id of volatile phase
+            const scalar YVolatileTot, // total volatile mass fraction
+            const scalarField& YVolatile, // volatile component mass fractions
+            bool& canCombust,          // 'can combust' flag
+            scalarField& dMassDV       // mass transfer - local to particle
         );
 
         //- Calculate surface reactions
@@ -227,15 +230,22 @@ protected:
         scalar calcSurfaceReactions
         (
             TrackData& td,
-            const scalar dt,
-            const label cellI,
-            const scalar T,
-            const scalarList& dMassVolatile,
-            scalarField& dMassSRGas,
-            scalarField& dMassSRLiquid,
-            scalarField& dMassSRSolid,
-            scalarList& dMassSRc,
-            scalar& dhTrans
+            const scalar dt,           // timestep
+            const label cellI,         // owner cell
+            const scalar d,            // diameter
+            const scalar T,            // temperature
+            const scalar mass,         // mass
+            const bool canCombust,     // 'can combust' flag
+            const scalarField& dMassVolatile, // mass transfer of volatiles
+            const scalarField& YMix,    // mixture mass fractions
+            const scalarField& YGas,   // gas-phase mass fractions
+            const scalarField& YLiquid,// liquid-phase mass fractions
+            const scalarField& YSolid, // solid-phase mass fractions
+            scalarField& dMassSRGas,   // gas-phase mass transfer - local
+            scalarField& dMassSRLiquid,// liquid-phase mass transfer - local
+            scalarField& dMassSRSolid, // solid-phase mass transfer - local
+            scalarField& dMassSRCarrier,// carrier phase mass transfer
+            scalar& dhTrans            // enthalpy transfer to carrier phase
         );
 
 
@@ -253,16 +263,16 @@ public:
         inline ReactingMultiphaseParcel
         (
             ReactingMultiphaseCloud<ParcelType>& owner,
-            const label typeId,
             const vector& position,
             const label cellI,
+            const label typeId,
+            const scalar nParticle0,
             const scalar d0,
             const vector& U0,
-            const scalar nParticle0,
+            const scalarField& Y0,
             const scalarField& YGas0,
             const scalarField& YLiquid0,
             const scalarField& YSolid0,
-            const scalarField& Y0,
             const constantProperties& constProps
         );
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
index 1018f1e3f2..d7d3cf1781 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
@@ -34,8 +34,8 @@ constantProperties
 )
 :
     ReactingParcel<ParcelType>::constantProperties(dict),
-    Ldevol_(dimensionedScalar(dict.lookup("Ldevol")).value()),
-    fh_(dimensionedScalar(dict.lookup("fh")).value())
+    LDevol_(dimensionedScalar(dict.lookup("LDevol")).value()),
+    hRetentionCoeff_(dimensionedScalar(dict.lookup("hRetentionCoeff")).value())
 {}
 
 
@@ -74,28 +74,28 @@ template<class ParcelType>
 inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
 (
     ReactingMultiphaseCloud<ParcelType>& owner,
-    const label typeId,
     const vector& position,
     const label cellI,
+    const label typeId,
+    const scalar nParticle0,
     const scalar d0,
     const vector& U0,
-    const scalar nParticle0,
+    const scalarField& Y0,
     const scalarField& YGas0,
     const scalarField& YLiquid0,
     const scalarField& YSolid0,
-    const scalarField& Y0,
     const constantProperties& constProps
 )
 :
     ReactingParcel<ParcelType>
     (
         owner,
-        typeId,
         position,
         cellI,
+        typeId,
+        nParticle0,
         d0,
         U0,
-        nParticle0,
         Y0,
         constProps
     ),
@@ -109,17 +109,18 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
 
 template<class ParcelType>
 inline Foam::scalar
-Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::Ldevol() const
+Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::LDevol() const
 {
-    return Ldevol_;
+    return LDevol_;
 }
 
 
 template<class ParcelType>
 inline Foam::scalar
-Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::fh() const
+Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
+hRetentionCoeff() const
 {
-    return fh_;
+    return hRetentionCoeff_;
 }
 
 
@@ -152,13 +153,6 @@ YGas() const
 }
 
 
-template<class ParcelType>
-inline Foam::scalarField& Foam::ReactingMultiphaseParcel<ParcelType>::YGas()
-{
-    return YGas_;
-}
-
-
 template<class ParcelType>
 inline const Foam::scalarField& Foam::ReactingMultiphaseParcel<ParcelType>::
 YLiquid() const
@@ -167,13 +161,6 @@ YLiquid() const
 }
 
 
-template<class ParcelType>
-inline Foam::scalarField& Foam::ReactingMultiphaseParcel<ParcelType>::YLiquid()
-{
-    return YLiquid_;
-}
-
-
 template<class ParcelType>
 inline const Foam::scalarField& Foam::ReactingMultiphaseParcel<ParcelType>::
 YSolid() const
@@ -182,6 +169,20 @@ YSolid() const
 }
 
 
+template<class ParcelType>
+inline Foam::scalarField& Foam::ReactingMultiphaseParcel<ParcelType>::YGas()
+{
+    return YGas_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalarField& Foam::ReactingMultiphaseParcel<ParcelType>::YLiquid()
+{
+    return YLiquid_;
+}
+
+
 template<class ParcelType>
 inline Foam::scalarField& Foam::ReactingMultiphaseParcel<ParcelType>::YSolid()
 {
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
index cbb952181c..64829732e8 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
@@ -45,7 +45,7 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
     if (readFields)
     {
         const ReactingMultiphaseCloud<ParcelType>& cR =
-            dynamic_cast<const ReactingMultiphaseCloud<ParcelType>& >(cloud);
+            dynamic_cast<const ReactingMultiphaseCloud<ParcelType>&>(cloud);
 
         const label idGas = cR.composition().idGas();
         const label nGas = cR.composition().componentNames(idGas).size();
@@ -68,7 +68,7 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
         }
 
         // scale the mass fractions
-        const scalarList& YMix = this->Y_;
+        const scalarField& YMix = this->Y_;
         YGas_ /= YMix[GAS] + ROOTVSMALL;
         YLiquid_ /= YMix[LIQUID] + ROOTVSMALL;
         YSolid_ /= YMix[SOLID] + ROOTVSMALL;
@@ -282,14 +282,14 @@ Foam::Ostream& Foam::operator<<
     scalarField YSolidLoc = p.YSolid()*p.Y()[2];
     if (os.format() == IOstream::ASCII)
     {
-        os  << static_cast<const ReactingMultiphaseParcel<ParcelType>& >(p)
+        os  << static_cast<const ReactingMultiphaseParcel<ParcelType>&>(p)
             << token::SPACE << YGasLoc
             << token::SPACE << YLiquidLoc
             << token::SPACE << YSolidLoc;
     }
     else
     {
-        os  << static_cast<const ReactingMultiphaseParcel<ParcelType>& >(p);
+        os  << static_cast<const ReactingMultiphaseParcel<ParcelType>&>(p);
         os << YGasLoc << YLiquidLoc << YSolidLoc;
     }
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 676e0f10b2..66bc6242b1 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -59,7 +59,7 @@ template<class ParcelType>
 void Foam::ReactingParcel<ParcelType>::updateMassFraction
 (
     const scalar mass0,
-    const scalarList& dMass,
+    const scalarField& dMass,
     scalarField& Y
 )
 {
@@ -83,9 +83,13 @@ void Foam::ReactingParcel<ParcelType>::calc
 {
     // Define local properties at beginning of time step
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const scalar mass0 = this->mass();
     const scalar np0 = this->nParticle_;
+    const scalar d0 = this->d_;
+    const vector& U0 = this->U_;
+    const scalar rho0 = this->rho_;
     const scalar T0 = this->T_;
+    const scalar cp0 = this->cp_;
+    const scalar mass0 = this->mass();
 
 
     // Intialise transfer terms
@@ -98,14 +102,15 @@ void Foam::ReactingParcel<ParcelType>::calc
     scalar dhTrans = 0.0;
 
     // Mass transfer due to phase change
-    scalarList dMassPC(td.cloud().gases().size(), 0.0);
+    scalarField dMassPC(td.cloud().gases().size(), 0.0);
 
 
     // Phase change
     // ~~~~~~~~~~~~
 
     // Return enthalpy source and calc mass transfer due to phase change
-    scalar ShPC = calcPhaseChange(td, dt, cellI, T0, 0, 1.0, Y_, dMassPC);
+    scalar ShPC =
+        calcPhaseChange(td, dt, cellI, d0, T0, U0, 0, 1.0, Y_, dMassPC);
 
     // Update particle component mass fractions
     updateMassFraction(mass0, dMassPC, Y_);
@@ -116,7 +121,21 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     // Calculate new particle temperature
     scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, dt, cellI, ShPC, htc, dhTrans);
+    scalar T1 =
+        calcHeatTransfer
+        (
+            td,
+            dt,
+            cellI,
+            d0,
+            U0,
+            rho0,
+            T0,
+            cp0,
+            ShPC,
+            htc,
+            dhTrans
+        );
 
 
     // Motion
@@ -131,7 +150,19 @@ void Foam::ReactingParcel<ParcelType>::calc
     // Calculate new particle velocity
     scalar Cud = 0.0;
     vector U1 =
-        calcVelocity(td, dt, cellI, Fx, 0.5*(mass0 + mass1), Cud, dUTrans);
+        calcVelocity
+        (
+            td,
+            dt,
+            cellI,
+            d0,
+            U0,
+            rho0,
+            0.5*(mass0 + mass1),
+            Fx,
+            Cud,
+            dUTrans
+        );
 
 
     // Accumulate carrier phase source terms
@@ -151,7 +182,7 @@ void Foam::ReactingParcel<ParcelType>::calc
         td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
 
         // Update enthalpy transfer
-        td.cloud().hTrans()[cellI] += np0*(dhTrans + ShPC);
+        td.cloud().hTrans()[cellI] += np0*dhTrans;
 
         // Coefficient to be applied in carrier phase enthalpy coupling
         td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
@@ -208,11 +239,13 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
     TrackData& td,
     const scalar dt,
     const label cellI,
+    const scalar d,
     const scalar T,
+    const vector& U,
     const label idPhase,
     const scalar YPhase,
     const scalarField& YComponents,
-    scalarList& dMass
+    scalarField& dMassPC
 )
 {
     if
@@ -229,16 +262,16 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
     (
         dt,
         cellI,
+        d,
         min(T, td.constProps().Tbp()), // Limiting to boiling temperature
         pc_,
-        this->d_,
         this->Tc_,
         this->muc_/this->rhoc_,
-        this->U_ - this->Uc_,
-        dMass
+        U - this->Uc_,
+        dMassPC
     );
 
-    scalar dMassTot = sum(dMass);
+    scalar dMassTot = sum(dMassPC);
 
     // Add to cumulative phase change mass
     td.cloud().addToMassPhaseChange(dMassTot);
@@ -246,10 +279,7 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
     // Effective latent heat of vaporisation
     scalar LEff = td.cloud().composition().L(idPhase, YComponents, pc_, T);
 
-    // Effective heat of vaporised components
-    scalar HEff = td.cloud().composition().H(idPhase, YComponents, pc_, T);
-
-    return dMassTot*(HEff - LEff);
+    return -dMassTot*LEff;
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index bc38527efb..52b9c881e7 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -169,7 +169,7 @@ protected:
             scalar mass0_;
 
             //- Mass fractions of mixture []
-            scalarList Y_;
+            scalarField Y_;
 
 
         // Cell-based quantities
@@ -185,20 +185,22 @@ protected:
         scalar calcPhaseChange
         (
             TrackData& td,
-            const scalar dt,
-            const label cellI,
-            const scalar T,
-            const label idPhase,
-            const scalar YPhase,
-            const scalarField& YComponents,
-            scalarList& dMass
+            const scalar dt,           // timestep
+            const label cellI,         // owner cell
+            const scalar d,            // diameter
+            const scalar T,            // temperature
+            const vector& U,           // velocity
+            const label idPhase,       // id of phase involved in phase change
+            const scalar YPhase,       // total mass fraction
+            const scalarField& YComponents, // component mass fractions
+            scalarField& dMassPC       // mass transfer - local to particle
         );
 
         //- Update mass fraction
         void updateMassFraction
         (
             const scalar mass0,
-            const scalarList& dMass,
+            const scalarField& dMass,
             scalarField& Y
         );
 
@@ -217,12 +219,12 @@ public:
         inline ReactingParcel
         (
             ReactingCloud<ParcelType>& owner,
-            const label typeId,
             const vector& position,
             const label cellI,
+            const label typeId,
+            const scalar nParticle0,
             const scalar d0,
             const vector& U0,
-            const scalar nParticle0,
             const scalarField& Y0,
             const constantProperties& constProps
         );
@@ -250,7 +252,7 @@ public:
             inline scalar mass0() const;
 
             //- Return const access to mass fractions of mixture
-            inline const scalarList& Y() const;
+            inline const scalarField& Y() const;
 
             //- Return the owner cell pressure
             inline scalar pc() const;
@@ -262,7 +264,7 @@ public:
             inline scalar& mass0();
 
             //- Return access to mass fractions of mixture
-            inline scalarList& Y();
+            inline scalarField& Y();
 
 
         // Main calculation loop
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
index aeb38c9c5c..da808caee8 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
@@ -74,12 +74,12 @@ template<class ParcelType>
 inline Foam::ReactingParcel<ParcelType>::ReactingParcel
 (
     ReactingCloud<ParcelType>& owner,
-    const label typeId,
     const vector& position,
     const label cellI,
+    const label typeId,
+    const scalar nParticle0,
     const scalar d0,
     const vector& U0,
-    const scalar nParticle0,
     const scalarField& Y0,
     const constantProperties& constProps
 )
@@ -87,12 +87,12 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
     ThermoParcel<ParcelType>
     (
         owner,
-        typeId,
         position,
         cellI,
+        typeId,
+        nParticle0,
         d0,
         U0,
-        nParticle0,
         constProps
     ),
     mass0_(0.0),
@@ -166,7 +166,7 @@ inline Foam::scalar Foam::ReactingParcel<ParcelType>::mass0() const
 
 
 template<class ParcelType>
-inline const Foam::scalarList& Foam::ReactingParcel<ParcelType>::Y() const
+inline const Foam::scalarField& Foam::ReactingParcel<ParcelType>::Y() const
 {
     return Y_;
 }
@@ -187,7 +187,7 @@ inline Foam::scalar& Foam::ReactingParcel<ParcelType>::mass0()
 
 
 template<class ParcelType>
-inline Foam::scalarList& Foam::ReactingParcel<ParcelType>::Y()
+inline Foam::scalarField& Foam::ReactingParcel<ParcelType>::Y()
 {
     return Y_;
 }
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
index eacf49b112..968c608d7e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
@@ -45,7 +45,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
     if (readFields)
     {
         const ReactingCloud<ParcelType>& cR =
-            dynamic_cast<const ReactingCloud<ParcelType>& >(cloud);
+            dynamic_cast<const ReactingCloud<ParcelType>&>(cloud);
 
         const label nMixture = cR.composition().phaseTypes().size();
         Y_.setSize(nMixture);
@@ -206,13 +206,13 @@ Foam::Ostream& Foam::operator<<
 {
     if (os.format() == IOstream::ASCII)
     {
-        os  << static_cast<const ThermoParcel<ParcelType>& >(p)
+        os  << static_cast<const ThermoParcel<ParcelType>&>(p)
             << token::SPACE << p.mass0()
             << token::SPACE << p.Y();
     }
     else
     {
-        os  << static_cast<const ThermoParcel<ParcelType>& >(p);
+        os  << static_cast<const ThermoParcel<ParcelType>&>(p);
         os.write
         (
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 276fcd0b59..e13420e56d 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -68,9 +68,13 @@ void Foam::ThermoParcel<ParcelType>::calc
 {
     // Define local properties at beginning of time step
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const vector U0 = this->U_;
-    const scalar mass0 = this->mass();
     const scalar np0 = this->nParticle_;
+    const scalar d0 = this->d_;
+    const vector U0 = this->U_;
+    const scalar rho0 = this->rho_;
+    const scalar T0 = this->T_;
+    const scalar cp0 = this->cp_;
+    const scalar mass0 = this->mass();
 
 
     // Initialise transfer terms
@@ -84,11 +88,25 @@ void Foam::ThermoParcel<ParcelType>::calc
     // ~~~~~~~~~~~~~
 
     // No additional enthalpy sources
-    vector Sh = 0.0;
+    scalar Sh = 0.0;
 
     // Calculate new particle velocity
     scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, dt, cellI, Sh, htc, dhTrans);
+    scalar T1 =
+        calcHeatTransfer
+        (
+            td,
+            dt,
+            cellI,
+            d0,
+            U0,
+            rho0,
+            T0,
+            cp0,
+            Sh,
+            htc,
+            dhTrans
+        );
 
 
     // Motion
@@ -99,7 +117,8 @@ void Foam::ThermoParcel<ParcelType>::calc
 
     // Calculate new particle velocity
     scalar Cud = 0.0;
-    vector U1 = calcVelocity(td, dt, cellI, Fx, Cud, mass0, dUTrans);
+    vector U1 =
+        calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx, Cud, dUTrans);
 
 
     //  Accumulate carrier phase source terms
@@ -133,6 +152,11 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     TrackData& td,
     const scalar dt,
     const label cellI,
+    const scalar d,
+    const vector& U,
+    const scalar rho,
+    const scalar T,
+    const scalar cp,
     const scalar Sh,
     scalar& htc,
     scalar& dhTrans
@@ -141,31 +165,31 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     if (!td.cloud().heatTransfer().active())
     {
         htc = 0.0;
-        return T_;
+        return T;
     }
 
     // Calc heat transfer coefficient
     htc = td.cloud().heatTransfer().h
     (
-        this->d_,
-        this->U_ - this->Uc_,
+        d,
+        U - this->Uc_,
         this->rhoc_,
-        this->rho_,
+        rho,
         cpc_,
-        cp_,
+        cp,
         this->muc_
     );
 
 
-    // Set new particle temperature
-    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Determine new particle temperature
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Assuming diameter = diameter at start of time step
-    scalar Ap = this->areasS();
+    // Particle area
+    scalar Ap = this->areaS(d);
 
     // Determine ap and bp coefficients
     scalar ap = Tc_ + Sh/(htc*Ap + ROOTVSMALL);
-    scalar bp = 6.0*htc/(this->rho_*this->d_*cp_);
+    scalar bp = 6.0*htc/(rho*d*cp);
     if (td.cloud().radiation())
     {
         // Carrier phase incident radiation field
@@ -177,7 +201,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
         // Helper variables
         const scalar sigma = radiation::sigmaSB.value();
         const scalar epsilon = td.constProps().epsilon0();
-        const scalar D = epsilon*sigma*pow3(T_)/(htc + ROOTVSMALL) + 1.0;
+        const scalar D = epsilon*sigma*pow3(T)/(htc + ROOTVSMALL) + 1.0;
         ap += 0.25*epsilon*G[cellI]/(htc + ROOTVSMALL);
         ap /= D;
         bp *= D;
@@ -185,7 +209,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     // Integrate to find the new parcel temperature
     IntegrationScheme<scalar>::integrationResult Tres =
-        td.cloud().TIntegrator().integrate(T_, dt, ap, bp);
+        td.cloud().TIntegrator().integrate(T, dt, ap, bp);
 
     // Enthalpy transfer
     // - Using average particle temperature
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index 8506bb7c0a..4bc5180eff 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -208,11 +208,16 @@ protected:
         scalar calcHeatTransfer
         (
             TrackData& td,
-            const scalar dt,
-            const label cellI,
-            const scalar Sh,
-            scalar& htc,
-            scalar& dhTrans
+            const scalar dt,           // timestep
+            const label cellI,         // owner cell
+            const scalar d,            // diameter
+            const vector& U,           // velocity
+            const scalar rho,          // density
+            const scalar T,            // temperature
+            const scalar cp,           // specific heat capacity
+            const scalar Sh,           // additional enthalpy sources
+            scalar& htc,               // heat transfer coeff
+            scalar& dhTrans            // enthalpy transfer to carrier phase
         );
 
 
@@ -230,12 +235,12 @@ public:
         inline ThermoParcel
         (
             ThermoCloud<ParcelType>& owner,
-            const label typeId,
             const vector& position,
             const label cellI,
+            const label typeId,
+            const scalar nParticle0,
             const scalar d0,
             const vector& U0,
-            const scalar nParticle0,
             const constantProperties& constProps
         );
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
index c6f0f2b816..c14f1891b4 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
@@ -73,24 +73,24 @@ template<class ParcelType>
 inline Foam::ThermoParcel<ParcelType>::ThermoParcel
 (
     ThermoCloud<ParcelType>& owner,
-    const label typeId,
     const vector& position,
     const label cellI,
+    const label typeId,
+    const scalar nParticle0,
     const scalar d0,
     const vector& U0,
-    const scalar nParticle0,
     const constantProperties& constProps
 )
 :
     KinematicParcel<ParcelType>
     (
         owner,
-        typeId,
         position,
         cellI,
+        typeId,
+        nParticle0,
         d0,
         U0,
-        nParticle0,
         constProps
     ),
     T_(constProps.T0()),
@@ -178,16 +178,16 @@ inline Foam::scalar Foam::ThermoParcel<ParcelType>::T() const
 
 
 template<class ParcelType>
-inline Foam::scalar& Foam::ThermoParcel<ParcelType>::T()
+inline Foam::scalar Foam::ThermoParcel<ParcelType>::cp() const
 {
-    return T_;
+    return cp_;
 }
 
 
 template<class ParcelType>
-inline Foam::scalar Foam::ThermoParcel<ParcelType>::cp() const
+inline Foam::scalar& Foam::ThermoParcel<ParcelType>::T()
 {
-    return cp_;
+    return T_;
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
index 08433679bb..d4f30ce917 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
@@ -140,13 +140,13 @@ Foam::Ostream& Foam::operator<<
 {
     if (os.format() == IOstream::ASCII)
     {
-        os  << static_cast<const KinematicParcel<ParcelType>& >(p)
+        os  << static_cast<const KinematicParcel<ParcelType>&>(p)
             << token::SPACE << p.T()
             << token::SPACE << p.cp();
     }
     else
     {
-        os  << static_cast<const KinematicParcel<ParcelType>& >(p);
+        os  << static_cast<const KinematicParcel<ParcelType>&>(p);
         os.write
         (
             reinterpret_cast<const char*>(&p.T_),
diff --git a/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.C b/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.C
index c0d438df24..d207887826 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.C
@@ -41,24 +41,24 @@ namespace Foam
 Foam::basicKinematicParcel::basicKinematicParcel
 (
     KinematicCloud<basicKinematicParcel>& owner,
-    const label typeId,
     const vector& position,
     const label cellI,
+    const label typeId,
+    const scalar nParticle0,
     const scalar d0,
     const vector& U0,
-    const scalar nParticle0,
     const constantProperties& constProps
 )
 :
     KinematicParcel<basicKinematicParcel>
     (
         owner,
-        typeId,
         position,
         cellI,
+        typeId,
+        nParticle0,
         d0,
         U0,
-        nParticle0,
         constProps
     )
 {}
diff --git a/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.H b/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.H
index b5da2267e7..0e5733870f 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.H
@@ -44,7 +44,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class basicKinematicParcel Declaration
+                   Class basicKinematicParcel Declaration
 \*---------------------------------------------------------------------------*/
 
 class basicKinematicParcel
@@ -64,12 +64,12 @@ public:
         basicKinematicParcel
         (
             KinematicCloud<basicKinematicParcel>& owner,
-            const label typeId,
             const vector& position,
             const label cellI,
+            const label typeId,
+            const scalar nParticle0,
             const scalar d0,
             const vector& U0,
-            const scalar nParticle0,
             const constantProperties& constProps
         );
 
@@ -90,8 +90,7 @@ public:
 
 
     //- Destructor
-
-        virtual ~basicKinematicParcel();
+    virtual ~basicKinematicParcel();
 };
 
 
diff --git a/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.C
index 0f7ce9173f..d200738319 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.C
@@ -41,12 +41,12 @@ namespace Foam
 Foam::basicReactingMultiphaseParcel::basicReactingMultiphaseParcel
 (
     ReactingMultiphaseCloud<basicReactingMultiphaseParcel>& owner,
-    const label typeId,
     const vector& position,
     const label cellI,
+    const label typeId,
+    const scalar nParticle0,
     const scalar d0,
     const vector& U0,
-    const scalar nParticle0,
     const scalarField& YGas0,
     const scalarField& YLiquid0,
     const scalarField& YSolid0,
@@ -57,12 +57,12 @@ Foam::basicReactingMultiphaseParcel::basicReactingMultiphaseParcel
     ReactingMultiphaseParcel<basicReactingMultiphaseParcel>
     (
         owner,
-        typeId,
         position,
         cellI,
+        typeId,
+        nParticle0,
         d0,
         U0,
-        nParticle0,
         YGas0,
         YLiquid0,
         YSolid0,
diff --git a/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.H
index bb5053ff74..cc537e5638 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.H
@@ -64,12 +64,12 @@ public:
         basicReactingMultiphaseParcel
         (
              ReactingMultiphaseCloud<basicReactingMultiphaseParcel>& owner,
-             const label typeId,
              const vector& position,
              const label cellI,
+             const label typeId,
+             const scalar nParticle0,
              const scalar d0,
              const vector& U0,
-             const scalar nParticle0,
              const scalarField& YGas0,
              const scalarField& YLiquid0,
              const scalarField& YSolid0,
@@ -93,7 +93,7 @@ public:
         }
 
 
-    // Destructor
+    //- Destructor
     virtual ~basicReactingMultiphaseParcel();
 };
 
diff --git a/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.C b/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.C
index 722c572e43..72ac1f4842 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.C
@@ -41,12 +41,12 @@ namespace Foam
 Foam::basicReactingParcel::basicReactingParcel
 (
     ReactingCloud<basicReactingParcel>& owner,
-    const label typeId,
     const vector& position,
     const label cellI,
+    const label typeId,
+    const scalar nParticle0,
     const scalar d0,
     const vector& U0,
-    const scalar nParticle0,
     const scalarField& Y0,
     const constantProperties& constProps
 )
@@ -54,12 +54,12 @@ Foam::basicReactingParcel::basicReactingParcel
     ReactingParcel<basicReactingParcel>
     (
         owner,
-        typeId,
         position,
         cellI,
+        typeId,
+        nParticle0,
         d0,
         U0,
-        nParticle0,
         Y0,
         constProps
     )
diff --git a/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.H b/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.H
index 008a33b25b..fae88a2da7 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.H
@@ -64,12 +64,12 @@ public:
         basicReactingParcel
         (
             ReactingCloud<basicReactingParcel>& owner,
-            const label typeId,
             const vector& position,
             const label cellI,
+            const label typeId,
+            const scalar nParticle0,
             const scalar d0,
             const vector& U0,
-            const scalar nParticle0,
             const scalarField& Y0,
             const constantProperties& constProps
         );
diff --git a/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.C b/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.C
index 76a575d1de..095b217aca 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.C
@@ -41,24 +41,24 @@ namespace Foam
 Foam::basicThermoParcel::basicThermoParcel
 (
     ThermoCloud<basicThermoParcel>& owner,
-    const label typeId,
     const vector position,
     const label cellI,
+    const label typeId,
+    const scalar nParticle0,
     const scalar d0,
     const vector U0,
-    const scalar nParticle0,
     const constantProperties& constProps
 )
 :
     ThermoParcel<basicThermoParcel>
     (
         owner,
-        typeId,
         position,
         cellI,
+        typeId,
+        nParticle0,
         d0,
         U0,
-        nParticle0,
         constProps
     )
 {}
diff --git a/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.H b/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.H
index c55471d3d1..ec74d2add8 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.H
@@ -44,7 +44,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class basicThermoParcel Declaration
+                      Class basicThermoParcel Declaration
 \*---------------------------------------------------------------------------*/
 
 class basicThermoParcel
@@ -63,12 +63,12 @@ public:
        basicThermoParcel
        (
             ThermoCloud<basicThermoParcel>& owner,
-            const label typeId,
             const vector position,
             const label cellI,
+            const label typeId,
+            const scalar nParticle0,
             const scalar d0,
             const vector U0,
-            const scalar nParticle0,
             const constantProperties& constProps
         );
 
@@ -87,9 +87,8 @@ public:
         }
 
 
-    // Destructors
-
-        virtual ~basicThermoParcel();
+    //- Destructor
+    virtual ~basicThermoParcel();
 };
 
 
diff --git a/src/lagrangian/intermediate/submodels/IO/DataEntry/polynomial/polynomial.H b/src/lagrangian/intermediate/submodels/IO/DataEntry/polynomial/polynomial.H
index f405b9c54e..6d3cd51f65 100644
--- a/src/lagrangian/intermediate/submodels/IO/DataEntry/polynomial/polynomial.H
+++ b/src/lagrangian/intermediate/submodels/IO/DataEntry/polynomial/polynomial.H
@@ -26,7 +26,7 @@ Class
     Foam::polynomial
 
 Description
-    Templated polynomial container data entry. Items are stored in a list of
+    Polynomial container data entry for scalars. Items are stored in a list of
     Tuple2's. Data is input in the form, e.g. for an entry <entryName> that
     describes y = x^2 + 2x^3
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
index 4bc06eaf84..d71e98d2b3 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
@@ -241,7 +241,7 @@ Foam::label Foam::CompositionModel<CloudType>::localId
 
 
 template<class CloudType>
-Foam::label Foam::CompositionModel<CloudType>::localToGlobalGaslId
+Foam::label Foam::CompositionModel<CloudType>::localToGlobalGasId
 (
     const label phaseI,
     const label id
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
index fc01732314..f6a41acf0b 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
@@ -188,7 +188,7 @@ public:
                 label localId(const label phaseI, const word& cmptName) const;
 
                 //- Return global gas id of component given local id
-                label localToGlobalGaslId
+                label localToGlobalGasId
                 (
                     const label phaseI,
                     const label id
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
index 830c4d7f1d..d1c3e8b8a8 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
@@ -135,13 +135,13 @@ void Foam::LiquidEvaporation<CloudType>::calculate
 (
     const scalar dt,
     const label cellI,
+    const scalar d,
     const scalar T,
     const scalar pc,
-    const scalar d,
     const scalar Tc,
     const scalar nuc,
     const vector& Ur,
-    scalarList& dMass
+    scalarField& dMassPC
 ) const
 {
     // construct carrier phase species volume fractions for cell, cellI
@@ -184,7 +184,7 @@ void Foam::LiquidEvaporation<CloudType>::calculate
         scalar Ni = max(kc*(Cs - Cinf), 0.0);
 
         // mass transfer [kg]
-        dMass[lid] += Ni*A*liquids_->properties()[lid].W()*dt;
+        dMassPC[lid] += Ni*A*liquids_->properties()[lid].W()*dt;
     }
 }
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
index c70c081e90..408fb3d06b 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
@@ -41,7 +41,7 @@ Description
 namespace Foam
 {
 /*---------------------------------------------------------------------------*\
-                       Class LiquidEvaporation Declaration
+                     Class LiquidEvaporation Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class CloudType>
@@ -105,13 +105,13 @@ public:
         (
             const scalar dt,
             const label cellI,
+            const scalar d,
             const scalar T,
             const scalar pc,
-            const scalar d,
             const scalar Tc,
             const scalar nuc,
             const vector& Ur,
-            scalarList& dMass
+            scalarField& dMassPC
         ) const;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C
index 7ecb03c773..ee3779b22e 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C
@@ -66,7 +66,7 @@ void Foam::NoPhaseChange<CloudType>::calculate
     const scalar,
     const scalar,
     const vector&,
-    scalarList&
+    scalarField&
 ) const
 {
     // Nothing to do...
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H
index e1ec25dd5e..6419f127ff 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H
@@ -74,13 +74,13 @@ public:
         (
             const scalar dt,
             const label cellI,
+            const scalar d,
             const scalar T,
             const scalar pc,
-            const scalar d,
             const scalar Tc,
             const scalar nuc,
             const vector& Ur,
-            scalarList& dMass
+            scalarField& dMassPC
         ) const;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
index 32e44a960f..b8ebc644b0 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
@@ -140,13 +140,13 @@ public:
         (
             const scalar dt,
             const label cellI,
+            const scalar d,
             const scalar T,
             const scalar pc,
-            const scalar d,
             const scalar Tc,
             const scalar nuc,
             const vector& Ur,
-            scalarList& dMass
+            scalarField& dMassPC
         ) const = 0;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
index 0630158e28..53eecf01f9 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
@@ -69,12 +69,12 @@ Foam::scalar Foam::NoSurfaceReaction<CloudType>::calculate
     const scalarField&,
     const scalarField&,
     const scalarField&,
-    const scalarList&,
-    const scalarList&,
+    const scalarField&,
+    const scalarField&,
     scalarField&,
     scalarField&,
     scalarField&,
-    scalarList&
+    scalarField&
 ) const
 {
     // do nothing
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
index 6566f60456..d3c99ebd5d 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
@@ -87,12 +87,12 @@ public:
             const scalarField& YGas,
             const scalarField& YLiquid,
             const scalarField& YSolid,
-            const scalarList& YMixture,
-            const scalarList& dMassVolatile,
+            const scalarField& YMixture,
+            const scalarField& dMassVolatile,
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,
-            scalarList& dMassSRCarrier
+            scalarField& dMassSRCarrier
         ) const;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
index 269c37de37..2deacb3465 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
@@ -42,7 +42,6 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 
 #include "scalarField.H"
-#include "scalarList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -146,12 +145,12 @@ public:
             const scalarField& YGas,
             const scalarField& YLiquid,
             const scalarField& YSolid,
-            const scalarList& YMixture,
-            const scalarList& dMassVolatile,
+            const scalarField& YMixture,
+            const scalarField& dMassVolatile,
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,
-            scalarList& dMassSRCarrier
+            scalarField& dMassSRCarrier
         ) const = 0;
 };