From 7677c87f1518ed3ac6c9edb12db4541c537656a6 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Thu, 5 Mar 2009 15:48:51 +0000
Subject: [PATCH] further updates en route to new evap model(s)

---
 src/lagrangian/intermediate/Make/files        |  14 +-
 .../evaporationProperties.C                   |  80 -------
 .../evaporationProperties.H                   | 117 ----------
 .../evaporationPropertiesIO.C                 | 130 -----------
 .../KinematicParcel/KinematicParcel.C         |  96 +++-----
 .../KinematicParcel/KinematicParcel.H         |  23 +-
 .../KinematicParcel/KinematicParcelI.H        |   4 +-
 .../ReactingMultiphaseParcel.C                | 170 +++-----------
 .../ReactingMultiphaseParcel.H                |  31 ++-
 .../ReactingMultiphaseParcelI.H               |  19 +-
 .../ReactingMultiphaseParcelIO.C              |  20 +-
 .../Templates/ReactingParcel/ReactingParcel.C | 219 ++++++------------
 .../Templates/ReactingParcel/ReactingParcel.H |  40 ++--
 .../ReactingParcel/ReactingParcelI.H          |  16 +-
 .../ReactingParcel/ReactingParcelIO.C         |  26 +--
 .../Templates/ThermoParcel/ThermoParcel.C     |  81 ++-----
 .../Templates/ThermoParcel/ThermoParcel.H     |  21 +-
 .../Templates/ThermoParcel/ThermoParcelI.H    |   4 +-
 .../basicKinematicParcel.C                    |   4 +-
 .../basicKinematicParcel.H                    |   2 +-
 .../basicReactingMultiphaseParcel.C           |   8 +-
 .../basicReactingMultiphaseParcel.H           |   4 +-
 .../basicReactingParcel/basicReactingParcel.C |   8 +-
 .../basicReactingParcel/basicReactingParcel.H |   4 +-
 .../basicThermoParcel/basicThermoParcel.C     |   4 +-
 .../basicThermoParcel/basicThermoParcel.H     |   2 +-
 .../IO/DataEntry/polynomial/polynomial.C      |   7 +
 .../LiquidEvaporation/LiquidEvaporation.C     |  66 ++++--
 .../LiquidEvaporation/LiquidEvaporation.H     |  20 +-
 .../NoPhaseChange/NoPhaseChange.C             |   8 +-
 .../NoPhaseChange/NoPhaseChange.H             |  10 +-
 .../PhaseChangeModel/PhaseChangeModel.H       |  10 +-
 32 files changed, 353 insertions(+), 915 deletions(-)
 delete mode 100644 src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationProperties.C
 delete mode 100644 src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationProperties.H
 delete mode 100644 src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationPropertiesIO.C

diff --git a/src/lagrangian/intermediate/Make/files b/src/lagrangian/intermediate/Make/files
index bc6b095488a..93572202fc6 100644
--- a/src/lagrangian/intermediate/Make/files
+++ b/src/lagrangian/intermediate/Make/files
@@ -7,18 +7,21 @@ $(DERIVEDPARCELS)/basicThermoParcel/basicThermoParcel.C
 $(DERIVEDPARCELS)/basicReactingParcel/basicReactingParcel.C
 $(DERIVEDPARCELS)/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.C
 
+
 /* Cloud base classes */
 clouds/baseClasses/kinematicCloud/kinematicCloud.C
 clouds/baseClasses/thermoCloud/thermoCloud.C
 clouds/baseClasses/reactingCloud/reactingCloud.C
 clouds/baseClasses/reactingMultiphaseCloud/reactingMultiphaseCloud.C
 
+
 /* Cloud container/injection mechanisms */
 clouds/derived/basicKinematicCloud/basicKinematicCloud.C
 clouds/derived/basicThermoCloud/basicThermoCloud.C
 clouds/derived/basicReactingCloud/basicReactingCloud.C
 clouds/derived/basicReactingMultiphaseCloud/basicReactingMultiphaseCloud.C
 
+
 /* kinematic parcel sub-models */
 KINEMATICPARCEL=$(DERIVEDPARCELS)/basicKinematicParcel
 $(KINEMATICPARCEL)/defineBasicKinematicParcel.C
@@ -27,6 +30,7 @@ $(KINEMATICPARCEL)/submodels/makeBasicKinematicParcelDragModels.C
 $(KINEMATICPARCEL)/submodels/makeBasicKinematicParcelInjectionModels.C
 $(KINEMATICPARCEL)/submodels/makeBasicKinematicParcelWallInteractionModels.C
 
+
 /* thermo parcel sub-models */
 THERMOPARCEL=$(DERIVEDPARCELS)/basicThermoParcel
 $(THERMOPARCEL)/defineBasicThermoParcel.C
@@ -36,6 +40,7 @@ $(THERMOPARCEL)/submodels/makeBasicThermoParcelHeatTransferModels.C
 $(THERMOPARCEL)/submodels/makeBasicThermoParcelInjectionModels.C
 $(THERMOPARCEL)/submodels/makeBasicThermoParcelWallInteractionModels.C
 
+
 /* reacting parcel sub-models */
 REACTINGPARCEL=$(DERIVEDPARCELS)/basicReactingParcel
 $(REACTINGPARCEL)/defineBasicReactingParcel.C
@@ -47,6 +52,7 @@ $(REACTINGPARCEL)/submodels/makeBasicReactingParcelInjectionModels.C
 $(REACTINGPARCEL)/submodels/makeBasicReactingParcelPhaseChangeModels.C
 $(REACTINGPARCEL)/submodels/makeBasicReactingParcelWallInteractionModels.C
 
+
 /* reacting multiphase parcel sub-models */
 REACTINGMPPARCEL=$(DERIVEDPARCELS)/basicReactingMultiphaseParcel
 $(REACTINGMPPARCEL)/defineBasicReactingMultiphaseParcel.C
@@ -60,25 +66,25 @@ $(REACTINGMPPARCEL)/submodels/makeBasicReactingMultiphaseParcelPhaseChangeModels
 $(REACTINGMPPARCEL)/submodels/makeBasicReactingMultiphaseParcelSurfaceReactionModels.C
 $(REACTINGMPPARCEL)/submodels/makeBasicReactingMultiphaseParcelWallInteractionModels.C
 
+
 /* bolt-on models */
 submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
 submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
 
+
 /* integration schemes */
 IntegrationScheme/makeIntegrationSchemes.C
 
+
 /* phase properties */
 phaseProperties/phaseProperties/phaseProperties.C
 phaseProperties/phaseProperties/phasePropertiesIO.C
 phaseProperties/phasePropertiesList/phasePropertiesList.C
 
-/* evaporation properties */
-evaporationProperties/evaporationProperties/evaporationProperties.C
-evaporationProperties/evaporationProperties/evaporationPropertiesIO.C
-
 
 /* data entries */
 submodels/IO/DataEntry/makeDataEntries.C
 submodels/IO/DataEntry/polynomial/polynomial.C
 
+
 LIB = $(FOAM_LIBBIN)/liblagrangianIntermediate
