diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
index 9a9830545e..d765e9f60b 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
@@ -97,6 +97,16 @@ void Foam::polyMesh::clearAddressing()
     geometricD_ = Vector<label>::zero;
     solutionD_ = Vector<label>::zero;
 
+    // Update zones
+    pointZones_.clearAddressing();
+    faceZones_.clearAddressing();
+    cellZones_.clearAddressing();
+
+    // Remove the stored tet base points
+    tetBasePtIsPtr_.clear();
+    // Remove the cell tree
+    cellTreePtr_.clear();
+
     pointMesh::Delete(*this);
 }
 
diff --git a/src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntryFwd.H b/src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntryFwd.H
new file mode 100644
index 0000000000..326a2a4259
--- /dev/null
+++ b/src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntryFwd.H
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DataEntryFws_H
+#define DataEntryFws_H
+
+#include "DataEntry.H"
+#include "vector.H"
+#include "symmTensor.H"
+#include "sphericalTensor.H"
+#include "tensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef DataEntry<label> labelDataEntry;
+    typedef DataEntry<scalar> scalarDataEntry;
+    typedef DataEntry<vector> vectorDataEntry;
+    typedef DataEntry<symmTensor> symmTensorDataEntry;
+    typedef DataEntry<sphericalTensor> sphericalTensorDataEntry;
+    typedef DataEntry<tensor> tensorDataEntry;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.C b/src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.C
index d015aede93..cfbaec7424 100644
--- a/src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.C
+++ b/src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.C
@@ -25,14 +25,14 @@ License
 
 #include "polynomial.H"
 #include "Time.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
     defineTypeNameAndDebug(polynomial, 0);
-    DataEntry<scalar>::adddictionaryConstructorToTable<polynomial>
-        addpolynomialConstructorToTable_;
+    addToRunTimeSelectionTable(scalarDataEntry, polynomial, dictionary);
 }
 
 
@@ -40,7 +40,7 @@ namespace Foam
 
 Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
 :
-    DataEntry<scalar>(entryName),
+    scalarDataEntry(entryName),
     coeffs_(),
     canIntegrate_(true),
     dimensions_(dimless)
@@ -52,15 +52,17 @@ Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
     is.putBack(firstToken);
     if (firstToken == token::BEGIN_SQR)
     {
-        is >> this->dimensions_;
+        is  >> this->dimensions_;
     }
 
     is  >> coeffs_;
 
     if (!coeffs_.size())
     {
-        FatalErrorIn("Foam::polynomial::polynomial(const word&, dictionary&)")
-            << "polynomial coefficients for entry " << this->name_
+        FatalErrorIn
+        (
+            "Foam::polynomial::polynomial(const word&, const dictionary&)"
+        )   << "polynomial coefficients for entry " << this->name_
             << " are invalid (empty)" << nl << exit(FatalError);
     }
 
@@ -77,8 +79,10 @@ Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
     {
         if (!canIntegrate_)
         {
-            WarningIn("Foam::polynomial::polynomial(const word&, dictionary&)")
-                << "Polynomial " << this->name_ << " cannot be integrated"
+            WarningIn
+            (
+                "Foam::polynomial::polynomial(const word&, const dictionary&)"
+            )   << "Polynomial " << this->name_ << " cannot be integrated"
                 << endl;
         }
     }
@@ -91,7 +95,7 @@ Foam::polynomial::polynomial
     const List<Tuple2<scalar, scalar> >& coeffs
 )
 :
-    DataEntry<scalar>(entryName),
+    scalarDataEntry(entryName),
     coeffs_(coeffs),
     canIntegrate_(true),
     dimensions_(dimless)
@@ -101,7 +105,7 @@ Foam::polynomial::polynomial
         FatalErrorIn
         (
             "Foam::polynomial::polynomial"
-            "(const word&, const List<Tuple2<scalar, scalar> >&&)"
+            "(const word&, const List<Tuple2<scalar, scalar> >&)"
         )   << "polynomial coefficients for entry " << this->name_
             << " are invalid (empty)" << nl << exit(FatalError);
     }
@@ -122,7 +126,7 @@ Foam::polynomial::polynomial
             WarningIn
             (
                 "Foam::polynomial::polynomial"
-                "(const word&, const List<Tuple2<scalar, scalar> >&&)"
+                "(const word&, const List<Tuple2<scalar, scalar> >&)"
             )   << "Polynomial " << this->name_ << " cannot be integrated"
                 << endl;
         }
@@ -132,7 +136,7 @@ Foam::polynomial::polynomial
 
 Foam::polynomial::polynomial(const polynomial& poly)
 :
-    DataEntry<scalar>(poly),
+    scalarDataEntry(poly),
     coeffs_(poly.coeffs_),
     canIntegrate_(poly.canIntegrate_),
     dimensions_(poly.dimensions_)
@@ -201,7 +205,8 @@ Foam::dimensioned<Foam::scalar> Foam::polynomial::dimValue
 
 Foam::dimensioned<Foam::scalar> Foam::polynomial::dimIntegrate
 (
-    const scalar x1, const scalar x2
+    const scalar x1,
+    const scalar x2
 ) const
 {
     return dimensioned<scalar>
@@ -212,4 +217,5 @@ Foam::dimensioned<Foam::scalar> Foam::polynomial::dimIntegrate
     );
 }
 
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.H b/src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.H
index dc4fc7986b..3e181c0787 100644
--- a/src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.H
+++ b/src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.H
@@ -48,6 +48,7 @@ SourceFiles
 #include "DataEntry.H"
 #include "Tuple2.H"
 #include "dimensionSet.H"
+#include "DataEntryFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -70,7 +71,7 @@ Ostream& operator<<
 
 class polynomial
 :