diff --git a/src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationProperties.C b/src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationProperties.C
deleted file mode 100644
index 49a353a18a2..00000000000
--- a/src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationProperties.C
+++ /dev/null
@@ -1,80 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "evaporationProperties.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::evaporationProperties::evaporationProperties()
-:
-    name_("unknownSpecie"),
-    Dab_(0.0),
-    TvsPSat_()
-{}
-
-
-Foam::evaporationProperties::evaporationProperties
-(
-    const evaporationProperties& pp
-)
-:
-    name_(pp.name_),
-    Dab_(pp.Dab_),
-    TvsPSat_(pp.TvsPSat_)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::evaporationProperties::~evaporationProperties()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-const Foam::word& Foam::evaporationProperties::name() const
-{
-    return name_;
-}
-
-
-Foam::scalar Foam::evaporationProperties::Dab() const
-{
-    return Dab_;
-}
-
-
-const Foam::DataEntry<Foam::scalar>&
-Foam::evaporationProperties::TvsPSat() const
-{
-    return TvsPSat_();
-}
-
-
-// ************************************************************************* //
-
diff --git a/src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationProperties.H b/src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationProperties.H
deleted file mode 100644
index 1ee2ed1c6ab..00000000000
--- a/src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationProperties.H
+++ /dev/null
@@ -1,117 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::evaporationProperties
-
-Description
-    Helper class to manage evaporation properties
-
-SourceFiles
-    evaporationProperties.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef evaporationProperties_H
-#define evaporationProperties_H
-
-#include "DataEntry.H"
-#include "autoPtr.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                    Class evaporationProperties Declaration
-\*---------------------------------------------------------------------------*/
-
-class evaporationProperties
-{
-   // Private data
-
-        //- Name of specie
-        word name_;
-
-        //- Diffusion coefficient
-        scalar Dab_;
-
-        //- Data entry for saturation pressure as a function of temperature
-        autoPtr<DataEntry<scalar> > TvsPSat_;
-
-
-    // Private member functions
-
-
-
-public:
-
-    // Constructors
-
-        //- Null constructor
-        evaporationProperties();
-
-        //- Construct from Istream
-        evaporationProperties(Istream&);
-
-        //- Construct as copy
-        evaporationProperties(const evaporationProperties&);
-
-
-    //- Destructor
-    ~evaporationProperties();
-
-
-    // Public member functions
-
-        // Access
-
-            //- Return const access to the specie name
-            const word& name() const;
-
-            //- Return const access to the diffusion coefficient
-            scalar Dab() const;
-
-            //- Return const access to the saturation pressure as a function
-            //  of temperature
-            const DataEntry<scalar>& TvsPSat() const;
-
-
-    // IOstream Operators
-
-        friend Istream& operator>>(Istream&, evaporationProperties&);
-        friend Ostream& operator<<(Ostream&, const evaporationProperties&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationPropertiesIO.C b/src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationPropertiesIO.C
deleted file mode 100644
index 3eed20d37bb..00000000000
--- a/src/lagrangian/intermediate/evaporationProperties/evaporationProperties/evaporationPropertiesIO.C
+++ /dev/null
@@ -1,130 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "evaporationProperties.H"
-#include "dictionaryEntry.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::evaporationProperties::evaporationProperties(Istream& is)
-:
-    name_("unknown"),
-    Dab_(0.0),
-    TvsPSat_(NULL)
-{
-    is.check
-    (
-        "Foam::evaporationProperties::evaporationProperties(Istream& is)"
-    );
-
-    dictionaryEntry evapInfo(dictionary::null, is);
-
-    if (!evapInfo.isDict())
-    {
-        FatalErrorIn
-        (
-            "Foam::evaporationProperties::evaporationProperties(Istream& is)"
-        )   << "Evaporation properties should be given in dictionary entries"
-            << nl << exit(FatalError);
-    }
-
-    name_ = evapInfo.keyword();
-
-    evapInfo.lookup("Dab") >> Dab_;
-    TvsPSat_.reset(DataEntry<scalar>::New("TvsPSat", evapInfo.dict()).ptr());
-}
-
-
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
-
-Foam::Istream& Foam::operator>>(Istream& is, evaporationProperties& ep)
-{
-    is.check
-    (
-        "Foam::Istream& Foam::operator>>"
-        "("
-            "Foam::Istream&,"
-            "Foam::evaporationProperties&"
-        ")"
-    );
-
-    dictionaryEntry evapInfo(dictionary::null, is);
-
-    if (!evapInfo.isDict())
-    {
-        FatalErrorIn
-        (
-            "Foam::Istream& Foam::operator>>"
-            "("
-                "Istream& is,"
-                "evaporationProperties& pp"
-            ")"
-        )   << "Evaporation properties should be given in dictionary entries"
-            << nl << exit(FatalError);
-    }
-
-    ep.name_ = evapInfo.keyword();
-
-    evapInfo.lookup("Dab") >> ep.Dab_;
-    ep.TvsPSat_.reset(DataEntry<scalar>::New("TvsPSat", evapInfo.dict()).ptr());
-
-    return is;
-}
-
-
-Foam::Ostream& Foam::operator<<(Ostream& os, const evaporationProperties& ep)
-{
-    os.check
-    (
-        "Foam::Ostream& Foam::operator<<"
-        "("
-            "Foam::Ostream&,"
-            "const Foam::evaporationProperties&"
-        ")"
-    );
-
-    os  << ep.name_ << nl << token::BEGIN_BLOCK << nl
-        << incrIndent;
-
-    os.writeKeyword("Dab") << ep.Dab_ << token::END_STATEMENT << nl;
-//    os  <<  ep.TvsPSat_() << nl;
-
-    os  << decrIndent << token::END_BLOCK << nl;
-
-    os.check
-    (
-        "Foam::Ostream& Foam::operator<<"
-        "("
-            "Foam::Ostream&,"
-            "const Foam::evaporationProperties&"
-        ")"
-    );
-
-    return os;
-}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index bec040ca837..32d40fc5923 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -35,18 +35,18 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
 (
     TrackData& td,
     const scalar dt,
-    const label celli
+    const label cellI
 )
 {
-    rhoc_ = td.rhoInterp().interpolate(this->position(), celli);
-    Uc_ = td.UInterp().interpolate(this->position(), celli);
-    muc_ = td.muInterp().interpolate(this->position(), celli);
+    rhoc_ = td.rhoInterp().interpolate(this->position(), cellI);
+    Uc_ = td.UInterp().interpolate(this->position(), cellI);
+    muc_ = td.muInterp().interpolate(this->position(), cellI);
 
     // Apply dispersion components to carrier phase velocity
     Uc_ = td.cloud().dispersion().update
     (
         dt,
-        celli,
+        cellI,
         U_,
         Uc_,
         UTurb_,
@@ -57,44 +57,33 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
 
 template<class ParcelType>
 template<class TrackData>
-void Foam::KinematicParcel<ParcelType>::calcCoupled
+void Foam::KinematicParcel<ParcelType>::calc
 (
     TrackData& td,
     const scalar dt,
-    const label celli
+    const label cellI
 )
 {
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Define local properties at beginning of timestep
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//    const scalar mass0 = mass();
-//    const vector U0 = U_;
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Initialise transfer terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Momentum transfer from the particle to the carrier phase
-    vector dUTrans = vector::zero;
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate velocity - update U
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Calculate velocity - update U, dUTrans
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
+    vector dUTrans = vector::zero;
     const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~
-    // Accumulate source terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~
+    if (td.cloud().coupled())
+    {
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // Accumulate carrier phase source terms
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Update momentum transfer
-    td.cloud().UTrans()[celli] += nParticle_*dUTrans;
+        // Update momentum transfer
+        td.cloud().UTrans()[cellI] += nParticle_*dUTrans;
 
-    // Accumulate coefficient to be applied in carrier phase momentum coupling
-    td.cloud().UCoeff()[celli] += nParticle_*mass()*Cud;
+        // Coefficient to be applied in carrier phase momentum coupling
+        td.cloud().UCoeff()[cellI] += nParticle_*mass()*Cud;
+    }
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -106,42 +95,18 @@ void Foam::KinematicParcel<ParcelType>::calcCoupled
 
 template<class ParcelType>
 template<class TrackData>
-void Foam::KinematicParcel<ParcelType>::calcUncoupled
-(
-    TrackData& td,
-    const scalar dt,
-    const label
-)
-{
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Initialise transfer terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Momentum transfer from the particle to the carrier phase
-    vector dUTrans = vector::zero;
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate velocity - update U
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Cud = 0.0;
-    this->U() = calcVelocity(td, dt, Cud, dUTrans);
-}
-
-
-template<class ParcelType>
-template<class TrackData>
-Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
+const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 (
     TrackData& td,
     const scalar dt,
     scalar& Cud,
     vector& dUTrans
-)
+) const
 {
     // Return linearised term from drag model
     Cud = td.cloud().drag().Cu(U_ - Uc_, d_, rhoc_, rho_, muc_);
 
+
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Set new particle velocity
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -188,7 +153,7 @@ bool Foam::KinematicParcel<ParcelType>::move
 
         // Remember which cell the Parcel is in
         // since this will change if a face is hit
-        label celli = p.cell();
+        label cellI = p.cell();
 
         dt *= p.trackToFace(p.position() + dt*U_, td);
 
@@ -196,19 +161,12 @@ bool Foam::KinematicParcel<ParcelType>::move
         p.stepFraction() = 1.0 - tEnd/deltaT;
 
         // Update cell based properties
-        p.updateCellQuantities(td, dt, celli);
+        p.updateCellQuantities(td, dt, cellI);
 
         // Avoid problems with extremely small timesteps
         if (dt > ROOTVSMALL)
         {
-            if (td.cloud().coupled())
-            {
-                p.calcCoupled(td, dt, celli);
-            }
-            else
-            {
-                p.calcUncoupled(td, dt, celli);
-            }
+            p.calc(td, dt, cellI);
         }
 
         if (p.onBoundary() && td.keepParticle)
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index 39ba450cacd..c60a2164984 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -224,13 +224,13 @@ protected:
 
         //- Calculate new particle velocity
         template<class TrackData>
-        vector calcVelocity
+        const vector calcVelocity
         (
             TrackData& td,
             const scalar dt,
             scalar& Cud,
             vector& dUTrans
-        );
+        ) const;
 
 
 public:
@@ -249,7 +249,7 @@ public:
             KinematicCloud<ParcelType>& owner,
             const label typeId,
             const vector& position,
-            const label celli,
+            const label cellI,
             const scalar d0,
             const vector& U0,
             const scalar nParticle0,
@@ -355,25 +355,16 @@ public:
             (
                 TrackData& td,
                 const scalar dt,
-                const label celli
-            );
-
-            //- Coupled calculation with the continuous phase
-            template<class TrackData>
-            void calcCoupled
-            (
-                TrackData& td,
-                const scalar dt,
-                const label celli
+                const label cellI
             );
 
-            //- Uncoupled calculation with the continuous phase
+            //- Update parcel properties over the time interval
             template<class TrackData>
-            void calcUncoupled
+            void calc
             (
                 TrackData& td,
                 const scalar dt,
-                const label
+                const label cellI
             );
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index 4b4ec3cce77..3e8d49b7ad4 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -64,14 +64,14 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
     KinematicCloud<ParcelType>& owner,
     const label typeId,
     const vector& position,
-    const label celli,
+    const label cellI,
     const scalar d0,
     const vector& U0,
     const scalar nParticle0,
     const constantProperties& constProps
 )
 :
-    Particle<ParcelType>(owner, position, celli),
+    Particle<ParcelType>(owner, position, cellI),
     typeId_(typeId),
     d_(d0),
     U_(U0),
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index bc132bb6f34..4d96e786018 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -34,20 +34,20 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::updateCellQuantities
 (
     TrackData& td,
     const scalar dt,
-    const label celli
+    const label cellI
 )
 {
-    ReactingParcel<ParcelType>::updateCellQuantities(td, dt, celli);
+    ReactingParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
 }
 
 
 template<class ParcelType>
 template<class TrackData>
-void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
+void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 (
     TrackData& td,
     const scalar dt,
-    const label celli
+    const label cellI
 )
 {
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -58,7 +58,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
     const scalar cp0 = this->cp_;
     const scalar np0 = this->nParticle_;
     const scalar T0 = this->T_;
-    scalarField& YMix = this->YMixture_;
+    scalarField& YMix = this->Y_;
 
     label idGas = td.cloud().composition().idGas();
     label idLiquid = td.cloud().composition().idLiquid();
@@ -94,14 +94,13 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
     // Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
+    scalar T1 = calcHeatTransfer(td, dt, cellI, htc, dhTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~
     // Calculate phase change
     // ~~~~~~~~~~~~~~~~~~~~~~
-    scalarField X = td.cloud().composition().X(idLiquid, YLiquid_);
-    calcPhaseChange(td, dt, T0, X, dMassMT);
+    scalar dMassPC = calcPhaseChange(td, dt, cellI, T1, dMassMT);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -116,7 +115,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
     // Initialise enthalpy retention to zero
     scalar dhRet = 0.0;
 
-    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
+    calcSurfaceReactions(td, dt, cellI, T0, T1, dMassMT, dMassSR, dhRet);
 
     // New total mass
     const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
@@ -128,10 +127,10 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
        *sum(dMassMT)
        /(0.5*(mass0 + mass1)*cp0);
 
-    // Add retained enthalpy from surface reaction to particle and remove
-    // from gas
-    T1 += dhRet/(0.5*(mass0 + mass1)*cp0);
-    dhTrans -= dhRet;
+    // Enthalpy retention divided between particle and carrier phase by the
+    // fraction fh and (1 - fh)
+    T1 += td.constProps().fh()*dhRet/(0.5*(mass0 + mass1)*cp0);
+    dhTrans -= (1.0 - td.constProps().fh())*dhRet;
 
     // Correct dhTrans to account for enthalpy of evolved volatiles
     dhTrans +=
@@ -151,21 +150,21 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
     // Transfer mass lost from particle to carrier mass source
     forAll(dMassMT, i)
     {
-        td.cloud().rhoTrans(i)[celli] += np0*(dMassMT[i] + dMassSR[i]);
+        td.cloud().rhoTrans(i)[cellI] += np0*(dMassMT[i] + dMassSR[i]);
     }
 
     // Update momentum transfer
-    td.cloud().UTrans()[celli] += np0*dUTrans;
+    td.cloud().UTrans()[cellI] += np0*dUTrans;
 
     // Accumulate coefficient to be applied in carrier phase momentum coupling
-    td.cloud().UCoeff()[celli] += np0*mass0*Cud;
+    td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
 
     // Update enthalpy transfer
     // - enthalpy of lost solids already accounted for
-    td.cloud().hTrans()[celli] += np0*dhTrans;
+    td.cloud().hTrans()[cellI] += np0*dhTrans;
 
     // Accumulate coefficient to be applied in carrier phase enthalpy coupling
-    td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
+    td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -178,9 +177,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
         // Absorb particle(s) into carrier phase
         forAll(dMassMT, i)
         {
-            td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i];
+            td.cloud().rhoTrans(i)[cellI] += np0*dMassMT[i];
         }
-        td.cloud().hTrans()[celli] +=
+        td.cloud().hTrans()[cellI] +=
             np0*mass1
            *(
                 YMix[0]
@@ -188,7 +187,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
               + YMix[2]
                *td.cloud().composition().H(idSolid, YSolid_, this->pc_, T1)
             );
-        td.cloud().UTrans()[celli] += np0*mass1*U1;
+        td.cloud().UTrans()[cellI] += np0*mass1*U1;
     }
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Set new particle properties
@@ -215,123 +214,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
 }
 
 
-template<class ParcelType>
-template<class TrackData>
-void Foam::ReactingMultiphaseParcel<ParcelType>::calcUncoupled
-(
-    TrackData& td,
-    const scalar dt,
-    const label celli
-)
-{
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Define local properties at beginning of timestep
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const scalar T0 = this->T_;
-    const scalar mass0 = this->mass();
-    const scalar cp0 = this->cp_;
-    scalarField& YMix = this->YMixture_;
-
-    label idGas = td.cloud().composition().idGas();
-    label idLiquid = td.cloud().composition().idLiquid();
-    label idSolid = td.cloud().composition().idSolid();
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Initialise transfer terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Momentum transfer from the particle to the carrier phase
-    vector dUTrans = vector::zero;
-
-    // Enthalpy transfer from the particle to the carrier phase
-    scalar dhTrans = 0.0;
-
-    // Mass transfer from particle to carrier phase
-    // - components exist in particle already
-    scalarList dMassMT(td.cloud().gases().size(), 0.0);
-
-    // Mass transfer due to surface reactions from particle to carrier phase
-    // - components do not necessarily exist in particle already
-    scalarList dMassSR(td.cloud().gases().size(), 0.0);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate velocity - update U
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate heat transfer - update T
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate phase change
-    // ~~~~~~~~~~~~~~~~~~~~~~
-    scalarField X = td.cloud().composition().X(idLiquid, YLiquid_);
-    calcPhaseChange(td, dt, T0, YLiquid_, dMassMT);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate Devolatilisation
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-    calcDevolatilisation(td, dt, T0, T1, dMassMT);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate surface reactions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Initialise enthalpy retention to zero
-    scalar dhRet = 0.0;
-
-    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
-
-    // New total mass
-    const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
-
-    // New specific heat capacity
-    const scalar cp1 =
-        YMix[0]*td.cloud().composition().cp(idGas, YGas_, this->pc_, T1)
-      + YMix[1]*td.cloud().composition().cp(idLiquid, YLiquid_, this->pc_, T1)
-      + YMix[2]*td.cloud().composition().cp(idSolid, YSolid_, this->pc_, T1);
-
-    // Add retained enthalpy to particle
-    T1 += dhRet/(mass0*0.5*(cp0 + cp1));
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Remove the particle when mass falls below minimum threshold
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (mass1 < td.constProps().minParticleMass())
-    {
-        td.keepParticle = false;
-    }
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Set new particle properties
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    else
-    {
-        this->U_ = U1;
-        this->T_ = T1;
-        this->cp_ = cp1;
-
-        // Update particle density or diameter
-        if (td.constProps().constantVolume())
-        {
-            this->rho_ = mass1/this->volume();
-        }
-        else
-        {
-            this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
-        }
-    }
-}
-
-
 template<class ParcelType>
 template<class TrackData>
 void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
@@ -358,7 +240,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
             ")\n"
             "no treatment currently available for particles containing "
             "liquid species"
-        )
+        );
     }
 
     // Check that model is active, and that the parcel temperature is
@@ -375,7 +257,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 
     // Determine mass to add to carrier phase
     const scalar mass = this->mass();
-    scalarField& YMix = this->YMixture_;
+    scalarField& YMix = this->Y_;
     const scalar dMassTot = td.cloud().devolatilisation().calculate
     (
         dt,
@@ -413,7 +295,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
 (
     TrackData& td,
     const scalar dt,
-    const label celli,
+    const label cellI,
     const scalar T0,
     const scalar T1,
     const scalarList& dMassMT,
@@ -432,7 +314,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
     td.cloud().surfaceReaction().calculate
     (
         dt,
-        celli,
+        cellI,
         this->d_,
         T0,
         T1,
@@ -444,12 +326,12 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
         YGas_,
         YLiquid_,
         YSolid_,
-        this->YMixture_,
+        this->Y_,
         dMassSR,
         dhRet
     );
 
-    // TODO: td.cloud().addToMassDevolatilisation(sum(dMassSR));
+    // TODO: td.cloud().addToMassSurfaceReaction(sum(dMassSR));
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index a0d9ea528f3..d53fc2a6d02 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -79,6 +79,10 @@ public:
             //- Latent heat of devolatilisation [J/kg]
             const scalar Ldevol_;
 
+            //- Fraction of enthalpy retained by parcel due to surface
+            //  reactions
+            const scalar fh_;
+
 
     public:
 
@@ -89,6 +93,10 @@ public:
 
             //- Return const access to the latent heat of devolatilisation
             inline scalar Ldevol() const;
+
+            //- Return const access to the fraction of enthalpy retained by
+            // parcel due to surface reactions
+            inline scalar fh() const;
     };
 
 
@@ -174,7 +182,7 @@ protected:
         (
             TrackData& td,
             const scalar dt,
-            const label celli,
+            const label cellI,
             const scalar T0,
             const scalar T1,
             const scalarList& dMassMT,
@@ -199,14 +207,14 @@ public:
             ReactingMultiphaseCloud<ParcelType>& owner,
             const label typeId,
             const vector& position,
-            const label celli,
+            const label cellI,
             const scalar d0,
             const vector& U0,
             const scalar nParticle0,
             const scalarField& YGas0,
             const scalarField& YLiquid0,
             const scalarField& YSolid0,
-            const scalarField& YMixture0,
+            const scalarField& Y0,
             const constantProperties& constProps
         );
 
@@ -262,25 +270,16 @@ public:
             (
                 TrackData& td,
                 const scalar dt,
-                const label celli
-            );
-
-            //- Coupled calculation with the continuous phase
-            template<class TrackData>
-            void calcCoupled
-            (
-                TrackData& td,
-                const scalar dt,
-                const label celli
+                const label cellI
             );
 
-            //- Uncoupled calculation with the continuous phase
+            //- Update parcel properties over the time interval
             template<class TrackData>
-            void calcUncoupled
+            void calc
             (
                 TrackData& td,
                 const scalar dt,
-                const label
+                const label cellI
             );
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
index a87ad6032b9..1018f1e3f2b 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
@@ -34,7 +34,8 @@ constantProperties
 )
 :
     ReactingParcel<ParcelType>::constantProperties(dict),