-    public DataEntry<scalar>
+    public scalarDataEntry
 {
     // Private data
 
@@ -107,9 +108,9 @@ public:
         polynomial(const polynomial& poly);
 
         //- Construct and return a clone
-        virtual tmp<DataEntry<scalar> > clone() const
+        virtual tmp<scalarDataEntry> clone() const
         {
-            return tmp<DataEntry<scalar> >(new polynomial(*this));
+            return tmp<scalarDataEntry>(new polynomial(*this));
         }
 
 
diff --git a/src/combustionModels/Make/files b/src/combustionModels/Make/files
index 20de91ecc6..88e81afe8b 100644
--- a/src/combustionModels/Make/files
+++ b/src/combustionModels/Make/files
@@ -10,6 +10,7 @@ rhoCombustionModel/rhoCombustionModel/rhoCombustionModelNew.C
 rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.C
 rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C
 
+diffusion/diffusions.C
 infinitelyFastChemistry/infinitelyFastChemistrys.C
 
 PaSR/PaSRs.C
diff --git a/src/combustionModels/PaSR/PaSR.C b/src/combustionModels/PaSR/PaSR.C
index cd12d896e0..d90eb68357 100644
--- a/src/combustionModels/PaSR/PaSR.C
+++ b/src/combustionModels/PaSR/PaSR.C
@@ -138,7 +138,7 @@ void Foam::combustionModels::PaSR<Type>::correct()
 
 template<class Type>
 Foam::tmp<Foam::fvScalarMatrix>
-Foam::combustionModels::PaSR<Type>::R(const volScalarField& Y) const
+Foam::combustionModels::PaSR<Type>::R(volScalarField& Y) const
 {
     tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
 
diff --git a/src/combustionModels/PaSR/PaSR.H b/src/combustionModels/PaSR/PaSR.H
index 4c4c3347e5..fba886f8ac 100644
--- a/src/combustionModels/PaSR/PaSR.H
+++ b/src/combustionModels/PaSR/PaSR.H
@@ -103,7 +103,7 @@ public:
             virtual void correct();
 
             //- Fuel consumption rate matrix.
-            virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
             //- Heat release rate calculated from fuel consumption rate matrix
             virtual tmp<volScalarField> dQ() const;
diff --git a/src/combustionModels/combustionModel/combustionModel.H b/src/combustionModels/combustionModel/combustionModel.H
index 36322adf84..46077de9e4 100644
--- a/src/combustionModels/combustionModel/combustionModel.H
+++ b/src/combustionModels/combustionModel/combustionModel.H
@@ -135,7 +135,7 @@ public:
         virtual void correct() = 0;
 
         //- Fuel consumption rate matrix, i.e. source term for fuel equation
-        virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const = 0;
+        virtual tmp<fvScalarMatrix> R(volScalarField& Y) const = 0;
 
         //- Heat release rate calculated from fuel consumption rate matrix
         virtual tmp<volScalarField> dQ() const = 0;
diff --git a/src/combustionModels/diffusion/diffusion.C b/src/combustionModels/diffusion/diffusion.C
new file mode 100644
index 0000000000..44e413dfc0
--- /dev/null
+++ b/src/combustionModels/diffusion/diffusion.C
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "diffusion.H"
+#include "fvcGrad.H"
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+diffusion<CombThermoType, ThermoType>::diffusion
+(
+    const word& modelType, const fvMesh& mesh
+)
+:
+    singleStepCombustion<CombThermoType, ThermoType>(modelType, mesh),
+    C_(readScalar(this->coeffs().lookup("C"))),
+    oxidantName_(this->coeffs().template lookupOrDefault<word>("oxidant", "O2"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+diffusion<CombThermoType, ThermoType>::~diffusion()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+void diffusion<CombThermoType, ThermoType>::correct()
+{
+    this->wFuel_ ==
+        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0);
+
+    if (this->active())
+    {
+        this->singleMixturePtr_->fresCorrect();
+
+        const label fuelI = this->singleMixturePtr_->fuelIndex();
+
+        const volScalarField& YFuel =
+            this->thermoPtr_->composition().Y()[fuelI];
+
+        if (this->thermoPtr_->composition().contains(oxidantName_))
+        {
+            const volScalarField& YO2 =
+                this->thermoPtr_->composition().Y(oxidantName_);
+
+            this->wFuel_ ==
+                C_*this->turbulence().muEff()
+               *mag(fvc::grad(YFuel) & fvc::grad(YO2))
+               *pos(YFuel)*pos(YO2);
+        }
+    }
+}
+
+
+template<class CombThermoType, class ThermoType>
+bool diffusion<CombThermoType, ThermoType>::read()
+{
+    if (singleStepCombustion<CombThermoType, ThermoType>::read())
+    {
+        this->coeffs().lookup("C") >> C_ ;
+        this->coeffs().readIfPresent("oxidant", oxidantName_);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/combustionModels/diffusion/diffusion.H b/src/combustionModels/diffusion/diffusion.H
new file mode 100644
index 0000000000..c3002a894f
--- /dev/null
+++ b/src/combustionModels/diffusion/diffusion.H
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::diffusion
+
+Description
+    Simple diffusion-based combustion model based on the principle mixed is
+    burnt. Additional parameter C is used to distribute the heat release rate
+    in time.
+
+SourceFiles
+    diffusion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef diffusion_H
+#define diffusion_H
+
+#include "singleStepCombustion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class diffusion Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CombThermoType, class ThermoType>
+class diffusion
+:
+    public singleStepCombustion<CombThermoType, ThermoType>
+{
+    // Private data
+
+        //- Model constant
+        scalar C_;
+
+        //- Name of oxidant - default is "O2"
+        word oxidantName_;
+
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        diffusion(const diffusion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const diffusion&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("diffusion");
+
+
+    // Constructors
+
+        //- Construct from components
+        diffusion(const word& modelType, const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~diffusion();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct combustion rate
+            virtual void correct();
+
+
+        // I-O
+
+            //- Update properties
+            virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "diffusion.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/combustionModels/diffusion/diffusions.C b/src/combustionModels/diffusion/diffusions.C
new file mode 100644
index 0000000000..73c201f0fc
--- /dev/null
+++ b/src/combustionModels/diffusion/diffusions.C
@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "thermoPhysicsTypes.H"
+#include "psiThermoCombustion.H"
+#include "rhoThermoCombustion.H"
+#include "diffusion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        psiThermoCombustion,
+        gasThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        psiThermoCombustion,
+        constGasThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        rhoThermoCombustion,
+        gasThermoPhysics,
+        rhoCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        rhoThermoCombustion,
+        constGasThermoPhysics,
+        rhoCombustionModel
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/combustionModels/noCombustion/noCombustion.C b/src/combustionModels/noCombustion/noCombustion.C
index 3abdedaf29..aa38aad70f 100644
--- a/src/combustionModels/noCombustion/noCombustion.C
+++ b/src/combustionModels/noCombustion/noCombustion.C
@@ -59,7 +59,7 @@ template<class CombThermoType>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::combustionModels::noCombustion<CombThermoType>::R
 (
-    const volScalarField& Y
+    volScalarField& Y
 ) const
 {
     tmp<fvScalarMatrix> tSu
diff --git a/src/combustionModels/noCombustion/noCombustion.H b/src/combustionModels/noCombustion/noCombustion.H
index 8e845d43f1..d8f8f39912 100644
--- a/src/combustionModels/noCombustion/noCombustion.H
+++ b/src/combustionModels/noCombustion/noCombustion.H
@@ -84,8 +84,8 @@ public:
             //- Correct combustion rate
             virtual void correct();
 
-            //- Fuel consumption rate matrix.
-            virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
+            //- Fuel consumption rate matrix
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
             //- Heat release rate calculated from fuel consumption rate matrix
             virtual tmp<volScalarField> dQ() const;
diff --git a/src/combustionModels/singleStepCombustion/singleStepCombustion.C b/src/combustionModels/singleStepCombustion/singleStepCombustion.C
index aab73f7abf..6d9b08c72d 100644
--- a/src/combustionModels/singleStepCombustion/singleStepCombustion.C
+++ b/src/combustionModels/singleStepCombustion/singleStepCombustion.C
@@ -53,7 +53,8 @@ singleStepCombustion<CombThermoType, ThermoType>::singleStepCombustion
         ),
         this->mesh(),
         dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
-    )
+    ),
+    semiImplicit_(readBool(this->coeffs_.lookup("semiImplicit")))
 {
     if (isA<singleStepReactingMixture<ThermoType> >(this->thermo()))
     {
@@ -79,6 +80,15 @@ singleStepCombustion<CombThermoType, ThermoType>::singleStepCombustion
             << "Please select a thermo package based on "
             << "singleStepReactingMixture" << exit(FatalError);
     }
+
+    if (semiImplicit_)
+    {
+        Info<< "Combustion mode: semi-implicit" << endl;
+    }
+    else
+    {
+        Info<< "Combustion mode: explicit" << endl;
+    }
 }
 
 
@@ -95,17 +105,28 @@ template<class CombThermoType, class ThermoType>
 Foam::tmp<Foam::fvScalarMatrix>
 singleStepCombustion<CombThermoType, ThermoType>::R
 (
-    const volScalarField& Y
+    volScalarField& Y
 ) const
 {
     const label specieI = this->thermoPtr_->composition().species()[Y.name()];
 
-    const volScalarField wSpecie
+    volScalarField wSpecie
     (
         wFuel_*singleMixturePtr_->specieStoichCoeffs()[specieI]
     );
 
-    return wSpecie + fvm::Sp(0.0*wSpecie, Y);
+    if (semiImplicit_)
+    {
+        const label fNorm = singleMixturePtr_->specieProd()[specieI];
+        const volScalarField fres = singleMixturePtr_->fres(specieI);
+        wSpecie /= max(fNorm*(Y - fres), 1e-2);
+
+        return -fNorm*wSpecie*fres + fNorm*fvm::Sp(wSpecie, Y);
+    }
+    else
+    {
+        return wSpecie + fvm::Sp(0.0*wSpecie, Y);
+    }
 }
 
 
@@ -114,7 +135,8 @@ Foam::tmp<Foam::volScalarField>
 singleStepCombustion< CombThermoType, ThermoType>::Sh() const
 {
     const label fuelI = singleMixturePtr_->fuelIndex();
-    const volScalarField& YFuel = this->thermoPtr_->composition().Y(fuelI);
+    volScalarField& YFuel =
+        const_cast<volScalarField&>(this->thermoPtr_->composition().Y(fuelI));
 
     return -singleMixturePtr_->qFuel()*(R(YFuel) & YFuel);
 }
diff --git a/src/combustionModels/singleStepCombustion/singleStepCombustion.H b/src/combustionModels/singleStepCombustion/singleStepCombustion.H
index fd743247b3..7d691396ae 100644
--- a/src/combustionModels/singleStepCombustion/singleStepCombustion.H
+++ b/src/combustionModels/singleStepCombustion/singleStepCombustion.H
@@ -71,6 +71,9 @@ protected:
         //- Fuel consumption rate
         volScalarField wFuel_;
 
+        //- Semi-implicit (true) or explicit (false) treatment
+        bool semiImplicit_;
+
 
 public:
 
@@ -89,7 +92,7 @@ public:
         // Evolution
 
             //- Fuel consumption rate matrix
-            virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
             //- Heat release rate calculated from fuel consumption rate matrix
             virtual tmp<volScalarField> dQ() const;
diff --git a/src/lagrangian/coalCombustion/Make/files b/src/lagrangian/coalCombustion/Make/files
index 019781420d..f3d07eaa1f 100644
--- a/src/lagrangian/coalCombustion/Make/files
+++ b/src/lagrangian/coalCombustion/Make/files
@@ -1,4 +1,5 @@
 /* Coal parcel and sub-models */
 coalParcel/makeCoalParcelSubmodels.C
+coalCloudList/coalCloudList.C
 
 LIB = $(FOAM_LIBBIN)/libcoalCombustion
diff --git a/src/lagrangian/coalCombustion/coalCloudList/coalCloudList.C b/src/lagrangian/coalCombustion/coalCloudList/coalCloudList.C
new file mode 100644
index 0000000000..548e81d216
--- /dev/null
+++ b/src/lagrangian/coalCombustion/coalCloudList/coalCloudList.C
@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "coalCloudList.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::coalCloudList::coalCloudList
+(
+    const volScalarField& rho,
+    const volVectorField& U,
+    const dimensionedVector& g,
+    const SLGThermo& slgThermo
+)
+:
+    PtrList<coalCloud>(),
+    mesh_(rho.mesh())
+{
+    IOdictionary props
+    (
+        IOobject
+        (
+            "coalCloudList",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ
+        )
+    );
+
+    const wordHashSet cloudNames(wordList(props.lookup("clouds")));
+
+    setSize(cloudNames.size());
+
+    label i = 0;
+    forAllConstIter(wordHashSet, cloudNames, iter)
+    {
+        const word& name = iter.key();
+
+        Info<< "creating cloud: " << name << endl;
+
+        set
+        (
+            i++,
+            new coalCloud
+            (
+                name,
+                rho,
+                U,
+                g,
+                slgThermo
+            )
+        );
+
+        Info<< endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::coalCloudList::evolve()
+{
+    forAll(*this, i)
+    {
+        operator[](i).evolve();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/coalCombustion/coalCloudList/coalCloudList.H b/src/lagrangian/coalCombustion/coalCloudList/coalCloudList.H
new file mode 100644
index 0000000000..1a023e27ba
--- /dev/null
+++ b/src/lagrangian/coalCombustion/coalCloudList/coalCloudList.H
@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef coalCloudList_H
+#define coalCloudList_H
+
+#include "coalCloud.H"
+#include "volFieldsFwd.H"
+#include "fvMatricesFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class coalCloudList Declaration
+\*---------------------------------------------------------------------------*/
+
+class coalCloudList
+:
+    public PtrList<coalCloud>
+{
+private:
+
+    //- Reference to the mesh
+    const fvMesh& mesh_;
+
+
+public:
+
+    // Constructor
+
+        coalCloudList
+        (
+            const volScalarField& rho,
+            const volVectorField& U,
+            const dimensionedVector& g,
+            const SLGThermo& slgThermo
+        );
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Evolve the cloud collection
+            void evolve();
+
+
+        // Source terms
+
+            //- Return const reference to momentum source
+            inline tmp<DimensionedField<vector, volMesh> > UTrans() const;
+
+            //- Return tmp momentum source term
+            inline tmp<fvVectorMatrix> SU(volVectorField& U) const;
+
+            //- Sensible enthalpy transfer [J/kg]
+            inline tmp<DimensionedField<scalar, volMesh> > hsTrans() const;
+
+            //- Return sensible enthalpy source term [J/kg/m3/s]
+            inline tmp<fvScalarMatrix> Sh(volScalarField& hs) const;
+
+
+            //- Return mass source term for specie i - specie eqn
+            inline tmp<fvScalarMatrix> SYi
+            (
+                const label i,
+                volScalarField& Yi
+            ) const;
+
+            //- Return total mass transfer [kg/m3]
+            inline tmp<DimensionedField<scalar, volMesh> > rhoTrans() const;
+
+            //- Return tmp total mass source for carrier phase
+            //  - fully explicit
+            inline tmp<DimensionedField<scalar, volMesh> > Srho() const;
+
+            //- Return tmp total mass source for carrier phase specie i
+            //  - fully explicit
+            inline tmp<DimensionedField<scalar, volMesh> > Srho
+            (
+                const label i
+            ) const;
+
+            //- Return total mass source term [kg/m3/s]
+            inline tmp<fvScalarMatrix> Srho(volScalarField& rho) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "coalCloudListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H b/src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H
new file mode 100644
index 0000000000..005d9ef604
--- /dev/null
+++ b/src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H
@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvMatrices.H"
+#include "volFields.H"
+#include "DimensionedField.H"
+
+Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
+Foam::coalCloudList::UTrans() const
+{
+    tmp<DimensionedField<vector, volMesh> > tfld
+    (
+        new DimensionedField<vector, volMesh>
+        (
+            IOobject
+            (
+                "UTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedVector("zero", dimMass*dimVelocity, vector::zero)
+        )
+    );
+
+    DimensionedField<vector, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        fld += operator[](i).UTrans();
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix> Foam::coalCloudList::SU
+(
+    volVectorField& U
+) const
+{
+    tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
+    fvVectorMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).SU(U);
+    }
+
+    return tfvm;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::hsTrans() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "hsTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimEnergy, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        fld += operator[](i).hsTrans();
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Sh
+(
+    volScalarField& hs
+) const
+{
+    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(hs, dimEnergy/dimTime));
+    fvScalarMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).Sh(hs);
+    }
+
+    return tfvm;
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::SYi
+(
+    const label ii,
+    volScalarField& Yi
+) const
+{
+    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(Yi, dimMass/dimTime));
+    fvScalarMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).SYi(ii, Yi);
+    }
+
+    return tfvm;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::rhoTrans() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "rhoTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimMass, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        forAll(operator[](i).rhoTrans(), j)
+        {
+            fld += operator[](i).rhoTrans()[j];
+        }
+    }
+
+    return tfld;
+}
+
+
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::Srho() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "rhoTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimDensity/dimTime, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        fld += operator[](i).Srho();
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::Srho
+(
+    const label i
+) const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "rhoTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimDensity/dimTime, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, j)
+    {
+        fld += operator[](j).Srho(i);
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Srho
+(
+    volScalarField& rho
+) const
+{
+    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(rho, dimMass/dimTime));
+    fvScalarMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).Srho(rho);
+    }
+
+    return tfvm;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H b/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H
index 37627758ff..2df26d9a3c 100644
--- a/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H
+++ b/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,7 @@ License
 #include "NoSurfaceReaction.H"
 #include "COxidationDiffusionLimitedRate.H"
 #include "COxidationKineticDiffusionLimitedRate.H"
+#include "COxidationHurtMitchell.H"
 #include "COxidationMurphyShaddix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -43,6 +44,7 @@ License
         COxidationKineticDiffusionLimitedRate,                                \
         CloudType                                                             \
     );                                                                        \
+    makeSurfaceReactionModelType(COxidationHurtMitchell, CloudType);          \
     makeSurfaceReactionModelType(COxidationMurphyShaddix, CloudType);
 
 
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.C b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.C
new file mode 100644
index 0000000000..c3ecbbffaf
--- /dev/null
+++ b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.C
@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "COxidationHurtMitchell.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::COxidationHurtMitchell<CloudType>::COxidationHurtMitchell
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    SurfaceReactionModel<CloudType>(dict, owner, typeName),
+    Sb_(readScalar(this->coeffDict().lookup("Sb"))),
+    CsLocalId_(-1),
+    ashLocalId_(-1),
+    O2GlobalId_(owner.composition().globalCarrierId("O2")),
+    CO2GlobalId_(owner.composition().globalCarrierId("CO2")),
+    WC_(0.0),
+    WO2_(0.0),
+    HcCO2_(0.0),
+    heatOfReaction_(-1.0)
+{
+    // Determine Cs and ash ids
+    label idSolid = owner.composition().idSolid();
+    CsLocalId_ = owner.composition().localId(idSolid, "C");
+    ashLocalId_ = owner.composition().localId(idSolid, "ash", true);
+
+    // Set local copies of thermo properties
+    WO2_ = owner.thermo().carrier().W(O2GlobalId_);
+    const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
+    WC_ = WCO2 - WO2_;
+
+    HcCO2_ = owner.thermo().carrier().Hc(CO2GlobalId_);
+
+    const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
+    const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
+    Info<< "    C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
+
+    if (this->coeffDict().readIfPresent("heatOfReaction", heatOfReaction_))
+    {
+        Info<< "    Using user specified heat of reaction: "
+            << heatOfReaction_ << " [J/kg]" << endl;
+    }
+}
+
+
+template<class CloudType>
+Foam::COxidationHurtMitchell<CloudType>::COxidationHurtMitchell
+(
+    const COxidationHurtMitchell<CloudType>& srm
+)
+:
+    SurfaceReactionModel<CloudType>(srm),
+    Sb_(srm.Sb_),
+    CsLocalId_(srm.CsLocalId_),
+    ashLocalId_(srm.ashLocalId_),
+    O2GlobalId_(srm.O2GlobalId_),
+    CO2GlobalId_(srm.CO2GlobalId_),
+    WC_(srm.WC_),
+    WO2_(srm.WO2_),
+    HcCO2_(srm.HcCO2_),
+    heatOfReaction_(srm.heatOfReaction_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::COxidationHurtMitchell<CloudType>::~COxidationHurtMitchell()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::COxidationHurtMitchell<CloudType>::calculate
+(
+    const scalar dt,
+    const label cellI,
+    const scalar d,
+    const scalar T,
+    const scalar Tc,
+    const scalar pc,
+    const scalar rhoc,
+    const scalar mass,
+    const scalarField& YGas,
+    const scalarField& YLiquid,
+    const scalarField& YSolid,
+    const scalarField& YMixture,
+    const scalar N,
+    scalarField& dMassGas,
+    scalarField& dMassLiquid,
+    scalarField& dMassSolid,
+    scalarField& dMassSRCarrier
+) const
+{
+    const label idGas = CloudType::parcelType::GAS;
+    const label idSolid = CloudType::parcelType::SLD;
+    const scalar Ychar = YMixture[idSolid]*YSolid[CsLocalId_];
+
+    // Surface combustion until combustible fraction is consumed
+    if (Ychar < SMALL)
+    {
+        return 0.0;
+    }
+
+    const SLGThermo& thermo = this->owner().thermo();
+
+    // Local mass fraction of O2 in the carrier phase
+    const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
+
+    // No combustion if no oxygen present
+    if (YO2 < SMALL)
+    {
+        return 0.0;
+    }
+
+    // Conversion from [g/cm^2) to [kg/m^2]
+    const scalar convSI = 1000.0/10000.0;
+
+    // Universal gas constant in [kcal/mol/K]
+    const scalar RRcal = 1985.877534;
+
+    // Dry mass fraction
+    scalar Ydaf = YMixture[idGas] + YMixture[idSolid];
+    if (ashLocalId_ != -1)
+    {
+        Ydaf -= YMixture[idSolid]*YSolid[ashLocalId_];
+    }
+
+    // Char percentage
+    const scalar charPrc = Ychar/Ydaf*100.0;
+
+    // Particle surface area
+    const scalar Ap = constant::mathematical::pi*sqr(d);
+
+    // Far field partial pressure O2 [Pa]
+    // Note: Should really use the surface partial pressure
+    const scalar ppO2 = max(0.0, rhoc*YO2/WO2_*specie::RR*Tc);
+
+    // Activation energy [kcal/mol]
+    const scalar E = -5.94 + 0.355*charPrc;
+
+    // Pre-exponential factor [g/(cm^2.s.atm^0.5)]
+    const scalar lnK1750 = 2.8 - 0.0758*charPrc;
+    const scalar A = exp(lnK1750 + E/RRcal/1750.0);
+
+    // Kinetic rate of char oxidation [g/(cm^2.s.atm^0.5)]
+    const scalar Rk = A*exp(-E/(RRcal*T));
+
+    // Molar reaction rate per unit surface area [kmol/(m^2.s)]
+    const scalar qCsLim = mass*Ychar/(WC_*Ap*dt);
+    const scalar qCs = min(convSI*Rk*Foam::sqrt(ppO2/101325.0), qCsLim);
+
+    // Calculate the number of molar units reacted [kmol]
+    const scalar dOmega = qCs*Ap*dt;
+
+    // Add to carrier phase mass transfer
+    dMassSRCarrier[O2GlobalId_] += -dOmega*Sb_*WO2_;
+    dMassSRCarrier[CO2GlobalId_] += dOmega*(WC_ + Sb_*WO2_);
+
+    // Add to particle mass transfer
+    dMassSolid[CsLocalId_] += dOmega*WC_;
+
+
+    // Return the heat of reaction [J]
+    // note: carrier sensible enthalpy exchange handled via change in mass
+    if (heatOfReaction_ < 0)
+    {
+        const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
+        return dOmega*(WC_*HsC - (WC_ + Sb_*WO2_)*HcCO2_);
+    }
+    else
+    {
+        return dOmega*WC_*heatOfReaction_;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.H b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.H
new file mode 100644
index 0000000000..cd7800328f
--- /dev/null
+++ b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.H
@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    COxidationHurtMitchell
+
+Description
+    Char oxidation model given by Hurt and Mitchell:
+
+    Based on the reference:
+        Hurt R. and Mitchell R., "Unified high-temperature char combustion
+        kinetics for a suite of coals of various rank", 24th Symposium in
+        Combustion, The Combustion Institute, 1992, p 1243-1250
+
+    Model specifies the rate of char combustion.
+
+        C(s) + Sb*O2 -> CO2
+
+    where Sb is the stoichiometry of the reaction
+
+    Model validity:
+        Gas temperature: Tc > 1500 K
+        Particle sizes:  75 um -> 200 um
+        Pox > 0.3 atm        
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef COxidationHurtMitchell_H
+#define COxidationHurtMitchell_H
+
+#include "SurfaceReactionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward class declarations
+template<class CloudType>
+class COxidationHurtMitchell;
+
+/*---------------------------------------------------------------------------*\
+                  Class COxidationHurtMitchell Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class COxidationHurtMitchell
+:
+    public SurfaceReactionModel<CloudType>
+{
+    // Private data
+
+        // Model constants
+
+            //- Stoichiometry of reaction
+            const scalar Sb_;
+
+
+        // Addressing
+
+            //- Cs position in global/local lists
+            label CsLocalId_;
+
+            //- Ash position in global/local lists
+            label ashLocalId_;
+
+            //- O2 position in global list
+            label O2GlobalId_;
+
+            //- CO2 positions in global list
+            label CO2GlobalId_;
+
+
+        // Local copies of thermo properties
+
+            //- Molecular weight of C [kg/kmol]
+            scalar WC_;
+
+            //- Molecular weight of O2 [kg/kmol]
+            scalar WO2_;
+
+            //- Formation enthalpy for CO2 [J/kg]
+            scalar HcCO2_;
+
+            //- Heat of reaction [J/kg] (optional)
+            scalar heatOfReaction_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("COxidationHurtMitchell");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        COxidationHurtMitchell
+        (
+            const dictionary& dict,
+            CloudType& owner
+        );
+
+        //- Construct copy
+        COxidationHurtMitchell
+        (
+            const COxidationHurtMitchell<CloudType>& srm
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<SurfaceReactionModel<CloudType> > clone() const
+        {
+            return autoPtr<SurfaceReactionModel<CloudType> >
+            (
+                new COxidationHurtMitchell<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~COxidationHurtMitchell();
+
+
+    // Member Functions
+
+        //- Update surface reactions
+        virtual scalar calculate
+        (
+            const scalar dt,
+            const label cellI,
+            const scalar d,
+            const scalar T,
+            const scalar Tc,
+            const scalar pc,
+            const scalar rhoc,
+            const scalar mass,
+            const scalarField& YGas,
+            const scalarField& YLiquid,
+            const scalarField& YSolid,
+            const scalarField& YMixture,
+            const scalar N,
+            scalarField& dMassGas,
+            scalarField& dMassLiquid,
+            scalarField& dMassSolid,
+            scalarField& dMassSRCarrier
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "COxidationHurtMitchell.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index b43c3fc690..c350d21017 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -661,7 +661,7 @@ void Foam::KinematicCloud<CloudType>::evolve()
 
 template<class CloudType>
 template<class TrackData>
-void  Foam::KinematicCloud<CloudType>::motion(TrackData& td)
+void Foam::KinematicCloud<CloudType>::motion(TrackData& td)
 {
     td.part() = TrackData::tpLinearTrack;
     CloudType::move(td,  solution_.trackTime());
@@ -670,6 +670,152 @@ void  Foam::KinematicCloud<CloudType>::motion(TrackData& td)
 }
 
 
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::patchData
+(
+    const parcelType& p,
+    const polyPatch& pp,
+    const scalar trackFraction,
+    const tetIndices& tetIs,
+    vector& nw,
+    vector& Up
+) const
+{
+    label patchI = pp.index();
+    label patchFaceI = pp.whichFace(p.face());
+
+    vector n = tetIs.faceTri(mesh_).normal();
+    n /= mag(n);
+
+    vector U = U_.boundaryField()[patchI][patchFaceI];
+
+    // Unless the face is rotating, the required normal is n;
+    nw = n;
+
+    if (!mesh_.moving())
+    {
+        // Only wall patches may have a non-zero wall velocity from
+        // the velocity field when the mesh is not moving.
+
+        if (isA<wallPolyPatch>(pp))
+        {
+            Up = U;
+        }
+        else
+        {
+            Up = vector::zero;
+        }
+    }
+    else
+    {
+        vector U00 = U_.oldTime().boundaryField()[patchI][patchFaceI];
+
+        vector n00 = tetIs.oldFaceTri(mesh_).normal();
+
+        // Difference in normal over timestep
+        vector dn = vector::zero;
+
+        if (mag(n00) > SMALL)
+        {
+            // If the old normal is zero (for example in layer
+            // addition) then use the current normal, meaning that the
+            // motion can only be translational, and dn remains zero,
+            // otherwise, calculate dn:
+
+            n00 /= mag(n00);
+
+            dn = n - n00;
+        }
+
+        // Total fraction through the timestep of the motion,
+        // including stepFraction before the current tracking step
+        // and the current trackFraction
+        // i.e.
+        // let s = stepFraction, t = trackFraction
+        // Motion of x in time:
+        // |-----------------|---------|---------|
+        // x00               x0        xi        x
+        //
+        // where xi is the correct value of x at the required
+        // tracking instant.
+        //
+        // x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
+        //
+        // i.e. the motion covered by previous tracking portions
+        // within this timestep, and
+        //
+        // xi = x0 + t*(x - x0)
+        //    = t*x + (1 - t)*x0
+        //    = t*x + (1 - t)*(s*x + (1 - s)*x00)
+        //    = (s + t - s*t)*x + (1 - (s + t - s*t))*x00
+        //
+        // let m = (s + t - s*t)
+        //
+        // xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
+        //
+        // In the same form as before.
+
+        scalar m =
+            p.stepFraction()
+          + trackFraction
+          - (p.stepFraction()*trackFraction);
+
+        // When the mesh is moving, the velocity field on wall patches
+        // will contain the velocity associated with the motion of the
+        // mesh, in which case it is interpolated in time using m.
+        // For other patches the face velocity will need to be
+        // reconstructed from the face centre motion.
+
+        const vector& Cf = mesh_.faceCentres()[p.face()];
+
+        vector Cf00 = mesh_.faces()[p.face()].centre(mesh_.oldPoints());
+
+        if (isA<wallPolyPatch>(pp))
+        {
+            Up = U00 + m*(U - U00);
+        }
+        else
+        {
+            Up = (Cf - Cf00)/mesh_.time().deltaTValue();
+        }
+
+        if (mag(dn) > SMALL)
+        {
+            // Rotational motion, nw requires interpolation and a
+            // rotational velocity around face centre correction to Up
+            // is required.
+
+            nw = n00 + m*dn;
+
+            // Cf at tracking instant
+            vector Cfi = Cf00 + m*(Cf - Cf00);
+
+            // Normal vector cross product
+            vector omega = (n00 ^ n);
+
+            scalar magOmega = mag(omega);
+
+            // magOmega = sin(angle between unit normals)
+            // Normalise omega vector by magOmega, then multiply by
+            // angle/dt to give the correct angular velocity vector.
+            omega *= Foam::asin(magOmega)/(magOmega*mesh_.time().deltaTValue());
+
+            // Project position onto face and calculate this position
+            // relative to the face centre.
+            vector facePos =
+                p.position()
+              - ((p.position() - Cfi) & nw)*nw
+              - Cfi;
+
+            Up += (omega ^ facePos);
+        }
+
+        // No further action is required if the motion is
+        // translational only, nw and Up have already been set.
+    }
+}
+
+
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
 {
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index 4cd80923da..ed373d0f6b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -554,6 +554,17 @@ public:
             template<class TrackData>
             void motion(TrackData& td);
 
+            //- Calculate the patch normal and velocity to interact with,
+            //  accounting for patch motion if required.
+            void patchData
+            (
+                const parcelType& p,
+                const polyPatch& pp,
+                const scalar trackFraction,
+                const tetIndices& tetIs,
+                vector& normal,
+                vector& Up
+            ) const;
 
         // Mapping
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index 4a83808fa2..00fb7b5a96 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -87,8 +87,6 @@ void Foam::ReactingCloud<CloudType>::cloudReset(ReactingCloud<CloudType>& c)
 
     compositionModel_.reset(c.compositionModel_.ptr());
     phaseChangeModel_.reset(c.phaseChangeModel_.ptr());
-
-    dMassPhaseChange_ = c.dMassPhaseChange_;
 }
 
 
@@ -111,8 +109,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     constProps_(this->particleProperties(), this->solution().active()),
     compositionModel_(NULL),
     phaseChangeModel_(NULL),
-    rhoTrans_(thermo.carrier().species().size()),
-    dMassPhaseChange_(0.0)
+    rhoTrans_(thermo.carrier().species().size())
 {
     if (this->solution().active())
     {
@@ -167,8 +164,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     constProps_(c.constProps_),
     compositionModel_(c.compositionModel_->clone()),
     phaseChangeModel_(c.phaseChangeModel_->clone()),
-    rhoTrans_(c.rhoTrans_.size()),
-    dMassPhaseChange_(c.dMassPhaseChange_)
+    rhoTrans_(c.rhoTrans_.size())
 {
     forAll(c.rhoTrans_, i)
     {
@@ -209,8 +205,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     compositionModel_(c.compositionModel_->clone()),
 //    compositionModel_(NULL),
     phaseChangeModel_(NULL),
-    rhoTrans_(0),
-    dMassPhaseChange_(0.0)
+    rhoTrans_(0)
 {}
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
index a16f1255af..2797250cc3 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
@@ -121,12 +121,6 @@ protected:
             PtrList<DimensionedField<scalar, volMesh> > rhoTrans_;
 
 
-        // Check
-
-            //- Total mass transferred to continuous phase via phase change
-            scalar dMassPhaseChange_;
-
-
     // Protected Member Functions
 
         // New parcel helper functions
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index 9eaa877057..874936990a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -51,7 +51,69 @@ void Foam::ThermoCloud<CloudType>::setModels()
         ).ptr()
     );
 
-    this->subModelProperties().lookup("radiation") >> radiation_;
+    if (this->solution().coupled())
+    {
+        this->subModelProperties().lookup("radiation") >> radiation_;
+    }
+
+    if (radiation_)
+    {
+        radAreaP_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaP",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar("zero", dimArea, 0.0)
+            )
+        );
+
+        radT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar("zero", pow4(dimTemperature), 0.0)
+            )
+        );
+
+        radAreaPT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaPT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar
+                (
+                    "zero",
+                    sqr(dimLength)*pow4(dimTemperature),
+                    0.0
+                )
+            )
+        );
+    }
 }
 
 
@@ -98,13 +160,16 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
     heatTransferModel_(NULL),
     TIntegrator_(NULL),
     radiation_(false),
+    radAreaP_(NULL),
+    radT4_(NULL),
+    radAreaPT4_(NULL),
     hsTrans_
     (
         new DimensionedField<scalar, volMesh>
         (
             IOobject
             (
-                this->name() + "hsTrans",
+                this->name() + "::hsTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::READ_IF_PRESENT,
@@ -120,7 +185,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
         (
             IOobject
             (
-                this->name() + "hsCoeff",
+                this->name() + "::hsCoeff",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::READ_IF_PRESENT,
@@ -165,13 +230,16 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
     heatTransferModel_(c.heatTransferModel_->clone()),
     TIntegrator_(c.TIntegrator_->clone()),
     radiation_(c.radiation_),
+    radAreaP_(NULL),
+    radT4_(NULL),
+    radAreaPT4_(NULL),
     hsTrans_
     (
         new DimensionedField<scalar, volMesh>
         (
             IOobject
             (
-                this->name() + "hsTrans",
+                this->name() + "::hsTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -187,7 +255,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
         (
             IOobject
             (
-                this->name() + "hsCoeff",
+                this->name() + "::hsCoeff",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -197,7 +265,61 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
             c.hsCoeff()
         )
     )
-{}
+{
+    if (radiation_)
+    {
+        radAreaP_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaP",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                c.radAreaP()
+            )
+        );
+
+        radT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                c.radT4()
+            )
+        );
+
+        radAreaPT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaPT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                c.radAreaPT4()
+            )
+        );
+    }
+}
 
 
 template<class CloudType>
@@ -218,6 +340,9 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
     heatTransferModel_(NULL),
     TIntegrator_(NULL),
     radiation_(false),
+    radAreaP_(NULL),
+    radT4_(NULL),
+    radAreaPT4_(NULL),
     hsTrans_(NULL),
     hsCoeff_(NULL)
 {}
@@ -285,6 +410,13 @@ void Foam::ThermoCloud<CloudType>::resetSourceTerms()
     CloudType::resetSourceTerms();
     hsTrans_->field() = 0.0;
     hsCoeff_->field() = 0.0;
+
+    if (radiation_)
+    {
+        radAreaP_->field() = 0.0;
+        radT4_->field() = 0.0;
+        radAreaPT4_->field() = 0.0;
+    }
 }
 
 
@@ -298,6 +430,13 @@ void Foam::ThermoCloud<CloudType>::relaxSources
 
     this->relax(hsTrans_(), cloudOldTime.hsTrans(), "h");
     this->relax(hsCoeff_(), cloudOldTime.hsCoeff(), "h");
+
+    if (radiation_)
+    {
+        this->relax(radAreaP_(), cloudOldTime.radAreaP(), "radiation");
+        this->relax(radT4_(), cloudOldTime.radT4(), "radiation");
+        this->relax(radAreaPT4_(), cloudOldTime.radAreaPT4(), "radiation");
+    }
 }
 
 
@@ -308,6 +447,13 @@ void Foam::ThermoCloud<CloudType>::scaleSources()
 
     this->scale(hsTrans_(), "h");
     this->scale(hsCoeff_(), "h");
+
+    if (radiation_)
+    {
+        this->scale(radAreaP_(), "radiation");
+        this->scale(radT4_(), "radiation");
+        this->scale(radAreaPT4_(), "radiation");
+    }
 }
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
index 4acd251621..88cfde3172 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
@@ -132,6 +132,15 @@ protected:
             //- Include radiation
             Switch radiation_;
 
+            //- Radiation sum of parcel projected areas
+            autoPtr<DimensionedField<scalar, volMesh> > radAreaP_;
+
+            //- Radiation sum of parcel temperature^4
+            autoPtr<DimensionedField<scalar, volMesh> > radT4_;
+
+            //- Radiation sum of parcel projected areas * temperature^4
+            autoPtr<DimensionedField<scalar, volMesh> > radAreaPT4_;
+
 
         // Sources
 
@@ -244,6 +253,26 @@ public:
                  //- Radiation flag
                 inline bool radiation() const;
 
+                //- Radiation sum of parcel projected areas [m2]
+                inline DimensionedField<scalar, volMesh>& radAreaP();
+
+                //- Radiation sum of parcel projected areas [m2]
+                inline const DimensionedField<scalar, volMesh>&
+                    radAreaP() const;
+
+                //- Radiation sum of parcel temperature^4 [K4]
+                inline DimensionedField<scalar, volMesh>& radT4();
+
+                //- Radiation sum of parcel temperature^4 [K4]
+                inline const DimensionedField<scalar, volMesh>& radT4() const;
+
+                //- Radiation sum of parcel projected area*temperature^4 [m2K4]
+                inline DimensionedField<scalar, volMesh>& radAreaPT4();
+
+                //- Radiation sum of parcel temperature^4 [m2K4]
+                inline const DimensionedField<scalar, volMesh>&
+                    radAreaPT4() const;
+
 
             // Sources
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
index bac95e49d5..764595fece 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
@@ -89,6 +89,114 @@ inline bool Foam::ThermoCloud<CloudType>::radiation() const
 }
 
 
+template<class CloudType>
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaP()
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaP()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaP_();
+}
+
+
+template<class CloudType>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaP() const
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaP()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaP_();
+}
+
+
+template<class CloudType>
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radT4()
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radT4_();
+}
+
+
+template<class CloudType>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radT4() const
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radT4_();
+}
+
+
+template<class CloudType>
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaPT4()
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaPT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaPT4_();
+}
+
+
+template<class CloudType>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaPT4() const
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaPT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaPT4_();
+}
+
+
 template<class CloudType>
 inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
 Foam::ThermoCloud<CloudType>::hsTrans()
@@ -182,21 +290,15 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::Ep() const
         )
     );
 
-    // Need to check if coupled as field is created on-the-fly
-    if (radiation_ && this->solution().coupled())
+    if (radiation_)
     {
         scalarField& Ep = tEp().internalField();
+        const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
+        const scalarField& sumAreaPT4 = radAreaPT4_->field();
 
-        forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
-        {
-            const parcelType& p = iter();
-            const label cellI = p.cell();
-            Ep[cellI] += p.nParticle()*p.areaP()*pow4(p.T());
-        }
-
-        Ep *= epsilon*physicoChemical::sigma.value()/V;
+        Ep = sumAreaPT4*epsilon*physicoChemical::sigma.value()/V/dt;
     }
 
     return tEp;
@@ -224,21 +326,15 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::ap() const
         )
     );
 
-    // Need to check if coupled as field is created on-the-fly
-    if (radiation_ && this->solution().coupled())
+    if (radiation_)
     {
         scalarField& ap = tap().internalField();
+        const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
+        const scalarField& sumAreaP = radAreaP_->field();
 
-        forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
-        {
-            const parcelType& p = iter();
-            const label cellI = p.cell();
-            ap[cellI] += p.nParticle()*p.areaP();
-        }
-
-        ap *= epsilon/V;
+        ap = sumAreaP*epsilon/V/dt;
     }
 
     return tap;
@@ -267,23 +363,16 @@ Foam::ThermoCloud<CloudType>::sigmap() const
         )
     );
 
-    // Need to check if coupled as field is created on-the-fly
-    if (radiation_ && this->solution().coupled())
+    if (radiation_)
     {
         scalarField& sigmap = tsigmap().internalField();
-
+        const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
         const scalar f = constProps_.f0();
+        const scalarField& sumAreaP = radAreaP_->field();
 
-        forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
-        {
-            const parcelType& p = iter();
-            const label cellI = p.cell();
-            sigmap[cellI] += p.nParticle()*p.areaP();
-        }
-
-        sigmap *= (1.0 - f)*(1.0 - epsilon)/V;
+        sigmap *= sumAreaP*(1.0 - f)*(1.0 - epsilon)/V/dt;
     }
 
     return tsigmap;
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 9ae6b2f5b9..8fa42ac0ec 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -363,7 +363,7 @@ bool Foam::KinematicParcel<ParcelType>::hitPatch
         static_cast<typename TrackData::cloudType::parcelType&>(*this);
 
     // Invoke post-processing model
-    td.cloud().functions().postPatch(p, patchI, pp.whichFace(p.face()));
+    td.cloud().functions().postPatch(p, pp, trackFraction, tetIs);
 
     // Invoke surface film model
     if (td.cloud().surfaceFilm().transferParcel(p, pp, td.keepParticle))
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 4be1fcdd12..84942b0349 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -268,12 +268,15 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     (
         td,
         dt,
+        this->age_,
         Ts,
         d0,
         T0,
         mass0,
         this->mass0_,
         YMix[GAS]*YGas_,
+        YMix[LIQ]*YLiquid_,
+        YMix[SLD]*YSolid_,
         canCombust_,
         dMassDV,
         Sh,
@@ -466,6 +469,16 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         // Update sensible enthalpy transfer
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
         td.cloud().hsCoeff()[cellI] += np0*Sph;
+
+        // Update radiation fields
+        if (td.cloud().radiation())
+        {
+            const scalar ap = this->areaP();
+            const scalar T4 = pow4(this->T_);
+            td.cloud().radAreaP()[cellI] += dt*np0*ap;
+            td.cloud().radT4()[cellI] += dt*np0*T4;
+            td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
+        }
     }
 }
 
@@ -476,12 +489,15 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 (
     TrackData& td,
     const scalar dt,
+    const scalar age,
     const scalar Ts,
     const scalar d,
     const scalar T,
     const scalar mass,
     const scalar mass0,
     const scalarField& YGasEff,
+    const scalarField& YLiquidEff,
+    const scalarField& YSolidEff,
     bool& canCombust,
     scalarField& dMassDV,
     scalar& Sh,
@@ -510,10 +526,13 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     td.cloud().devolatilisation().calculate
     (
         dt,
+        age,
         mass0,
         mass,
         T,
         YGasEff,
+        YLiquidEff,
+        YSolidEff,
         canCombust,
         dMassDV
     );
@@ -638,7 +657,8 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
     ParcelType(p),
     YGas_(p.YGas_),
     YLiquid_(p.YLiquid_),
-    YSolid_(p.YSolid_)
+    YSolid_(p.YSolid_),
+    canCombust_(p.canCombust_)
 {}
 
 
@@ -652,7 +672,8 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
     ParcelType(p, mesh),
     YGas_(p.YGas_),
     YLiquid_(p.YLiquid_),
-    YSolid_(p.YSolid_)
+    YSolid_(p.YSolid_),
+    canCombust_(p.canCombust_)
 {}
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index 6707b2a52d..db20b8d9c7 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -196,12 +196,15 @@ protected:
         (
             TrackData& td,
             const scalar dt,           // timestep
-            const scalar Ts,           // Surface temperature
+            const scalar age,          // age
+            const scalar Ts,           // surface temperature
             const scalar d,            // diameter
             const scalar T,            // temperature
             const scalar mass,         // mass
             const scalar mass0,        // mass (initial on injection)
-            const scalarField& YGasEff,// Gas component mass fractions
+            const scalarField& YGasEff,// gas component mass fractions
+            const scalarField& YLiquidEff,// liquid component mass fractions
+            const scalarField& YSolidEff,// solid component mass fractions
             bool& canCombust,          // 'can combust' flag
             scalarField& dMassDV,      // mass transfer - local to particle
             scalar& Sh,                // explicit particle enthalpy source
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index cfd21ceed0..e9bc815304 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -458,6 +458,16 @@ void Foam::ReactingParcel<ParcelType>::calc
         // Update sensible enthalpy transfer
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
         td.cloud().hsCoeff()[cellI] += np0*Sph;
+
+        // Update radiation fields
+        if (td.cloud().radiation())
+        {
+            const scalar ap = this->areaP();
+            const scalar T4 = pow4(this->T_);
+            td.cloud().radAreaP()[cellI] += dt*np0*ap;
+            td.cloud().radT4()[cellI] += dt*np0*T4;
+            td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
+        }
     }
 }
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index b64b1e8f99..cb9b714204 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -144,7 +144,7 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
         Ts = td.cloud().constProps().TMin();
     }
 
-    // Assuming thermo props vary linearly with T for small dT
+    // Assuming thermo props vary linearly with T for small d(T)
     const scalar TRatio = Tc_/Ts;
 
     rhos = this->rhoc_*TRatio;
@@ -252,6 +252,16 @@ void Foam::ThermoParcel<ParcelType>::calc
 
         // Update sensible enthalpy coefficient
         td.cloud().hsCoeff()[cellI] += np0*Sph;
+
+        // Update radiation fields
+        if (td.cloud().radiation())
+        {
+            const scalar ap = this->areaP();
+            const scalar T4 = pow4(this->T_);
+            td.cloud().radAreaP()[cellI] += dt*np0*ap;
+            td.cloud().radT4()[cellI] += dt*np0*T4;
+            td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
+        }
     }
 }
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C
index b8409231ff..dbc359ae9a 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,8 +107,9 @@ template<class CloudType>
 void Foam::CloudFunctionObject<CloudType>::postPatch
 (
     const typename CloudType::parcelType&,
-    const label,
-    const label
+    const polyPatch&,
+    const scalar,
+    const tetIndices&
 )
 {
     // do nothing
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H
index b0818aa708..5baf0f7ae9 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -141,8 +141,9 @@ public:
             virtual void postPatch
             (
                 const typename CloudType::parcelType& p,
-                const label patchI,
-                const label patchFaceI
+                const polyPatch& pp,
+                const scalar trackFraction,
+                const tetIndices& testIs
             );
 
             //- Post-face hook
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C
index 556c72908c..2f97468e36 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -146,13 +146,14 @@ template<class CloudType>
 void Foam::CloudFunctionObjectList<CloudType>::postPatch
 (
     const typename CloudType::parcelType& p,
-    const label patchI,
-    const label patchFaceI
+    const polyPatch& pp,
+    const scalar trackFraction,
+    const tetIndices& tetIs
 )
 {
     forAll(*this, i)
     {
-        this->operator[](i).postPatch(p, patchI, patchFaceI);
+        this->operator[](i).postPatch(p, pp, trackFraction, tetIs);
     }
 }
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.H
index b8b296fe21..5daeab9f29 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -121,8 +121,9 @@ public:
             virtual void postPatch
             (
                 const typename CloudType::parcelType& p,
-                const label patchI,
-                const label patchFaceI
+                const polyPatch& pp,
+                const scalar trackFraction,
+                const tetIndices& tetIs
             );
 
             //- Post-face hook
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
index 87d03eade7..81772debfe 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
@@ -166,29 +166,32 @@ template<class CloudType>
 void Foam::ParticleErosion<CloudType>::postPatch
 (
     const parcelType& p,
-    const label patchI,
-    const label patchFaceI
+    const polyPatch& pp,
+    const scalar trackFraction,
+    const tetIndices& tetIs
 )
 {
+    const label patchI = pp.index();
+
     const label localPatchI = applyToPatch(patchI);
 
     if (localPatchI != -1)
     {
-        const fvMesh& mesh = this->owner().mesh();
+        vector nw;
+        vector Up;
 
         // patch-normal direction
-        vector nw = p.currentTetIndices().faceTri(mesh).normal();
+        this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
 
-        // particle direction of travel
-        const vector& U = p.U();
+        // particle velocity reletive to patch
+        const vector& U = p.U() - Up;
 
         // quick reject if particle travelling away from the patch
-        if ((-nw & U) < 0)
+        if ((nw & U) < 0)
         {
             return;
         }
 
-        nw /= mag(nw);
         const scalar magU = mag(U);
         const vector Udir = U/magU;
 
@@ -197,6 +200,7 @@ void Foam::ParticleErosion<CloudType>::postPatch
 
         const scalar coeff = p.nParticle()*p.mass()*sqr(magU)/(p_*psi_*K_);
 
+        const label patchFaceI = pp.whichFace(p.face());
         scalar& Q = QPtr_->boundaryField()[patchI][patchFaceI];
         if (tan(alpha) < K_/6.0)
         {
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.H
index c23dfc4c40..c741bd4004 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -126,8 +126,9 @@ public:
             virtual void postPatch
             (
                 const parcelType& p,
-                const label patchI,
-                const label patchFaceI
+                const polyPatch& pp,
+                const scalar trackFraction,
+                const tetIndices& tetIs
             );
 };
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.C
index afbd9bfbd9..f66a872cfc 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -218,11 +218,14 @@ template<class CloudType>
 void Foam::PatchPostProcessing<CloudType>::postPatch
 (
     const parcelType& p,
-    const label patchI,
-    const label
+    const polyPatch& pp,
+    const scalar,
+    const tetIndices& tetIs
 )
 {
+    const label patchI = pp.index();
     const label localPatchI = applyToPatch(patchI);
+
     if (localPatchI != -1 && patchData_[localPatchI].size() < maxStoredParcels_)
     {
         times_[localPatchI].append(this->owner().time().value());
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.H
index 9945d24b58..f7d618e6ce 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,8 +127,9 @@ public:
             virtual void postPatch
             (
                 const parcelType& p,
-                const label patchI,
-                const label patchFaceI
+                const polyPatch& pp,
+                const scalar trackFraction,
+                const tetIndices& tetIs
             );
 };
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
index 5319b75b27..6315eac827 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
@@ -306,8 +306,9 @@ void Foam::PairCollision<CloudType>::wallInteraction()
                         this->owner().functions().postPatch
                         (
                             p,
-                            patchI,
-                            patchFaceI
+                            mesh.boundaryMesh()[patchI],
+                            1.0,
+                            p.currentTetIndices()
                         );
                      }
                 }
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index f1c6d2e524..2fa553a679 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -284,6 +284,7 @@ Foam::InjectionModel<CloudType>::InjectionModel(CloudType& owner)
     SOI_(0.0),
     volumeTotal_(0.0),
     massTotal_(0.0),
+    massFlowRate_(owner.db().time(), "massFlowRate"),
     massInjected_(this->template getModelProperty<scalar>("massInjected")),
     nInjections_(this->template getModelProperty<label>("nInjections")),
     parcelsAddedTotal_
@@ -310,6 +311,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
     SOI_(0.0),
     volumeTotal_(0.0),
     massTotal_(0.0),
+    massFlowRate_(owner.db().time(), "massFlowRate"),
     massInjected_(this->template getModelProperty<scalar>("massInjected")),
     nInjections_(this->template getModelProperty<scalar>("nInjections")),
     parcelsAddedTotal_
@@ -335,7 +337,8 @@ Foam::InjectionModel<CloudType>::InjectionModel
     }
     else
     {
-        this->coeffDict().lookup("massFlowRate") >> massTotal_;
+        massFlowRate_.reset(this->coeffDict());
+        massTotal_ = massFlowRate_.value(owner.db().time().value());
     }
 
     const word parcelBasisType = this->coeffDict().lookup("parcelBasisType");
@@ -384,6 +387,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
     SOI_(im.SOI_),
     volumeTotal_(im.volumeTotal_),
     massTotal_(im.massTotal_),
+    massFlowRate_(im.massFlowRate_),
     massInjected_(im.massInjected_),
     nInjections_(im.nInjections_),
     parcelsAddedTotal_(im.parcelsAddedTotal_),
@@ -458,7 +462,16 @@ Foam::scalar Foam::InjectionModel<CloudType>::volumeToInject
 template<class CloudType>
 Foam::scalar Foam::InjectionModel<CloudType>::averageParcelMass()
 {
-    label nTotal = parcelsToInject(0.0, timeEnd() - timeStart());
+    label nTotal = 0.0;
+    if (this->owner().solution().transient())
+    {
+        nTotal = parcelsToInject(0.0, timeEnd() - timeStart());
+    }
+    else
+    {
+        nTotal = parcelsToInject(0.0, 1.0);
+    }
+
     return massTotal_/nTotal;
 }
 
@@ -594,6 +607,8 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
     const polyMesh& mesh = this->owner().mesh();
     typename TrackData::cloudType& cloud = td.cloud();
 
+    massTotal_ = massFlowRate_.value(mesh.time().value());
+
     // Reset counters
     time0_ = 0.0;
     label parcelsAdded = 0;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
index a3d67947b2..c3e0534ff4 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
@@ -53,6 +53,7 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "SubModelBase.H"
 #include "vector.H"
+#include "TimeDataEntry.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -102,6 +103,9 @@ protected:
             //- Total mass to inject [kg]
             scalar massTotal_;
 
+            //- Mass flow rate profile for steady calculations
+            TimeDataEntry<scalar> massFlowRate_;
+
             //- Total mass injected to date [kg]
             scalar massInjected_;
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C
index 2350ee62f8..8ffae4f75f 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C
@@ -233,7 +233,7 @@ bool Foam::LocalInteraction<CloudType>::correct
                 vector nw;
                 vector Up;
 
-                this->patchData(p, pp, trackFraction, tetIs, nw, Up);
+                this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
 
                 // Calculate motion relative to patch velocity
                 U -= Up;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C
index 9d2566ab36..6ffabc3715 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -179,159 +179,6 @@ bool Foam::PatchInteractionModel<CloudType>::correct
 }
 
 
-template<class CloudType>
-void Foam::PatchInteractionModel<CloudType>::patchData
-(
-    typename CloudType::parcelType& p,
-    const polyPatch& pp,
-    const scalar trackFraction,
-    const tetIndices& tetIs,
-    vector& nw,
-    vector& Up
-) const
-{
-    const fvMesh& mesh = this->owner().mesh();
-
-    const volVectorField& Ufield =
-        mesh.objectRegistry::lookupObject<volVectorField>(UName_);
-
-    label patchI = pp.index();
-    label patchFaceI = pp.whichFace(p.face());
-
-    vector n = tetIs.faceTri(mesh).normal();
-    n /= mag(n);
-
-    vector U = Ufield.boundaryField()[patchI][patchFaceI];
-
-    // Unless the face is rotating, the required normal is n;
-    nw = n;
-
-    if (!mesh.moving())
-    {
-        // Only wall patches may have a non-zero wall velocity from
-        // the velocity field when the mesh is not moving.
-
-        if (isA<wallPolyPatch>(pp))
-        {
-            Up = U;
-        }
-        else
-        {
-            Up = vector::zero;
-        }
-    }
-    else
-    {
-        vector U00 = Ufield.oldTime().boundaryField()[patchI][patchFaceI];
-
-        vector n00 = tetIs.oldFaceTri(mesh).normal();
-
-        // Difference in normal over timestep
-        vector dn = vector::zero;
-
-        if (mag(n00) > SMALL)
-        {
-            // If the old normal is zero (for example in layer
-            // addition) then use the current normal, meaning that the
-            // motion can only be translational, and dn remains zero,
-            // otherwise, calculate dn:
-
-            n00 /= mag(n00);
-
-            dn = n - n00;
-        }
-
-        // Total fraction thought the timestep of the motion,
-        // including stepFraction before the current tracking step
-        // and the current trackFraction
-        // i.e.
-        // let s = stepFraction, t = trackFraction
-        // Motion of x in time:
-        // |-----------------|---------|---------|
-        // x00               x0        xi        x
-        //
-        // where xi is the correct value of x at the required
-        // tracking instant.
-        //
-        // x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
-        //
-        // i.e. the motion covered by previous tracking portions
-        // within this timestep, and
-        //
-        // xi = x0 + t*(x - x0)
-        //    = t*x + (1 - t)*x0
-        //    = t*x + (1 - t)*(s*x + (1 - s)*x00)
-        //    = (s + t - s*t)*x + (1 - (s + t - s*t))*x00
-        //
-        // let m = (s + t - s*t)
-        //
-        // xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
-        //
-        // In the same form as before.
-
-        scalar m =
-            p.stepFraction()
-          + trackFraction
-          - (p.stepFraction()*trackFraction);
-
-        // When the mesh is moving, the velocity field on wall patches
-        // will contain the velocity associated with the motion of the
-        // mesh, in which case it is interpolated in time using m.
-        // For other patches the face velocity will need to be
-        // reconstructed from the face centre motion.
-
-        const vector& Cf = mesh.faceCentres()[p.face()];
-
-        vector Cf00 = mesh.faces()[p.face()].centre(mesh.oldPoints());
-
-        if (isA<wallPolyPatch>(pp))
-        {
-            Up = U00 + m*(U - U00);
-        }
-        else
-        {
-            Up = (Cf - Cf00)/this->owner().time().deltaTValue();
-        }
-
-        if (mag(dn) > SMALL)
-        {
-            // Rotational motion, nw requires interpolation and a
-            // rotational velocity around face centre correction to Up
-            // is required.
-
-            nw = n00 + m*dn;
-
-            // Cf at tracking instant
-            vector Cfi = Cf00 + m*(Cf - Cf00);
-
-            // Normal vector cross product
-            vector omega = (n00 ^ n);
-
-            scalar magOmega = mag(omega);
-
-            // magOmega = sin(angle between unit normals)
-            // Normalise omega vector by magOmega, then multiply by
-            // angle/dt to give the correct angular velocity vector.
-            omega *=
-                Foam::asin(magOmega)
-               /(magOmega*this->owner().time().deltaTValue());
-
-            // Project position onto face and calculate this position
-            // relative to the face centre.
-            vector facePos =
-                p.position()
-              - ((p.position() - Cfi) & nw)*nw
-              - Cfi;
-
-            Up += (omega ^ facePos);
-        }
-
-        // No further action is required if the motion is
-        // translational only, nw and Up have already been set.
-    }
-}
-
-
 template<class CloudType>
 void Foam::PatchInteractionModel<CloudType>::info(Ostream& os)
 {
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.H b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.H
index 293d2cdb45..7f3ea38198 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -164,18 +164,6 @@ public:
             const tetIndices& tetIs
         );
 
-        //- Calculate the patch normal and velocity to interact with,
-        //  accounting for patch motion if required.
-        void patchData
-        (
-            typename CloudType::parcelType& p,
-            const polyPatch& pp,
-            const scalar trackFraction,
-            const tetIndices& tetIs,
-            vector& normal,
-            vector& Up
-        ) const;
-
 
         // I-O
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C
index 147408204b..f6c1fda460 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,7 +74,7 @@ bool Foam::Rebound<CloudType>::correct
     vector nw;
     vector Up;
 
-    this->patchData(p, pp, trackFraction, tetIs, nw, Up);
+    this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
 
     // Calculate motion relative to patch velocity
     U -= Up;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
index b3052a4ec2..d9b38daccf 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
@@ -148,7 +148,7 @@ bool Foam::StandardWallInteraction<CloudType>::correct
                 vector nw;
                 vector Up;
 
-                this->patchData(p, pp, trackFraction, tetIs, nw, Up);
+                this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
 
                 // Calculate motion relative to patch velocity
                 U -= Up;
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
index 76a22ffe44..72e56bf884 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
@@ -156,18 +156,20 @@ Foam::CompositionModel<CloudType>::componentNames(const label phaseI) const
 template<class CloudType>
 Foam::label Foam::CompositionModel<CloudType>::globalCarrierId
 (
-    const word& cmptName
+    const word& cmptName,
+    const bool allowNotFound
 ) const
 {
     label id = thermo_.carrierId(cmptName);
 
-    if (id < 0)
+    if (id < 0 && !allowNotFound)
     {
         FatalErrorIn
         (
             "Foam::label Foam::CompositionModel<CloudType>::globalCarrierId"
             "("
-                "const word&"
+                "const word&, "
+                "const bool"
             ") const"
         )   << "Unable to determine global id for requested component "
             << cmptName << ". Available components are " << nl
@@ -182,19 +184,21 @@ template<class CloudType>
 Foam::label Foam::CompositionModel<CloudType>::globalId
 (
     const label phaseI,
-    const word& cmptName
+    const word& cmptName,
+    const bool allowNotFound
 ) const
 {
     label id = phaseProps_[phaseI].globalId(cmptName);
 
-    if (id < 0)
+    if (id < 0 && !allowNotFound)
     {
         FatalErrorIn
         (
             "Foam::label Foam::CompositionModel<CloudType>::globalId"
             "("
                 "const label, "
-                "const word&"
+                "const word&, "
+                "const bool"
             ") const"
         )   << "Unable to determine global id for requested component "
             << cmptName << abort(FatalError);
@@ -218,19 +222,21 @@ template<class CloudType>
 Foam::label Foam::CompositionModel<CloudType>::localId
 (
     const label phaseI,
-    const word& cmptName
+    const word& cmptName,
+    const bool allowNotFound
 ) const
 {
     label id = phaseProps_[phaseI].id(cmptName);
 
-    if (id < 0)
+    if (id < 0 && !allowNotFound)
     {
         FatalErrorIn
         (
             "Foam::label Foam::CompositionModel<CloudType>::localId"
             "("
                 "const label, "
-                "const word&"
+                "const word&, "
+                "const bool"
             ") const"
         )   << "Unable to determine local id for component " << cmptName
             << abort(FatalError);
@@ -244,12 +250,13 @@ template<class CloudType>
 Foam::label Foam::CompositionModel<CloudType>::localToGlobalCarrierId
 (
     const label phaseI,
-    const label id
+    const label id,
+    const bool allowNotFound
 ) const
 {
     label gid = phaseProps_[phaseI].globalCarrierIds()[id];
 
-    if (gid < 0)
+    if (gid < 0 && !allowNotFound)
     {
         FatalErrorIn
         (
@@ -257,7 +264,8 @@ Foam::label Foam::CompositionModel<CloudType>::localToGlobalCarrierId
             "Foam::CompositionModel<CloudType>::localToGlobalCarrierId"
             "("
                 "const label, "
-                "const label"
+                "const label, "
+                "const bool"
             ") const"
         )   << "Unable to determine global carrier id for phase "
             << phaseI << " with local id " << id
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
index 724354902c..cabdd91eac 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -166,22 +166,37 @@ public:
                 const wordList& componentNames(const label phaseI) const;
 
                 //- Return global id of component cmptName in carrier thermo
-                label globalCarrierId(const word& cmptName) const;
+                label globalCarrierId
+                (
+                    const word& cmptName,
+                    const bool allowNotFound = false
+                ) const;
 
                 //- Return global id of component cmptName in phase phaseI
-                label globalId(const label phaseI, const word& cmptName) const;
+                label globalId
+                (
+                    const label phaseI,
+                    const word& cmptName,
+                    const bool allowNotFound = false
+                ) const;
 
                 //- Return global ids of for phase phaseI
                 const labelList& globalIds(const label phaseI) const;
 
                 //- Return local id of component cmptName in phase phaseI
-                label localId(const label phaseI, const word& cmptName) const;
+                label localId
+                (
+                    const label phaseI,
+                    const word& cmptName,
+                    const bool allowNotFound = false
+                ) const;
 
                 //- Return global carrier id of component given local id
                 label localToGlobalCarrierId
                 (
                     const label phaseI,
-                    const label id
+                    const label id,
+                    const bool allowNotFound = false
                 ) const;
 
                 //- Return the list of phase phaseI mass fractions
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C
index eb4cfe6efe..0de0772381 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,10 +102,13 @@ template<class CloudType>
 void Foam::ConstantRateDevolatilisation<CloudType>::calculate
 (
     const scalar dt,
+    const scalar age,
     const scalar mass0,
     const scalar mass,
     const scalar T,
     const scalarField& YGasEff,
+    const scalarField& YLiquidEff,
+    const scalarField& YSolidEff,
     bool& canCombust,
     scalarField& dMassDV
 ) const
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H
index d84810e81e..7faffb990d 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,10 +104,13 @@ public:
         virtual void calculate
         (
             const scalar dt,
+            const scalar age,
             const scalar mass0,
             const scalar mass,
             const scalar T,
             const scalarField& YGasEff,
+            const scalarField& YLiquidEff,
+            const scalarField& YSolidEff,
             bool& canCombust,
             scalarField& dMassDV
         ) const;
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.C
index 331ae8d78e..6b4325db4f 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.C
@@ -78,6 +78,9 @@ void Foam::DevolatilisationModel<CloudType>::calculate
     const scalar,
     const scalar,
     const scalar,
+    const scalar,
+    const scalarField&,
+    const scalarField&,
     const scalarField&,
     bool&,
     scalarField&
@@ -91,6 +94,9 @@ void Foam::DevolatilisationModel<CloudType>::calculate
             "const scalar, "
             "const scalar, "
             "const scalar, "
+            "const scalar, "
+            "const scalarField&, "
+            "const scalarField&, "
             "const scalarField&, "
             "bool&, "
             "scalarField&"
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.H
index 5b72c0a36d..7cff8d426a 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                      Class DevolatilisationModel Declaration
+                    Class DevolatilisationModel Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class CloudType>
@@ -126,10 +126,13 @@ public:
         virtual void calculate
         (
             const scalar dt,
+            const scalar age,
             const scalar mass0,
             const scalar mass,
             const scalar T,
             const scalarField& YGasEff,
+            const scalarField& YLiquidEff,
+            const scalarField& YSolidEff,
             bool& canCombust,
             scalarField& dMassDV
         ) const;
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.C
index f2b5819715..238cb4a02e 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,6 +71,9 @@ void Foam::NoDevolatilisation<CloudType>::calculate
     const scalar,
     const scalar,
     const scalar,
+    const scalar,
+    const scalarField&,
+    const scalarField&,
     const scalarField&,
     bool& canCombust,
     scalarField&
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.H
index 1b905d42fd..b3a617ba86 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,10 +85,13 @@ public:
         virtual void calculate
         (
             const scalar dt,
+            const scalar age,
             const scalar mass0,
             const scalar mass,
             const scalar T,
             const scalarField& YGasEff,
+            const scalarField& YLiquidEff,
+            const scalarField& YSolidEff,
             bool& canCombust,
             scalarField& dMassDV
         ) const;
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.C
index fd0f1429e3..fe4b856cfb 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,10 +105,13 @@ template<class CloudType>
 void Foam::SingleKineticRateDevolatilisation<CloudType>::calculate
 (
     const scalar dt,
+    const scalar age,
     const scalar mass0,
     const scalar mass,
     const scalar T,
     const scalarField& YGasEff,
+    const scalarField& YLiquidEff,
+    const scalarField& YSolidEff,
     bool& canCombust,
     scalarField& dMassDV
 ) const
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.H
index 535c7dfe1d..d9c3064ce9 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -229,10 +229,13 @@ public:
         virtual void calculate
         (
             const scalar dt,
+            const scalar age,
             const scalar mass0,
             const scalar mass,
             const scalar T,
             const scalarField& YGasEff,
+            const scalarField& YLiquidEff,
+            const scalarField& YSolidEff,
             bool& canCombust,
             scalarField& dMassDV
         ) const;
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index 0b0e3800f7..7e85e6723e 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
@@ -25,7 +25,6 @@ License
 
 #include "AMIInterpolation.H"
 #include "meshTools.H"
-#include "mergePoints.H"
 #include "mapDistribute.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -204,455 +203,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::resetTree
 }
 
 
-template<class SourcePatch, class TargetPatch>
-Foam::label Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcDistribution
-(
-    const SourcePatch& srcPatch,
-    const TargetPatch& tgtPatch
-) const
-{
-    label procI = 0;
-
-    if (Pstream::parRun())
-    {
-        List<label> facesPresentOnProc(Pstream::nProcs(), 0);
-        if ((srcPatch.size() > 0) || (tgtPatch.size() > 0))
-        {
-            facesPresentOnProc[Pstream::myProcNo()] = 1;
-        }
-        else
-        {
-            facesPresentOnProc[Pstream::myProcNo()] = 0;
-        }
-
-        Pstream::gatherList(facesPresentOnProc);
-        Pstream::scatterList(facesPresentOnProc);
-
-        label nHaveFaces = sum(facesPresentOnProc);
-
-        if (nHaveFaces > 1)
-        {
-            procI = -1;
-            if (debug)
-            {
-                Info<< "AMIInterpolation::calcDistribution: "
-                    << "AMI split across multiple processors" << endl;
-            }
-        }
-        else if (nHaveFaces == 1)
-        {
-            procI = findIndex(facesPresentOnProc, 1);
-            if (debug)
-            {
-                Info<< "AMIInterpolation::calcDistribution: "
-                    << "AMI local to processor" << procI << endl;
-            }
-        }
-    }
-
-
-    // Either not parallel or no faces on any processor
-    return procI;
-}
-
-
-template<class SourcePatch, class TargetPatch>
-Foam::label
-Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcOverlappingProcs
-(
-    const List<treeBoundBoxList>& procBb,
-    const treeBoundBox& bb,
-    boolList& overlaps
-) const
-{
-    overlaps.setSize(procBb.size());
-    overlaps = false;
-
-    label nOverlaps = 0;
-
-    forAll(procBb, procI)
-    {
-        const List<treeBoundBox>& bbs = procBb[procI];
-
-        forAll(bbs, bbI)
-        {
-            if (bbs[bbI].overlaps(bb))
-            {
-                overlaps[procI] = true;
-                nOverlaps++;
-                break;
-            }
-        }
-    }
-    return nOverlaps;
-}
-
-
-template<class SourcePatch, class TargetPatch>
-void Foam::AMIInterpolation<SourcePatch, TargetPatch>::distributePatches
-(
-    const mapDistribute& map,
-    const TargetPatch& pp,
-    const globalIndex& gi,
-    List<faceList>& faces,
-    List<pointField>& points,
-    List<labelList>& faceIDs
-) const
-{
-    PstreamBuffers pBufs(Pstream::nonBlocking);
-
-    for (label domain = 0; domain < Pstream::nProcs(); domain++)
-    {
-        const labelList& sendElems = map.subMap()[domain];
-
-        if (domain != Pstream::myProcNo() && sendElems.size())
-        {
-            labelList globalElems(sendElems.size());
-            forAll(sendElems, i)
-            {
-                globalElems[i] = gi.toGlobal(sendElems[i]);
-            }
-
-            faceList subFaces(UIndirectList<face>(pp, sendElems));
-            primitivePatch subPatch
-            (
-                SubList<face>(subFaces, subFaces.size()),
-                pp.points()
-            );
-
-            if (debug & 2)
-            {
-                Pout<< "distributePatches: to processor " << domain
-                    << " sending faces " << subPatch.faceCentres() << endl;
-            }
-
-            UOPstream toDomain(domain, pBufs);
-            toDomain
-                << subPatch.localFaces() << subPatch.localPoints()
-                << globalElems;
-        }
-    }
-
-    // Start receiving
-    pBufs.finishedSends();
-
-    faces.setSize(Pstream::nProcs());
-    points.setSize(Pstream::nProcs());
-    faceIDs.setSize(Pstream::nProcs());
-
-    {
-        // Set up 'send' to myself
-        const labelList& sendElems = map.subMap()[Pstream::myProcNo()];
-        faceList subFaces(UIndirectList<face>(pp, sendElems));
-        primitivePatch subPatch
-        (
-            SubList<face>(subFaces, subFaces.size()),
-            pp.points()
-        );
-
-        // Receive
-        if (debug & 2)
-        {
-            Pout<< "distributePatches: to processor " << Pstream::myProcNo()
-                << " sending faces " << subPatch.faceCentres() << endl;
-        }
-
-        faces[Pstream::myProcNo()] = subPatch.localFaces();
-        points[Pstream::myProcNo()] = subPatch.localPoints();
-
-        faceIDs[Pstream::myProcNo()].setSize(sendElems.size());
-        forAll(sendElems, i)
-        {
-            faceIDs[Pstream::myProcNo()][i] = gi.toGlobal(sendElems[i]);
-        }
-    }
-
-    // Consume
-    for (label domain = 0; domain < Pstream::nProcs(); domain++)
-    {
-        const labelList& recvElems = map.constructMap()[domain];
-
-        if (domain != Pstream::myProcNo() && recvElems.size())
-        {
-            UIPstream str(domain, pBufs);
-
-            str >> faces[domain]
-                >> points[domain]
-                >> faceIDs[domain];
-        }
-    }
-}
-
-
-template<class SourcePatch, class TargetPatch>
-void Foam::AMIInterpolation<SourcePatch, TargetPatch>::
-distributeAndMergePatches
-(
-    const mapDistribute& map,
-    const TargetPatch& tgtPatch,
-    const globalIndex& gi,
-    faceList& tgtFaces,
-    pointField& tgtPoints,
-    labelList& tgtFaceIDs
-) const
-{
-    // Exchange per-processor data
-    List<faceList> allFaces;
-    List<pointField> allPoints;
-    List<labelList> allTgtFaceIDs;
-    distributePatches(map, tgtPatch, gi, allFaces, allPoints, allTgtFaceIDs);
-
-    // Renumber and flatten
-    label nFaces = 0;
-    label nPoints = 0;
-    forAll(allFaces, procI)
-    {
-        nFaces += allFaces[procI].size();
-        nPoints += allPoints[procI].size();
-    }
-
-    tgtFaces.setSize(nFaces);
-    tgtPoints.setSize(nPoints);
-    tgtFaceIDs.setSize(nFaces);
-
-    nFaces = 0;
-    nPoints = 0;
-
-    // My own data first
-    {
-        const labelList& faceIDs = allTgtFaceIDs[Pstream::myProcNo()];
-        SubList<label>(tgtFaceIDs, faceIDs.size()).assign(faceIDs);
-
-        const faceList& fcs = allFaces[Pstream::myProcNo()];
-        forAll(fcs, i)
-        {
-            const face& f = fcs[i];
-            face& newF = tgtFaces[nFaces++];
-            newF.setSize(f.size());
-            forAll(f, fp)
-            {
-                newF[fp] = f[fp] + nPoints;
-            }
-        }
-
-        const pointField& pts = allPoints[Pstream::myProcNo()];
-        forAll(pts, i)
-        {
-            tgtPoints[nPoints++] = pts[i];
-        }
-    }
-
-
-    // Other proc data follows
-    forAll(allFaces, procI)
-    {
-        if (procI != Pstream::myProcNo())
-        {
-            const labelList& faceIDs = allTgtFaceIDs[procI];
-            SubList<label>(tgtFaceIDs, faceIDs.size(), nFaces).assign(faceIDs);
-
-            const faceList& fcs = allFaces[procI];
-            forAll(fcs, i)
-            {
-                const face& f = fcs[i];
-                face& newF = tgtFaces[nFaces++];
-                newF.setSize(f.size());
-                forAll(f, fp)
-                {
-                    newF[fp] = f[fp] + nPoints;
-                }
-            }
-
-            const pointField& pts = allPoints[procI];
-            forAll(pts, i)
-            {
-                tgtPoints[nPoints++] = pts[i];
-            }
-        }
-    }
-
-    // Merge
-    labelList oldToNew;
-    pointField newTgtPoints;
-    bool hasMerged = mergePoints
-    (
-        tgtPoints,
-        SMALL,
-        false,
-        oldToNew,
-        newTgtPoints
-    );
-
-    if (hasMerged)
-    {
-        if (debug)
-        {
-            Pout<< "Merged from " << tgtPoints.size()
-                << " down to " << newTgtPoints.size() << " points" << endl;
-        }
-
-        tgtPoints.transfer(newTgtPoints);
-        forAll(tgtFaces, i)
-        {
-            inplaceRenumber(oldToNew, tgtFaces[i]);
-        }
-    }
-}
-
-
-template<class SourcePatch, class TargetPatch>
-Foam::autoPtr<Foam::mapDistribute>
-Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcProcMap
-(
-    const SourcePatch& srcPatch,
-    const TargetPatch& tgtPatch
-) const
-{
-    // Get decomposition of patch
-    List<treeBoundBoxList> procBb(Pstream::nProcs());
-
-    if (srcPatch.size())
-    {
-        procBb[Pstream::myProcNo()] = treeBoundBoxList
-        (
-            1,  // For now single bounding box per proc
-            treeBoundBox
-            (
-                srcPatch.points(),
-                srcPatch.meshPoints()
-            )
-        );
-    }
-    else
-    {
-        procBb[Pstream::myProcNo()] = treeBoundBoxList();
-    }
-
-    // slightly increase size of bounding boxes to allow for cases where
-    // bounding boxes are perfectly alligned
-    forAll(procBb[Pstream::myProcNo()], bbI)
-    {
-        treeBoundBox& bb = procBb[Pstream::myProcNo()][bbI];
-        bb.inflate(0.01);
-    }
-
-    Pstream::gatherList(procBb);
-    Pstream::scatterList(procBb);
-
-
-    if (debug)
-    {
-        Info<< "Determining extent of srcPatch per processor:" << nl
-            << "\tproc\tbb" << endl;
-        forAll(procBb, procI)
-        {
-            Info<< '\t' << procI << '\t' << procBb[procI] << endl;
-        }
-    }
-
-
-    // Determine which faces of tgtPatch overlaps srcPatch per proc
-    const faceList& faces = tgtPatch.localFaces();
-    const pointField& points = tgtPatch.localPoints();
-
-    labelListList sendMap;
-
-    {
-        // Per processor indices into all segments to send
-        List<DynamicList<label> > dynSendMap(Pstream::nProcs());
-
-        // Work array - whether processor bb overlaps the face bounds
-        boolList procBbOverlaps(Pstream::nProcs());
-
-        forAll(faces, faceI)
-        {
-            if (faces[faceI].size())
-            {
-                treeBoundBox faceBb(points, faces[faceI]);
-
-                // Find the processor this face overlaps
-                calcOverlappingProcs(procBb, faceBb, procBbOverlaps);
-
-                forAll(procBbOverlaps, procI)
-                {
-                    if (procBbOverlaps[procI])
-                    {
-                        dynSendMap[procI].append(faceI);
-                    }
-                }
-            }
-        }
-
-        // Convert dynamicList to labelList
-        sendMap.setSize(Pstream::nProcs());
-        forAll(sendMap, procI)
-        {
-            sendMap[procI].transfer(dynSendMap[procI]);
-        }
-    }
-
-    // Debug printing
-    if (debug)
-    {
-        Pout<< "Of my " << faces.size() << " I need to send to:" << nl
-            << "\tproc\tfaces" << endl;
-        forAll(sendMap, procI)
-        {
-            Pout<< '\t' << procI << '\t' << sendMap[procI].size() << endl;
-        }
-    }
-
-
-    // Send over how many faces I need to receive
-    labelListList sendSizes(Pstream::nProcs());
-    sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
-    forAll(sendMap, procI)
-    {
-        sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
-    }
-    Pstream::gatherList(sendSizes);
-    Pstream::scatterList(sendSizes);
-
-
-    // Determine order of receiving
-    labelListList constructMap(Pstream::nProcs());
-
-    // My local segment first
-    constructMap[Pstream::myProcNo()] = identity
-    (
-        sendMap[Pstream::myProcNo()].size()
-    );
-
-    label segmentI = constructMap[Pstream::myProcNo()].size();
-    forAll(constructMap, procI)
-    {
-        if (procI != Pstream::myProcNo())
-        {
-            // What I need to receive is what other processor is sending to me
-            label nRecv = sendSizes[procI][Pstream::myProcNo()];
-            constructMap[procI].setSize(nRecv);
-
-            for (label i = 0; i < nRecv; i++)
-            {
-                constructMap[procI][i] = segmentI++;
-            }
-        }
-    }
-
-    autoPtr<mapDistribute> mapPtr
-    (
-        new mapDistribute
-        (
-            segmentI,       // size after construction
-            sendMap.xfer(),
-            constructMap.xfer()
-        )
-    );
-
-    return mapPtr;
-}
-
-
 template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::projectPointsToSurface
 (
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
index 5b91afc38a..abf4523afd 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
@@ -554,6 +554,7 @@ public:
 
 #ifdef NoRepository
 #   include "AMIInterpolation.C"
+#   include "AMIInterpolationParallelOps.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationParallelOps.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationParallelOps.C
new file mode 100644
index 0000000000..5abe0778eb
--- /dev/null
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationParallelOps.C
@@ -0,0 +1,481 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "AMIInterpolation.H"
+#include "mergePoints.H"
+#include "mapDistribute.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class SourcePatch, class TargetPatch>
+Foam::label Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcDistribution
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch
+) const
+{
+    label procI = 0;
+
+    if (Pstream::parRun())
+    {
+        List<label> facesPresentOnProc(Pstream::nProcs(), 0);
+        if ((srcPatch.size() > 0) || (tgtPatch.size() > 0))
+        {
+            facesPresentOnProc[Pstream::myProcNo()] = 1;
+        }
+        else
+        {
+            facesPresentOnProc[Pstream::myProcNo()] = 0;
+        }
+
+        Pstream::gatherList(facesPresentOnProc);
+        Pstream::scatterList(facesPresentOnProc);
+
+        label nHaveFaces = sum(facesPresentOnProc);
+
+        if (nHaveFaces > 1)
+        {
+            procI = -1;
+            if (debug)
+            {
+                Info<< "AMIInterpolation::calcDistribution: "
+                    << "AMI split across multiple processors" << endl;
+            }
+        }
+        else if (nHaveFaces == 1)
+        {
+            procI = findIndex(facesPresentOnProc, 1);
+            if (debug)
+            {
+                Info<< "AMIInterpolation::calcDistribution: "
+                    << "AMI local to processor" << procI << endl;
+            }
+        }
+    }
+
+
+    // Either not parallel or no faces on any processor
+    return procI;
+}
+
+
+template<class SourcePatch, class TargetPatch>
+Foam::label
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcOverlappingProcs
+(
+    const List<treeBoundBoxList>& procBb,
+    const treeBoundBox& bb,
+    boolList& overlaps
+) const
+{
+    overlaps.setSize(procBb.size());
+    overlaps = false;
+
+    label nOverlaps = 0;
+
+    forAll(procBb, procI)
+    {
+        const List<treeBoundBox>& bbs = procBb[procI];
+
+        forAll(bbs, bbI)
+        {
+            if (bbs[bbI].overlaps(bb))
+            {
+                overlaps[procI] = true;
+                nOverlaps++;
+                break;
+            }
+        }
+    }
+    return nOverlaps;
+}
+
+
+template<class SourcePatch, class TargetPatch>
+void Foam::AMIInterpolation<SourcePatch, TargetPatch>::distributePatches
+(
+    const mapDistribute& map,
+    const TargetPatch& pp,
+    const globalIndex& gi,
+    List<faceList>& faces,
+    List<pointField>& points,
+    List<labelList>& faceIDs
+) const
+{
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& sendElems = map.subMap()[domain];
+
+        if (domain != Pstream::myProcNo() && sendElems.size())
+        {
+            labelList globalElems(sendElems.size());
+            forAll(sendElems, i)
+            {
+                globalElems[i] = gi.toGlobal(sendElems[i]);
+            }
+
+            faceList subFaces(UIndirectList<face>(pp, sendElems));
+            primitivePatch subPatch
+            (
+                SubList<face>(subFaces, subFaces.size()),
+                pp.points()
+            );
+
+            if (debug & 2)
+            {
+                Pout<< "distributePatches: to processor " << domain
+                    << " sending faces " << subPatch.faceCentres() << endl;
+            }
+
+            UOPstream toDomain(domain, pBufs);
+            toDomain
+                << subPatch.localFaces() << subPatch.localPoints()
+                << globalElems;
+        }
+    }
+
+    // Start receiving
+    pBufs.finishedSends();
+
+    faces.setSize(Pstream::nProcs());
+    points.setSize(Pstream::nProcs());
+    faceIDs.setSize(Pstream::nProcs());
+
+    {
+        // Set up 'send' to myself
+        const labelList& sendElems = map.subMap()[Pstream::myProcNo()];
+        faceList subFaces(UIndirectList<face>(pp, sendElems));
+        primitivePatch subPatch
+        (
+            SubList<face>(subFaces, subFaces.size()),
+            pp.points()
+        );
+
+        // Receive
+        if (debug & 2)
+        {
+            Pout<< "distributePatches: to processor " << Pstream::myProcNo()
+                << " sending faces " << subPatch.faceCentres() << endl;
+        }
+
+        faces[Pstream::myProcNo()] = subPatch.localFaces();
+        points[Pstream::myProcNo()] = subPatch.localPoints();
+
+        faceIDs[Pstream::myProcNo()].setSize(sendElems.size());
+        forAll(sendElems, i)
+        {
+            faceIDs[Pstream::myProcNo()][i] = gi.toGlobal(sendElems[i]);
+        }
+    }
+
+    // Consume
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& recvElems = map.constructMap()[domain];
+
+        if (domain != Pstream::myProcNo() && recvElems.size())
+        {
+            UIPstream str(domain, pBufs);
+
+            str >> faces[domain]
+                >> points[domain]
+                >> faceIDs[domain];
+        }
+    }
+}
+
+
+template<class SourcePatch, class TargetPatch>
+void Foam::AMIInterpolation<SourcePatch, TargetPatch>::
+distributeAndMergePatches
+(
+    const mapDistribute& map,
+    const TargetPatch& tgtPatch,
+    const globalIndex& gi,
+    faceList& tgtFaces,
+    pointField& tgtPoints,
+    labelList& tgtFaceIDs
+) const
+{
+    // Exchange per-processor data
+    List<faceList> allFaces;
+    List<pointField> allPoints;
+    List<labelList> allTgtFaceIDs;
+    distributePatches(map, tgtPatch, gi, allFaces, allPoints, allTgtFaceIDs);
+
+    // Renumber and flatten
+    label nFaces = 0;
+    label nPoints = 0;
+    forAll(allFaces, procI)
+    {
+        nFaces += allFaces[procI].size();
+        nPoints += allPoints[procI].size();
+    }
+
+    tgtFaces.setSize(nFaces);
+    tgtPoints.setSize(nPoints);
+    tgtFaceIDs.setSize(nFaces);
+
+    nFaces = 0;
+    nPoints = 0;
+
+    // My own data first
+    {
+        const labelList& faceIDs = allTgtFaceIDs[Pstream::myProcNo()];
+        SubList<label>(tgtFaceIDs, faceIDs.size()).assign(faceIDs);
+
+        const faceList& fcs = allFaces[Pstream::myProcNo()];
+        forAll(fcs, i)
+        {
+            const face& f = fcs[i];
+            face& newF = tgtFaces[nFaces++];
+            newF.setSize(f.size());
+            forAll(f, fp)
+            {
+                newF[fp] = f[fp] + nPoints;
+            }
+        }
+
+        const pointField& pts = allPoints[Pstream::myProcNo()];
+        forAll(pts, i)
+        {
+            tgtPoints[nPoints++] = pts[i];
+        }
+    }
+
+
+    // Other proc data follows
+    forAll(allFaces, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            const labelList& faceIDs = allTgtFaceIDs[procI];
+            SubList<label>(tgtFaceIDs, faceIDs.size(), nFaces).assign(faceIDs);
+
+            const faceList& fcs = allFaces[procI];
+            forAll(fcs, i)
+            {
+                const face& f = fcs[i];
+                face& newF = tgtFaces[nFaces++];
+                newF.setSize(f.size());
+                forAll(f, fp)
+                {
+                    newF[fp] = f[fp] + nPoints;
+                }
+            }
+
+            const pointField& pts = allPoints[procI];
+            forAll(pts, i)
+            {
+                tgtPoints[nPoints++] = pts[i];
+            }
+        }
+    }
+
+    // Merge
+    labelList oldToNew;
+    pointField newTgtPoints;
+    bool hasMerged = mergePoints
+    (
+        tgtPoints,
+        SMALL,
+        false,
+        oldToNew,
+        newTgtPoints
+    );
+
+    if (hasMerged)
+    {
+        if (debug)
+        {
+            Pout<< "Merged from " << tgtPoints.size()
+                << " down to " << newTgtPoints.size() << " points" << endl;
+        }
+
+        tgtPoints.transfer(newTgtPoints);
+        forAll(tgtFaces, i)
+        {
+            inplaceRenumber(oldToNew, tgtFaces[i]);
+        }
+    }
+}
+
+
+template<class SourcePatch, class TargetPatch>
+Foam::autoPtr<Foam::mapDistribute>
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcProcMap
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch
+) const
+{
+    // Get decomposition of patch
+    List<treeBoundBoxList> procBb(Pstream::nProcs());
+
+    if (srcPatch.size())
+    {
+        procBb[Pstream::myProcNo()] = treeBoundBoxList
+        (
+            1,  // For now single bounding box per proc
+            treeBoundBox
+            (
+                srcPatch.points(),
+                srcPatch.meshPoints()
+            )
+        );
+    }
+    else
+    {
+        procBb[Pstream::myProcNo()] = treeBoundBoxList();
+    }
+
+    // slightly increase size of bounding boxes to allow for cases where
+    // bounding boxes are perfectly alligned
+    forAll(procBb[Pstream::myProcNo()], bbI)
+    {
+        treeBoundBox& bb = procBb[Pstream::myProcNo()][bbI];
+        bb.inflate(0.01);
+    }
+
+    Pstream::gatherList(procBb);
+    Pstream::scatterList(procBb);
+
+
+    if (debug)
+    {
+        Info<< "Determining extent of srcPatch per processor:" << nl
+            << "\tproc\tbb" << endl;
+        forAll(procBb, procI)
+        {
+            Info<< '\t' << procI << '\t' << procBb[procI] << endl;
+        }
+    }
+
+
+    // Determine which faces of tgtPatch overlaps srcPatch per proc
+    const faceList& faces = tgtPatch.localFaces();
+    const pointField& points = tgtPatch.localPoints();
+
+    labelListList sendMap;
+
+    {
+        // Per processor indices into all segments to send
+        List<DynamicList<label> > dynSendMap(Pstream::nProcs());
+
+        // Work array - whether processor bb overlaps the face bounds
+        boolList procBbOverlaps(Pstream::nProcs());
+
+        forAll(faces, faceI)
+        {
+            if (faces[faceI].size())
+            {
+                treeBoundBox faceBb(points, faces[faceI]);
+
+                // Find the processor this face overlaps
+                calcOverlappingProcs(procBb, faceBb, procBbOverlaps);
+
+                forAll(procBbOverlaps, procI)
+                {
+                    if (procBbOverlaps[procI])
+                    {
+                        dynSendMap[procI].append(faceI);
+                    }
+                }
+            }
+        }
+
+        // Convert dynamicList to labelList
+        sendMap.setSize(Pstream::nProcs());
+        forAll(sendMap, procI)
+        {
+            sendMap[procI].transfer(dynSendMap[procI]);
+        }
+    }
+
+    // Debug printing
+    if (debug)
+    {
+        Pout<< "Of my " << faces.size() << " I need to send to:" << nl
+            << "\tproc\tfaces" << endl;
+        forAll(sendMap, procI)
+        {
+            Pout<< '\t' << procI << '\t' << sendMap[procI].size() << endl;
+        }
+    }
+
+
+    // Send over how many faces I need to receive
+    labelListList sendSizes(Pstream::nProcs());
+    sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
+    forAll(sendMap, procI)
+    {
+        sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
+    }
+    Pstream::gatherList(sendSizes);
+    Pstream::scatterList(sendSizes);
+
+
+    // Determine order of receiving
+    labelListList constructMap(Pstream::nProcs());
+
+    // My local segment first
+    constructMap[Pstream::myProcNo()] = identity
+    (
+        sendMap[Pstream::myProcNo()].size()
+    );
+
+    label segmentI = constructMap[Pstream::myProcNo()].size();
+    forAll(constructMap, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            // What I need to receive is what other processor is sending to me
+            label nRecv = sendSizes[procI][Pstream::myProcNo()];
+            constructMap[procI].setSize(nRecv);
+
+            for (label i = 0; i < nRecv; i++)
+            {
+                constructMap[procI][i] = segmentI++;
+            }
+        }
+    }
+
+    autoPtr<mapDistribute> mapPtr
+    (
+        new mapDistribute
+        (
+            segmentI,       // size after construction
+            sendMap.xfer(),
+            constructMap.xfer()
+        )
+    );
+
+    return mapPtr;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleDyMFoam/propeller/Allclean b/tutorials/incompressible/pimpleDyMFoam/propeller/Allclean
index f87f062960..7be04a3bab 100755
--- a/tutorials/incompressible/pimpleDyMFoam/propeller/Allclean
+++ b/tutorials/incompressible/pimpleDyMFoam/propeller/Allclean
@@ -4,9 +4,11 @@ cd ${0%/*} || exit 1    # run from this directory
 # Source tutorial clean functions
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
+# remove surface
+\rm -f constant/triSurface/propellerTip.obj.gz
+
 cleanCase
 
 \rm -rf 0
 \rm -rf constant/extendedFeatureEdgeMesh/
 \rm -f constant/triSurface/*.eMesh*
-