-    Ldevol_(dimensionedScalar(dict.lookup("Ldevol")).value())
+    Ldevol_(dimensionedScalar(dict.lookup("Ldevol")).value()),
+    fh_(dimensionedScalar(dict.lookup("fh")).value())
 {}
 
 
@@ -75,14 +76,14 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
     ReactingMultiphaseCloud<ParcelType>& owner,
     const label typeId,
     const vector& position,
-    const label celli,
+    const label cellI,
     const scalar d0,
     const vector& U0,
     const scalar nParticle0,
     const scalarField& YGas0,
     const scalarField& YLiquid0,
     const scalarField& YSolid0,
-    const scalarField& YMixture0,
+    const scalarField& Y0,
     const constantProperties& constProps
 )
 :
@@ -91,11 +92,11 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
         owner,
         typeId,
         position,
-        celli,
+        cellI,
         d0,
         U0,
         nParticle0,
-        YMixture0,
+        Y0,
         constProps
     ),
     YGas_(YGas0),
@@ -114,6 +115,14 @@ Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::Ldevol() const
 }
 
 
+template<class ParcelType>
+inline Foam::scalar
+Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::fh() const
+{
+    return fh_;
+}
+
+
 // * * * * * * * * * * * trackData Member Functions  * * * * * * * * * * * * //
 
 template<class ParcelType>
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
index bb53f10a205..e94d3b9177a 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
@@ -68,7 +68,7 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
         }
 
         // scale the mass fractions
-        const scalarField& YMix = this->YMixture_;
+        const scalarField& YMix = this->Y_;
         YGas_ /= YMix[0] + ROOTVSMALL;
         YLiquid_ /= YMix[1] + ROOTVSMALL;
         YSolid_ /= YMix[2] + ROOTVSMALL;
@@ -129,7 +129,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         forAllIter(typename Cloud<ParcelType>, c, iter)
         {
             ReactingMultiphaseParcel<ParcelType>& p = iter();
-            p.YGas_[j] = YGas[i++]/(p.YMixture()[0] + ROOTVSMALL);
+            p.YGas_[j] = YGas[i++]/(p.Y()[0] + ROOTVSMALL);
         }
     }
     // Populate YLiquid for each parcel
@@ -144,7 +144,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         forAllIter(typename Cloud<ParcelType>, c, iter)
         {
             ReactingMultiphaseParcel<ParcelType>& p = iter();
-            p.YLiquid_[j] = YLiquid[i++]/(p.YMixture()[1] + ROOTVSMALL);
+            p.YLiquid_[j] = YLiquid[i++]/(p.Y()[1] + ROOTVSMALL);
         }
     }
     // Populate YSolid for each parcel
@@ -159,7 +159,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         forAllIter(typename Cloud<ParcelType>, c, iter)
         {
             ReactingMultiphaseParcel<ParcelType>& p = iter();
-            p.YSolid_[j] = YSolid[i++]/(p.YMixture()[2] + ROOTVSMALL);
+            p.YSolid_[j] = YSolid[i++]/(p.Y()[2] + ROOTVSMALL);
         }
     }
 }
@@ -192,7 +192,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
             forAllConstIter(typename Cloud<ParcelType>, c, iter)
             {
                 const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
-                YGas[i++] = p0.YGas()[j]*p0.YMixture()[0];
+                YGas[i++] = p0.YGas()[j]*p0.Y()[0];
             }
 
             YGas.write();
@@ -211,7 +211,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
             forAllConstIter(typename Cloud<ParcelType>, c, iter)
             {
                 const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
-                YLiquid[i++] = p0.YLiquid()[j]*p0.YMixture()[1];
+                YLiquid[i++] = p0.YLiquid()[j]*p0.Y()[1];
             }
 
             YLiquid.write();
@@ -230,7 +230,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
             forAllConstIter(typename Cloud<ParcelType>, c, iter)
             {
                 const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
-                YSolid[i++] = p0.YSolid()[j]*p0.YMixture()[2];
+                YSolid[i++] = p0.YSolid()[j]*p0.Y()[2];
             }
 
             YSolid.write();
@@ -248,9 +248,9 @@ Foam::Ostream& Foam::operator<<
     const ReactingMultiphaseParcel<ParcelType>& p
 )
 {
-    scalarField YGasLoc = p.YGas()*p.YMixture()[0];
-    scalarField YLiquidLoc = p.YLiquid()*p.YMixture()[1];
-    scalarField YSolidLoc = p.YSolid()*p.YMixture()[2];
+    scalarField YGasLoc = p.YGas()*p.Y()[0];
+    scalarField YLiquidLoc = p.YLiquid()*p.Y()[1];
+    scalarField YSolidLoc = p.YSolid()*p.Y()[2];
     if (os.format() == IOstream::ASCII)
     {
         os  << static_cast<const ReactingMultiphaseParcel<ParcelType>& >(p)
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 01fa099a925..6240d503bae 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -34,22 +34,38 @@ void Foam::ReactingParcel<ParcelType>::updateCellQuantities
 (
     TrackData& td,
     const scalar dt,
-    const label celli
+    const label cellI
 )
 {
-    ThermoParcel<ParcelType>::updateCellQuantities(td, dt, celli);
+    ThermoParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
 
-    pc_ = td.pInterp().interpolate(this->position(), celli);
+    pc_ = td.pInterp().interpolate(this->position(), cellI);
+}
+
+
+template<class ParcelType>
+void Foam::ReactingParcel<ParcelType>::updateMassFraction
+(
+    const scalar mass0,
+    const scalar mass1,
+    const scalarList& dMass,
+    scalarField& Y
+)
+{
+    forAll(Y, i)
+    {
+        Y[i] = (Y[i]*mass0 - dMass[i])/mass1;
+    }
 }
 
 
 template<class ParcelType>
 template<class TrackData>
-void Foam::ReactingParcel<ParcelType>::calcCoupled
+void Foam::ReactingParcel<ParcelType>::calc
 (
     TrackData& td,
     const scalar dt,
-    const label celli
+    const label cellI
 )
 {
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -61,7 +77,6 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     const scalar np0 = this->nParticle_;
     const scalar T0 = this->T_;
 
-
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
     // Initialise transfer terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -73,7 +88,6 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     scalar dhTrans = 0.0;
 
     // Mass transfer from particle to carrier phase
-    // - components exist in particle already
     scalarList dMassMT(td.cloud().gases().size(), 0.0);
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -87,147 +101,58 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
+    scalar T1 = calcHeatTransfer(td, dt, cellI, htc, dhTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate phase change
-    // ~~~~~~~~~~~~~~~~~~~~~~
-    scalarField X = td.cliud().composition().X(0, YMixture_);
-    calcPhaseChange(td, dt, T, X, dMassMT);
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Calculate phase change - update mass, Y, cp, T, dhTrans
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    const scalar dMassPC = calcPhaseChange(td, dt, cellI, T1, dMassMT);
 
-    // New total mass
-    const scalar mass1 = mass0 - sum(dMassMT);
+    // Update particle mass
+    scalar mass1 = mass0 - dMassPC;
 
+    // Update particle component mass fractions
+    updateMassFraction(mass0, mass1, dMassMT, Y_);
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~
-    // Accumulate source terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Transfer mass lost from particle to carrier mass source
-    forAll(dMassMT, i)
-    {
-        td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i];
-    }
+    // New specific heat capacity of
+    scalar cp1 = td.cloud().composition().cp(0, Y_, pc_, T1);
 
-    // Update momentum transfer
-    td.cloud().UTrans()[celli] += np0*dUTrans;
+    // Correct temperature due to evaporated components
+    // TODO: use hl function in liquidMixture???
+    T1 -= td.constProps().Lvap()*dMassPC/(0.5*(mass0 + mass1)*cp1);
 
-    // Accumulate coefficient to be applied in carrier phase momentum coupling
-    td.cloud().UCoeff()[celli] += np0*mass0*Cud;
+    // Correct dhTrans to account for the change in enthalpy due to the
+    // phase change
+    scalar HMix = td.cloud().composition().H(0, Y_, pc, 0.5*(T0 + T1));
+    dhTrans += dMassPC*HMix;
 
-    // Update enthalpy transfer
-    td.cloud().hTrans()[celli] += np0*dhTrans;
 
-    // Accumulate coefficient to be applied in carrier phase enthalpy coupling
-    td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Remove the particle when mass falls below minimum threshold
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (mass1 < td.constProps().minParticleMass())
+    if (td.cloud().coupled())
     {
-        td.keepParticle = false;
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // Accumulate carrier phase source terms
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        // Absorb particle(s) into carrier phase
+        // Transfer mass lost from particle to carrier mass source
         forAll(dMassMT, i)
         {
-            td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i];
+            td.cloud().rhoTrans(i)[cellI] += np0*dMassMT[i];
         }
-        td.cloud().UTrans()[celli] += np0*mass1*U1;
-    }
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Set new particle properties
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    else
-    {
-        this->U_ = U1;
-        this->T_ = T1;
-        //        this->cp_ = ??? // TODO:
 
-        // Update particle density or diameter
-        if (td.constProps().constantVolume())
-        {
-            this->rho_ = mass1/this->volume();
-        }
-        else
-        {
-            this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
-        }
-    }
-}
+        // Update momentum transfer
+        td.cloud().UTrans()[cellI] += np0*dUTrans;
 
+        // Coefficient to be applied in carrier phase momentum coupling
+        td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
 
-template<class ParcelType>
-template<class TrackData>
-void Foam::ReactingParcel<ParcelType>::calcUncoupled
-(
-    TrackData& td,
-    const scalar dt,
-    const label celli
-)
-{
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Define local properties at beginning of timestep
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const scalar T0 = this->T_;
-    const scalar mass0 = this->mass();
-    const scalar cp0 = this->cp_;
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Initialise transfer terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Momentum transfer from the particle to the carrier phase
-    vector dUTrans = vector::zero;
-
-    // Enthalpy transfer from the particle to the carrier phase
-    scalar dhTrans = 0.0;
-
-    // Mass transfer from particle to carrier phase
-    // - components exist in particle already
-    scalarList dMassMT(td.cloud().gases().size(), 0.0);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate velocity - update U
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
+        // Update enthalpy transfer
+        td.cloud().hTrans()[cellI] += np0*dhTrans;
 
+        // Coefficient to be applied in carrier phase enthalpy coupling
+        td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
+    }
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate heat transfer - update T
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate phase change
-    // ~~~~~~~~~~~~~~~~~~~~~~
-    scalarField X = td.cloud().composition().X(0, YMixture_);
-    scalar dMassPC = calcPhaseChange(td, dt, T, X, dMassMT);
-    T1 -= td.constProps().Lvap()*dMassPC/(0.5*mass0*cp0);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate surface reactions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Initialise enthalpy retention to zero
-    scalar dhRet = 0.0;
-
-    // New total mass
-    const scalar mass1 = mass0 - sum(dMassMT);
-
-    // New specific heat capacity
-    const scalar cp1 = cp0; // TODO: new cp1
-
-    // Add retained enthalpy to particle
-    T1 += dhRet/(mass0*0.5*(cp0 + cp1));
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Remove the particle when mass falls below minimum threshold
@@ -235,6 +160,17 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     if (mass1 < td.constProps().minParticleMass())
     {
         td.keepParticle = false;
+
+        if (td.cloud().coupled())
+        {
+            // Absorb particle(s) into carrier phase
+            forAll(dMassMT, i)
+            {
+                td.cloud().rhoTrans(i)[cellI] += np0*mass1*Y_[i];
+            }
+            td.cloud().UTrans()[cellI] += np0*mass1*U1;
+            td.cloud().hTrans()[cellI] += np0*mass1*HMix;
+        }
     }
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Set new particle properties
@@ -260,39 +196,36 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
 
 template<class ParcelType>
 template<class TrackData>
-void Foam::ReactingParcel<ParcelType>::calcPhaseChange
+Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
 (
     TrackData& td,
     const scalar dt,
+    const label cellI,
     const scalar T,
-    scalarField& X,
     scalarList& dMassMT
 )
 {
     if (!td.cloud().phaseChange().active())
     {
-        return;
+        return 0.0;
     }
 
-    // TODO: separate treatment for boiling
-
-    scalar dMassTot = td.cloud().phaseChange().calculate
+    scalar dMass = td.cloud().phaseChange().calculate
     (
+        dt,
+        cellI,
         T,
+        pc_,
         this->d_,
-        X,
-        dMassMT,
-        this->U_ - this->Uc_,
         this->Tc_,
-        pc_,
         this->muc_/this->rhoc_,
-        dt
+        this->U_ - this->Uc_,
+        dMassMT
     );
 
-    td.cloud().addToMassPhaseChange(dMassTot);
-    // TODO: Re-calculate mass fractions
-
+    td.cloud().addToMassPhaseChange(dMass);
 
+    return dMass;
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index a4e19fbbae1..6848a58d24e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -175,7 +175,7 @@ protected:
             scalar mass0_;
 
             //- Mass fractions of mixture []
-            scalarField YMixture_;
+            scalarField Y_;
 
 
         // Cell-based quantities
@@ -188,15 +188,24 @@ protected:
 
         //- Calculate Phase change
         template<class TrackData>
-        void calcPhaseChange
+        scalar calcPhaseChange
         (
             TrackData& td,
             const scalar dt,
+            const label cellI,
             const scalar T,
-            scalarField& X,
             scalarList& dMassMT
         );
 
+        //- Update mass fraction
+        void updateMassFraction
+        (
+            const scalar mass0,
+            const scalar mass1,
+            const scalarList& dMass,
+            scalarField& Y
+        );
+
 
 public:
 
@@ -214,11 +223,11 @@ public:
             ReactingCloud<ParcelType>& owner,
             const label typeId,
             const vector& position,
-            const label celli,
+            const label cellI,
             const scalar d0,
             const vector& U0,
             const scalar nParticle0,
-            const scalarField& YMixture0,
+            const scalarField& Y0,
             const constantProperties& constProps
         );
 
@@ -245,7 +254,7 @@ public:
             inline scalar mass0() const;
 
             //- Return const access to mass fractions of mixture
-            inline const scalarField& YMixture() const;
+            inline const scalarField& Y() const;
 
             //- Return the owner cell pressure
             inline scalar pc() const;
@@ -257,7 +266,7 @@ public:
             inline scalar& mass0();
 
             //- Return access to mass fractions of mixture
-            inline scalarField& YMixture();
+            inline scalarField& Y();
 
 
         // Main calculation loop
@@ -268,25 +277,16 @@ public:
             (
                 TrackData& td,
                 const scalar dt,
-                const label celli
-            );
-
-            //- Coupled calculation with the continuous phase
-            template<class TrackData>
-            void calcCoupled
-            (
-                TrackData& td,
-                const scalar dt,
-                const label celli
+                const label cellI
             );
 
-            //- Uncoupled calculation with the continuous phase
+            //- Update parcel properties over the time interval
             template<class TrackData>
-            void calcUncoupled
+            void calc
             (
                 TrackData& td,
                 const scalar dt,
-                const label
+                const label cellI
             );
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
index 5e8c48025e3..b3d6fd0b4a6 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
@@ -77,11 +77,11 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
     ReactingCloud<ParcelType>& owner,
     const label typeId,
     const vector& position,
-    const label celli,
+    const label cellI,
     const scalar d0,
     const vector& U0,
     const scalar nParticle0,
-    const scalarField& YMixture0,
+    const scalarField& Y0,
     const constantProperties& constProps
 )
 :
@@ -90,14 +90,14 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
         owner,
         typeId,
         position,
-        celli,
+        cellI,
         d0,
         U0,
         nParticle0,
         constProps
     ),
     mass0_(0.0),
-    YMixture_(YMixture0),
+    Y_(Y0),
     pc_(0.0)
 {
     // Set initial parcel mass
@@ -176,9 +176,9 @@ inline Foam::scalar Foam::ReactingParcel<ParcelType>::mass0() const
 
 template<class ParcelType>
 inline const Foam::scalarField& Foam::ReactingParcel<ParcelType>::
-YMixture() const
+Y() const
 {
-    return YMixture_;
+    return Y_;
 }
 
 
@@ -197,9 +197,9 @@ inline Foam::scalar& Foam::ReactingParcel<ParcelType>::mass0()
 
 
 template<class ParcelType>
-inline Foam::scalarField& Foam::ReactingParcel<ParcelType>::YMixture()
+inline Foam::scalarField& Foam::ReactingParcel<ParcelType>::Y()
 {
-    return YMixture_;
+    return Y_;
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
index 98076d4b27d..585f688eef2 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
@@ -39,7 +39,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
 :
     ThermoParcel<ParcelType>(cloud, is, readFields),
     mass0_(0.0),
-    YMixture_(0),
+    Y_(0),
     pc_(0.0)
 {
     if (readFields)
@@ -48,11 +48,11 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
             dynamic_cast<const ReactingCloud<ParcelType>& >(cloud);
 
         const label nMixture = cR.composition().phaseTypes().size();
-        YMixture_.setSize(nMixture);
+        Y_.setSize(nMixture);
 
         if (is.format() == IOstream::ASCII)
         {
-            is >> mass0_ >> YMixture_;
+            is >> mass0_ >> Y_;
         }
         else
         {
@@ -61,7 +61,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
                 reinterpret_cast<char*>(&mass0_),
               + sizeof(mass0_)
             );
-            is >> YMixture_;
+            is >> Y_;
         }
     }
 
@@ -109,13 +109,13 @@ void Foam::ReactingParcel<ParcelType>::readFields
     forAllIter(typename Cloud<ParcelType>, c, iter)
     {
         ReactingParcel<ParcelType>& p = iter();
-        p.YMixture_.setSize(nPhases, 0.0);
+        p.Y_.setSize(nPhases, 0.0);
     }
 
-    // Populate YMixture for each parcel
+    // Populate Y for each parcel
     forAll(phaseTypes, j)
     {
-        IOField<scalar> YMixture
+        IOField<scalar> Y
         (
             c.fieldIOobject("Y" + phaseTypes[j], IOobject::MUST_READ)
         );
@@ -124,7 +124,7 @@ void Foam::ReactingParcel<ParcelType>::readFields
         forAllIter(typename Cloud<ParcelType>, c, iter)
         {
             ReactingParcel<ParcelType>& p = iter();
-            p.YMixture_[j] = YMixture[i++];
+            p.Y_[j] = Y[i++];
         }
     }
 }
@@ -157,7 +157,7 @@ void Foam::ReactingParcel<ParcelType>::writeFields
         const wordList phaseTypes = c.composition().phaseTypes();
         forAll(phaseTypes, j)
         {
-            IOField<scalar> YMixture
+            IOField<scalar> Y
             (
                 c.fieldIOobject("Y" + phaseTypes[j], IOobject::NO_READ),
                 np
@@ -167,10 +167,10 @@ void Foam::ReactingParcel<ParcelType>::writeFields
             forAllConstIter(typename Cloud<ParcelType>, c, iter)
             {
                 const ReactingParcel<ParcelType>& p0 = iter();
-                YMixture[i++] = p0.YMixture()[j];
+                Y[i++] = p0.Y()[j];
             }
 
-            YMixture.write();
+            Y.write();
         }
     }
 }
@@ -189,7 +189,7 @@ Foam::Ostream& Foam::operator<<
     {
         os  << static_cast<const ThermoParcel<ParcelType>& >(p)
             << token::SPACE << p.mass0()
-            << token::SPACE << p.YMixture();
+            << token::SPACE << p.Y();
     }
     else
     {
@@ -200,7 +200,7 @@ Foam::Ostream& Foam::operator<<
             reinterpret_cast<const char*>(&p.mass0_),
             sizeof(p.mass0())
         );
-        os << p.YMixture();
+        os << p.Y();
     }
 
     // Check state of Ostream
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 16854519d04..99bf172920f 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -34,23 +34,23 @@ void Foam::ThermoParcel<ParcelType>::updateCellQuantities
 (
     TrackData& td,
     const scalar dt,
-    const label celli
+    const label cellI
 )
 {
-    KinematicParcel<ParcelType>::updateCellQuantities(td, dt, celli);
+    KinematicParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
 
-    Tc_ = td.TInterp().interpolate(this->position(), celli);
-    cpc_ = td.cpInterp().interpolate(this->position(), celli);
+    Tc_ = td.TInterp().interpolate(this->position(), cellI);
+    cpc_ = td.cpInterp().interpolate(this->position(), cellI);
 }
 
 
 template<class ParcelType>
 template<class TrackData>
-void Foam::ThermoParcel<ParcelType>::calcCoupled
+void Foam::ThermoParcel<ParcelType>::calc
 (
     TrackData& td,
     const scalar dt,
-    const label celli
+    const label cellI
 )
 {
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -59,8 +59,6 @@ void Foam::ThermoParcel<ParcelType>::calcCoupled
     const vector U0 = this->U_;
     const scalar mass0 = this->mass();
     const scalar np0 = this->nParticle_;
-//    const scalar T0 = T_;
-//    const scalar cp0 = cp_;
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -85,25 +83,27 @@ void Foam::ThermoParcel<ParcelType>::calcCoupled
     // Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    const scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
+    const scalar T1 = calcHeatTransfer(td, dt, cellI, htc, dhTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~
-    // Accumulate source terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Update momentum transfer
-    td.cloud().UTrans()[celli] += np0*dUTrans;
+    if (td.cloud().coupled())
+    {
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // Accumulate carrier phase source terms
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Accumulate coefficient to be applied in carrier phase momentum coupling
-    td.cloud().UCoeff()[celli] += np0*mass0*Cud;
+        // Update momentum transfer
+        td.cloud().UTrans()[cellI] += np0*dUTrans;
 
-    // Update enthalpy transfer
-    td.cloud().hTrans()[celli] += np0*dhTrans;
+        // Coefficient to be applied in carrier phase momentum coupling
+        td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
 
-    // Accumulate coefficient to be applied in carrier phase enthalpy coupling
-    td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
+        // Update enthalpy transfer
+        td.cloud().hTrans()[cellI] += np0*dhTrans;
 
+        // Coefficient to be applied in carrier phase enthalpy coupling
+        td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
+    }
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Set new particle properties
@@ -113,48 +113,13 @@ void Foam::ThermoParcel<ParcelType>::calcCoupled
 }
 
 
-template<class ParcelType>
-template<class TrackData>
-void Foam::ThermoParcel<ParcelType>::calcUncoupled
-(
-    TrackData& td,
-    const scalar dt,
-    const label celli
-)
-{
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Initialise transfer terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Momentum transfer from the particle to the carrier phase
-    vector dUTrans = vector::zero;
-
-    // Enthalpy transfer from the particle to the carrier phase
-    scalar dhTrans = 0.0;
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate velocity - update U
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Cud = 0.0;
-    this->U_ = calcVelocity(td, dt, Cud, dUTrans);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate heat transfer - update T
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar htc = 0.0;
-    T_ = calcHeatTransfer(td, dt, celli, htc, dhTrans);
-}
-
-
 template<class ParcelType>
 template <class TrackData>
 Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 (
     TrackData& td,
     const scalar dt,
-    const label celli,
+    const label cellI,
     scalar& htc,
     scalar& dhTrans
 )
@@ -193,7 +158,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
         const scalar sigma = radiation::sigmaSB.value();
         const scalar epsilon = td.constProps().epsilon0();
         const scalar epsilonSigmaT3 = epsilon*sigma*pow3(T_);
-        ap = (htc*Tc_ + 0.25*epsilon*G[celli])/(htc + epsilonSigmaT3);
+        ap = (htc*Tc_ + 0.25*epsilon*G[cellI])/(htc + epsilonSigmaT3);
         bp += epsilonSigmaT3;
     }
     bp *= 6.0/(this->rho_*this->d_*cp_);
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index 0f2b48416d0..bc66a2a7d2c 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -209,7 +209,7 @@ protected:
         (
             TrackData& td,
             const scalar dt,
-            const label celli,
+            const label cellI,
             scalar& htc,
             scalar& dhTrans
         );
@@ -231,7 +231,7 @@ public:
             ThermoCloud<ParcelType>& owner,
             const label typeId,
             const vector& position,
-            const label celli,
+            const label cellI,
             const scalar d0,
             const vector& U0,
             const scalar nParticle0,
@@ -281,25 +281,16 @@ public:
             (
                 TrackData& td,
                 const scalar dt,
-                const label celli
+                const label cellI
             );
 
-            //- Coupled calculation with the continuous phase
+            //- Update parcel properties over the time interval
             template<class TrackData>
-            void calcCoupled
+            void calc
             (
                 TrackData& td,
                 const scalar dt,
-                const label celli
-            );
-
-            //- Uncoupled calculation with the continuous phase
-            template<class TrackData>
-            void calcUncoupled
-            (
-                TrackData& td,
-                const scalar dt,
-                const label
+                const label cellI
             );
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
index f1f6e1d8fde..c6f0f2b816f 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H
@@ -75,7 +75,7 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
     ThermoCloud<ParcelType>& owner,
     const label typeId,
     const vector& position,
-    const label celli,
+    const label cellI,
     const scalar d0,
     const vector& U0,
     const scalar nParticle0,
@@ -87,7 +87,7 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
         owner,
         typeId,
         position,
-        celli,
+        cellI,
         d0,
         U0,
         nParticle0,
diff --git a/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.C b/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.C
index 244f49ad907..c0d438df24c 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.C
@@ -43,7 +43,7 @@ Foam::basicKinematicParcel::basicKinematicParcel
     KinematicCloud<basicKinematicParcel>& owner,
     const label typeId,
     const vector& position,
-    const label celli,
+    const label cellI,
     const scalar d0,
     const vector& U0,
     const scalar nParticle0,
@@ -55,7 +55,7 @@ Foam::basicKinematicParcel::basicKinematicParcel
         owner,
         typeId,
         position,
-        celli,
+        cellI,
         d0,
         U0,
         nParticle0,
diff --git a/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.H b/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.H
index decc0b57e83..b5da2267e77 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.H
@@ -66,7 +66,7 @@ public:
             KinematicCloud<basicKinematicParcel>& owner,
             const label typeId,
             const vector& position,
-            const label celli,
+            const label cellI,
             const scalar d0,
             const vector& U0,
             const scalar nParticle0,
diff --git a/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.C
index 99d4ec05b0a..0f7ce9173f8 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.C
@@ -43,14 +43,14 @@ Foam::basicReactingMultiphaseParcel::basicReactingMultiphaseParcel
     ReactingMultiphaseCloud<basicReactingMultiphaseParcel>& owner,
     const label typeId,
     const vector& position,
-    const label celli,
+    const label cellI,
     const scalar d0,
     const vector& U0,
     const scalar nParticle0,
     const scalarField& YGas0,
     const scalarField& YLiquid0,
     const scalarField& YSolid0,
-    const scalarField& YMixture0,
+    const scalarField& Y0,
     const constantProperties& constProps
 )
 :
@@ -59,14 +59,14 @@ Foam::basicReactingMultiphaseParcel::basicReactingMultiphaseParcel
         owner,
         typeId,
         position,
-        celli,
+        cellI,
         d0,
         U0,
         nParticle0,
         YGas0,
         YLiquid0,
         YSolid0,
-        YMixture0,
+        Y0,
         constProps
     )
 {}
diff --git a/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.H
index c39fae262ee..bb5053ff748 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.H
@@ -66,14 +66,14 @@ public:
              ReactingMultiphaseCloud<basicReactingMultiphaseParcel>& owner,
              const label typeId,
              const vector& position,
-             const label celli,
+             const label cellI,
              const scalar d0,
              const vector& U0,
              const scalar nParticle0,
              const scalarField& YGas0,
              const scalarField& YLiquid0,
              const scalarField& YSolid0,
-             const scalarField& YMixture0,
+             const scalarField& Y0,
              const constantProperties& constProps
         );
 
diff --git a/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.C b/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.C
index bb3a5e5da6b..722c572e433 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.C
@@ -43,11 +43,11 @@ Foam::basicReactingParcel::basicReactingParcel
     ReactingCloud<basicReactingParcel>& owner,
     const label typeId,
     const vector& position,
-    const label celli,
+    const label cellI,
     const scalar d0,
     const vector& U0,
     const scalar nParticle0,
-    const scalarField& YMixture0,
+    const scalarField& Y0,
     const constantProperties& constProps
 )
 :
@@ -56,11 +56,11 @@ Foam::basicReactingParcel::basicReactingParcel
         owner,
         typeId,
         position,
-        celli,
+        cellI,
         d0,
         U0,
         nParticle0,
-        YMixture0,
+        Y0,
         constProps
     )
 {}
diff --git a/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.H b/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.H
index 5deef5d7291..008a33b25bc 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.H
@@ -66,11 +66,11 @@ public:
             ReactingCloud<basicReactingParcel>& owner,
             const label typeId,
             const vector& position,
-            const label celli,
+            const label cellI,
             const scalar d0,
             const vector& U0,
             const scalar nParticle0,
-            const scalarField& YMixture0,
+            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 edbeef83b38..76a575d1dec 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.C
@@ -43,7 +43,7 @@ Foam::basicThermoParcel::basicThermoParcel
     ThermoCloud<basicThermoParcel>& owner,
     const label typeId,
     const vector position,
-    const label celli,
+    const label cellI,
     const scalar d0,
     const vector U0,
     const scalar nParticle0,
@@ -55,7 +55,7 @@ Foam::basicThermoParcel::basicThermoParcel
         owner,
         typeId,
         position,
-        celli,
+        cellI,
         d0,
         U0,
         nParticle0,
diff --git a/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.H b/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.H
index 0eb00002df4..c55471d3d14 100644
--- a/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.H
@@ -65,7 +65,7 @@ public:
             ThermoCloud<basicThermoParcel>& owner,
             const label typeId,
             const vector position,
-            const label celli,
+            const label cellI,
             const scalar d0,
             const vector U0,
             const scalar nParticle0,
diff --git a/src/lagrangian/intermediate/submodels/IO/DataEntry/polynomial/polynomial.C b/src/lagrangian/intermediate/submodels/IO/DataEntry/polynomial/polynomial.C
index 56101256e0a..276d54146be 100644
--- a/src/lagrangian/intermediate/submodels/IO/DataEntry/polynomial/polynomial.C
+++ b/src/lagrangian/intermediate/submodels/IO/DataEntry/polynomial/polynomial.C
@@ -26,6 +26,13 @@ License
 
 #include "polynomial.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(polynomial, 0);
+}
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::polynomial::polynomial(const word& entryName, Istream& is)
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
index b72c34f9bcf..e3b9383ecd4 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
@@ -67,6 +67,26 @@ Foam::label Foam::LiquidEvaporation<CloudType>::carrierSpecieId
 }
 
 
+template<class CloudType>
+Foam::scalarField Foam::LiquidEvaporation<CloudType>::calcXc
+(
+    const label cellI
+) const
+{
+    scalarField Xc(this->owner().carrierThermo().composition().Y().size());
+
+    scalar Winv = 0.0;
+    forAll(Xc, i)
+    {
+        scalar Y = this->owner().carrierThermo().composition().Y()[i][cellI];
+        Winv += Y/this->owner().gases()[i].W();
+        Xc[i] = Y/this->owner().gases()[i].W();
+    }
+
+    return Xc/Winv;
+}
+
+
 template <class CloudType>
 Foam::scalar Foam::LiquidEvaporation<CloudType>::Sh
 (
@@ -99,10 +119,10 @@ Foam::LiquidEvaporation<CloudType>::LiquidEvaporation
         )
     ),
     Tvap_(readScalar(this->coeffDict().lookup("Tvap"))),
-    evapProps_(this->coeffDict().lookup("activeLiquids")),
-    liqToGasMap_(evapProps_.size())
+    activeLiquids_(this->coeffDict().lookup("activeLiquids")),
+    liqToGasMap_(activeLiquids_.size())
 {
-    if (evapProps_.size() == 0)
+    if (activeLiquids_.size() == 0)
     {
         WarningIn
         (
@@ -116,9 +136,9 @@ Foam::LiquidEvaporation<CloudType>::LiquidEvaporation
     }
 
     // Calculate mapping between liquid and carrier phase species
-    forAll(evapProps_, i)
+    forAll(activeLiquids_, i)
     {
-        liqToGasMap_[i] = carrierSpecieId(evapProps_[i].name());
+        liqToGasMap_[i] = carrierSpecieId(activeLiquids_[i]);
     }
 }
 
@@ -142,15 +162,15 @@ bool Foam::LiquidEvaporation<CloudType>::active() const
 template<class CloudType>
 Foam::scalar Foam::LiquidEvaporation<CloudType>::calculate
 (
+    const scalar dt,
+    const label cellI,
     const scalar T,
+    const scalar pc,
     const scalar d,
-    const scalarField& Xc,
-    scalarList& dMassMT,
-    const vector& Ur,
     const scalar Tc,
-    const scalar pc,
     const scalar nuc,
-    const scalar dt
+    const vector& Ur,
+    scalarList& dMassMT
 ) const
 {
     scalar dMassTot = 0.0;
@@ -162,20 +182,25 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::calculate
     }
     else
     {
-        // droplet area
+        // Construct carrier phase species volume fractions
+        scalarField Xc = calcXc(cellI);
+
+        // droplet surface area
         scalar A = mathematicalConstant::pi*sqr(d);
 
         // Reynolds number
         scalar Re = mag(Ur)*d/(nuc + ROOTVSMALL);
 
         // Calculate mass transfer of each specie in liquid
-        forAll(evapProps_, i)
+        forAll(activeLiquids_, i)
         {
-            // Diffusion coefficient for species i
-            scalar Dab = evapProps_[i].Dab();
+            label gid = liqToGasMap_[i];
+
+            // Vapour diffusivity [m2/s]
+            scalar Dab = liquids_->properties()[gid].D(pc, T);
 
-            // Saturation pressure for species i at temperature T
-            scalar pSat = evapProps_[i].TvsPSat().value(T);
+            // Saturation pressure for species i [pa]
+            scalar pSat = liquids_->properties()[gid].pv(pc, T);
 
             // Schmidt number
             scalar Sc = nuc/(Dab + ROOTVSMALL);
@@ -190,15 +215,14 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::calculate
             scalar Cs = pSat/(specie::RR*T);
 
             // vapour concentration in bulk gas [kgmol/m3]
-            scalar Cinf = Xc[i]*pc/(specie::RR*Tc);
+            scalar Cinf = Xc[gid]*pc/(specie::RR*Tc);
 
             // molar flux of vapour [kgmol/m2/s]
             scalar Ni = max(kc*(Cs - Cinf), 0.0);
 
-            // mass transfer
-            label globalLiqId = liqToGasMap_[i];
-            scalar dm = Ni*A*liquids_->properties()[globalLiqId].W()*dt;
-            dMassMT[globalLiqId] -= dm;
+            // mass transfer [kg]
+            scalar dm = Ni*A*liquids_->properties()[gid].W()*dt;
+            dMassMT[gid] -= dm;
             dMassTot += dm;
         }
     }
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
index 0224b1bd5da..85d3714df34 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
@@ -35,9 +35,6 @@ Description
 
 #include "PhaseChangeModel.H"
 #include "liquidMixture.H"
-#include "Tuple2.H"
-
-#include "evaporationProperties.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,8 +61,8 @@ protected:
         //  allowed
         scalar Tvap_;
 
-        //- List of evaporation properties
-        List<evaporationProperties> evapProps_;
+        //- List of active liquid names
+        List<word> activeLiquids_;
 
         //- Mapping between liquid and carrier gas species
         List<label> liqToGasMap_;
@@ -79,6 +76,9 @@ protected:
         //- Sherwood number as a function of Reynolds and Schmidt numbers
         scalar Sh(const scalar Re, const scalar Sc) const;
 
+        //- Calculate the carrier phase comonent volume fractions at cellI
+        scalarField calcXc(const label cellI) const;
+
 
 public:
 
@@ -108,15 +108,15 @@ public:
         //- Update model
         scalar calculate
         (
+            const scalar dt,
+            const label cellI,
             const scalar T,
+            const scalar pc,
             const scalar d,
-            const scalarField& Xc,
-            scalarList& dMassMT,
-            const vector& Ur,
             const scalar Tc,
-            const scalar pc,
             const scalar nuc,
-            const scalar dt
+            const vector& Ur,
+            scalarList& dMassMT
         ) const;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C
index 360f0b39262..aa4fbc4cee3 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C
@@ -58,15 +58,15 @@ bool Foam::NoPhaseChange<CloudType>::active() const
 template<class CloudType>
 Foam::scalar Foam::NoPhaseChange<CloudType>::calculate
 (
+    const scalar,
+    const label,
     const scalar,
     const scalar,
-    const scalarField&,
-    scalarList&,
-    const vector&,
     const scalar,
     const scalar,
     const scalar,
-    const scalar
+    const vector&,
+    scalarList&
 ) 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 d360e36d09e..1a58a5141f7 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H
@@ -72,15 +72,15 @@ public:
         //- Update model
         scalar calculate
         (
+            const scalar dt,
+            const label cellI,
             const scalar T,
+            const scalar pc,
             const scalar d,
-            const scalarField& Xc,
-            scalarList& dMassMT,
-            const vector& Ur,
             const scalar Tc,
-            const scalar pc,
             const scalar nuc,
-            const scalar dt
+            const vector& Ur,
+            scalarList& dMassMT
         ) const;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
index d0d40badc47..ff9a77e9b68 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
@@ -138,15 +138,15 @@ public:
         //- Update model
         virtual scalar calculate
         (
+            const scalar dt,
+            const label cellI,
             const scalar T,
+            const scalar pc,
             const scalar d,
-            const scalarField& Xc,
-            scalarList& dMassMT,
-            const vector& Ur,
             const scalar Tc,
-            const scalar pc,
             const scalar nuc,
-            const scalar dt
+            const vector& Ur,
+            scalarList& dMassMT
         ) const = 0;
 };
 
-- 
GitLab