diff --git a/applications/test/Polynomial/Make/files b/applications/test/Polynomial/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..598fccfecfa73a3dceee6885e9e3d5230af38a97
--- /dev/null
+++ b/applications/test/Polynomial/Make/files
@@ -0,0 +1,3 @@
+PolynomialTest.C
+
+EXE = $(FOAM_USER_APPBIN)/PolynomialTest
diff --git a/applications/test/Polynomial/Make/options b/applications/test/Polynomial/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..1b983985770833401ab9dbd2bbce72c2c867c0c8
--- /dev/null
+++ b/applications/test/Polynomial/Make/options
@@ -0,0 +1,3 @@
+EXE_INC = \
+
+EXE_LIBS = \
diff --git a/applications/test/Polynomial/PolynomialTest.C b/applications/test/Polynomial/PolynomialTest.C
new file mode 100644
index 0000000000000000000000000000000000000000..e3dfd2ea4c3372c14bb822accd0842ef30d7b5b0
--- /dev/null
+++ b/applications/test/Polynomial/PolynomialTest.C
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    PolynomialTest
+
+Description
+    Test application for the templated Polynomial class
+
+\*---------------------------------------------------------------------------*/
+
+#include "IFstream.H"
+#include "Polynomial.H"
+#include "Random.H"
+
+using namespace Foam;
+
+scalar polyValue(const scalar x)
+{
+    // Hard-coded polynomial 8 coeff (7th order)
+    return
+        0.11
+      + 0.45*x
+      - 0.94*sqr(x)
+      + 1.58*pow3(x)
+      - 2.58*pow4(x)
+      + 0.08*pow5(x)
+      + 3.15*pow6(x)
+      - 4.78*x*pow6(x);
+}
+
+
+scalar intPolyValue(const scalar x)
+{
+    // Hard-coded integrated form of above polynomial
+    return
+        0.11*x
+      + 0.45/2.0*sqr(x)
+      - 0.94/3.0*pow3(x)
+      + 1.58/4.0*pow4(x)
+      - 2.58/5.0*pow5(x)
+      + 0.08/6.0*pow6(x)
+      + 3.15/7.0*x*pow6(x)
+      - 4.78/8.0*x*x*pow6(x);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    IFstream is("polyTestInput");
+
+    Polynomial<8> poly("testPoly", is);
+    Polynomial<9> intPoly(poly.integrate(0.0));
+
+    Info<< "poly = " << poly << endl;
+    Info<< "intPoly = " << intPoly << nl << endl;
+
+    Info<< "2*poly = " << 2*poly << endl;
+    Info<< "poly+poly = " << poly + poly << nl << endl;
+
+    Info<< "3*poly = " << 3*poly << endl;
+    Info<< "poly+poly+poly = " << poly + poly + poly << nl << endl;
+
+    Info<< "3*poly - 2*poly = " << 3*poly - 2*poly << nl << endl;
+
+    Polynomial<8> polyCopy = poly;
+    Info<< "poly, polyCopy = " << poly << ", " << polyCopy << nl << endl;
+    polyCopy = 2.5*poly;
+    Info<< "2.5*polyCopy = " << polyCopy << nl << endl;
+
+    Random rnd(123456);
+    for (int i=0; i<10; i++)
+    {
+        scalar x = rnd.scalar01()*100;
+
+        scalar px = polyValue(x);
+        scalar ipx = intPolyValue(x);
+
+        scalar pxTest = poly.evaluate(x);
+        scalar ipxTest = intPoly.evaluate(x);
+
+        Info<<"\nx = " << x << endl;
+        Info<< "    px, pxTest = " << px << ", " << pxTest << endl;
+        Info<< "    ipx, ipxTest = " << ipx << ", " << ipxTest << endl;
+
+        if (mag(px - pxTest) > SMALL)
+        {
+            Info<< "    *** WARNING: px != pxTest: " << px - pxTest << endl;
+        }
+
+        if (mag(ipx - ipxTest) > SMALL)
+        {
+            Info<< "    *** WARNING: ipx != ipxTest: " << ipx - ipxTest << endl;
+        }
+
+        Info<< endl;
+    }
+
+    Info<< nl << "Done." << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/Polynomial/polyTestInput b/applications/test/Polynomial/polyTestInput
new file mode 100644
index 0000000000000000000000000000000000000000..4ba8f65e51f4dc8a68381d3e9af7bbf30a1d2a39
--- /dev/null
+++ b/applications/test/Polynomial/polyTestInput
@@ -0,0 +1,11 @@
+    testPoly
+    (
+        0.11
+        0.45
+       -0.94
+        1.58
+       -2.58
+        0.08
+        3.15
+       -4.78
+    )
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 6729cab4c2e9932765ee6bc60bcd9034cd19b160..524b4d00d5d103820851e8ee993c1185fa329935 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -344,7 +344,7 @@ void Foam::ReactingParcel<ParcelType>::calc
             td.cloud().hcTrans()[cellI] +=
                 np0
                *dMassPC[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].H(T0);
+               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
 
         // Update momentum transfer
@@ -423,6 +423,11 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     scalarField& Cs
 )
 {
+    typedef PhaseChangeModel
+    <
+        typename ReactingParcel<ParcelType>::trackData::cloudType
+    > phaseChangeModelType;
+
     if
     (
         !td.cloud().phaseChange().active()
@@ -464,17 +469,26 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
             td.cloud().composition().localToGlobalCarrierId(idPhase, i);
         const label idl = td.cloud().composition().globalIds(idPhase)[i];
 
-        const scalar hv = td.cloud().mcCarrierThermo().speciesData()[idc].H(Ts);
-        const scalar hl =
-            td.cloud().composition().liquids().properties()[idl].h(pc_, Ts);
+        if
+        (
+            td.cloud().phaseChange().enthalpyTransfer()
+         == phaseChangeModelType::etLatentHeat
+        )
+        {
+            scalar hlp =
+                td.cloud().composition().liquids().properties()[idl].hl(pc_, T);
 
-        // Enthalphy transfer to carrier phase - method 1 using enthalpy diff
-        Sh += dMassPC[i]*(hl - hv)/dt;
+            Sh -= dMassPC[i]*hlp/dt;
+        }
+        else
+        {
+            // Note: enthalpies of both phases must use the same reference
+            scalar hc = td.cloud().mcCarrierThermo().speciesData()[idc].H(T);
+            scalar hp =
+                td.cloud().composition().liquids().properties()[idl].h(pc_, T);
 
-        // Enthalphy transfer to carrier phase - method 2 using latent heat
-//        const scalar hl =
-//            td.cloud().composition().liquids().properties()[idl].hl(pc_, Ts);
-//        Sh -= dMassPC[i]*hl/dt;
+            Sh -= dMassPC[i]*(hc - hp)/dt;
+        }
 
         // Update particle surface thermo properties
         const scalar Dab =
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index 8219ec2ac7b8f49495c081074cc1e50cb94730da..41d7c6d30c224fe87011cee76cadb0a583d99eac 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -136,6 +136,9 @@ public:
 
     public:
 
+        typedef ReactingCloud<ParcelType> cloudType;
+
+
         // Constructors
 
             //- Construct from components
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C
index c60f1bf903871b32a9a6602f56e0267334c0ae5d..c88e07ae19622177aacd6dcb2ec3c083cf29f0d4 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C
@@ -26,6 +26,48 @@ License
 
 #include "PhaseChangeModel.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class CloudType>
+const Foam::wordList Foam::PhaseChangeModel<CloudType>::
+enthalpyTransferTypeNames
+(
+    IStringStream
+    (
+        "("
+            "latentHeat "
+            "enthalpyDifference"
+        ")"
+    )()
+);
+
+
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+template<class CloudType>
+typename Foam::PhaseChangeModel<CloudType>::enthalpyTransferType
+Foam::PhaseChangeModel<CloudType>::wordToEnthalpyTransfer(const word& etName)
+const
+{
+    forAll(enthalpyTransferTypeNames, i)
+    {
+        if (etName == enthalpyTransferTypeNames[i])
+        {
+            return enthalpyTransferType(i);
+        }
+    }
+
+    FatalErrorIn
+    (
+        "PhaseChangeModel<CloudType>::enthalpyTransferType"
+        "PhaseChangeModel<CloudType>::wordToEnthalpyTransfer(const word&) const"
+    )   << "Unknown enthalpyType " << etName << ". Valid selections are:" << nl
+        << enthalpyTransferTypeNames << exit(FatalError);
+
+    return enthalpyTransferType(0);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -36,7 +78,8 @@ Foam::PhaseChangeModel<CloudType>::PhaseChangeModel
 :
     dict_(dictionary::null),
     owner_(owner),
-    coeffDict_(dictionary::null)
+    coeffDict_(dictionary::null),
+    enthalpyTransfer_(etLatentHeat)
 {}
 
 
@@ -50,7 +93,11 @@ Foam::PhaseChangeModel<CloudType>::PhaseChangeModel
 :
     dict_(dict),
     owner_(owner),
-    coeffDict_(dict.subDict(type + "Coeffs"))
+    coeffDict_(dict.subDict(type + "Coeffs")),
+    enthalpyTransfer_
+    (
+        wordToEnthalpyTransfer(coeffDict_.lookup("enthalpyTransfer"))
+    )
 {}
 
 
@@ -83,6 +130,14 @@ const Foam::dictionary& Foam::PhaseChangeModel<CloudType>::coeffDict() const
 }
 
 
+template<class CloudType>
+const typename Foam::PhaseChangeModel<CloudType>::enthalpyTransferType&
+Foam::PhaseChangeModel<CloudType>::enthalpyTransfer() const
+{
+    return enthalpyTransfer_;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "NewPhaseChangeModel.C"
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
index e19856a0648c83c7cf70386ad7e3c60f0fd01f4d..27b92e2539119525271223a6dcb5e9ee6f33869d 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
@@ -53,6 +53,21 @@ namespace Foam
 template<class CloudType>
 class PhaseChangeModel
 {
+public:
+
+    // Public enumerations
+
+        //- Enthalpy transfer type
+        enum enthalpyTransferType
+        {
+            etLatentHeat,
+            etEnthalpyDifference
+        };
+
+        //- Name representations of enthalpy transfer types
+        static const Foam::wordList enthalpyTransferTypeNames;
+
+
 protected:
 
     // Protected data
@@ -66,9 +81,15 @@ protected:
         //- The coefficient dictionary
         const dictionary coeffDict_;
 
+        //- Enthalpy transfer type enumeration
+        enthalpyTransferType enthalpyTransfer_;
+
 
     // Protected member functions
 
+        //- Convert word to enthalpy transfer type
+        enthalpyTransferType wordToEnthalpyTransfer(const word& etName) const;
+
         //- Sherwood number
         scalar Sh() const;
 
@@ -129,6 +150,9 @@ public:
         //- Return the coefficient dictionary
         const dictionary& coeffDict() const;
 
+        //- Return the enthalpy transfer type enumeration
+        const enthalpyTransferType& enthalpyTransfer() const;
+
 
     // Member Functions
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C
index 3e204e65b1dea330e3840ebc99c12d172c9e1ddb..854ba279b6b4d0acddde2af1393f749be62b8e36 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C
@@ -449,7 +449,7 @@ template<class CompType, class ThermoType>
 Foam::tmp<Foam::volScalarField>
 Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
 {
-    scalar pf,cf,pr,cr;
+    scalar pf, cf, pr, cr;
     label lRef, rRef;
 
     const volScalarField rho
diff --git a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.C b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.C
index 9fd15c2986c5fe3365f6c6d65b8af2bb2198aa34..fcea63c7bbc43ebefa27d2ade428cd27cf9c3b93 100644
--- a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.C
+++ b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.C
@@ -39,7 +39,9 @@ icoPolynomial<PolySize>::icoPolynomial(Istream& is)
 :
     specie(is),
     rhoPolynomial_("rhoPolynomial", is)
-{}
+{
+    rhoPolynomial_ *= this->W();
+}
 
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
@@ -47,7 +49,8 @@ icoPolynomial<PolySize>::icoPolynomial(Istream& is)
 template<int PolySize>
 Ostream& operator<<(Ostream& os, const icoPolynomial<PolySize>& ip)
 {
-    os  << static_cast<const specie&>(ip);
+    os  << static_cast<const specie&>(ip) << tab
+        << "rhoPolynomial" << tab << ip.rhoPolynomial_/ip.W();
 
     os.check
     (
diff --git a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H
index 3c66e03d4232c1836b66151a49a8506edf11fc27..6dbea13cd72df244edc78286bb9be7fc6011e6fa 100644
--- a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H
+++ b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H
@@ -99,7 +99,8 @@ class icoPolynomial
 {
     // Private data
 
-        //- Density [kg/m^3]
+        //- Density
+        //  Note: input in [kg/m3], but internally uses [kg/m3/kmol]
         Polynomial<PolySize> rhoPolynomial_;
 
 
@@ -117,6 +118,9 @@ public:
         //- Construct from Istream
         icoPolynomial(Istream&);
 
+        //- Construct as copy
+        inline icoPolynomial(const icoPolynomial&);
+
         //- Construct as named copy
         inline icoPolynomial(const word& name, const icoPolynomial&);
 
@@ -141,6 +145,7 @@ public:
 
     // Member operators
 
+        inline icoPolynomial& operator=(const icoPolynomial&);
         inline void operator+=(const icoPolynomial&);
         inline void operator-=(const icoPolynomial&);
 
diff --git a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H
index d174c57192bac3d842188fc6a03ddda8eeb23aac..d85fd4a2a609fc38489342ab4ae7b19e32937f8f 100644
--- a/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H
+++ b/src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H
@@ -42,6 +42,17 @@ inline Foam::icoPolynomial<PolySize>::icoPolynomial
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+template<int PolySize>
+inline Foam::icoPolynomial<PolySize>::icoPolynomial
+(
+    const icoPolynomial<PolySize>& ip
+)
+:
+    specie(ip),
+    rhoPolynomial_(ip.rhoPolynomial_)
+{}
+
+
 template<int PolySize>
 inline Foam::icoPolynomial<PolySize>::icoPolynomial
 (
@@ -78,7 +89,7 @@ Foam::icoPolynomial<PolySize>::New(Istream& is)
 template<int PolySize>
 inline Foam::scalar Foam::icoPolynomial<PolySize>::rho(scalar, scalar T) const
 {
-    return rhoPolynomial_.evaluate(T);
+    return rhoPolynomial_.evaluate(T)/this->W();
 }
 
 
@@ -98,6 +109,20 @@ inline Foam::scalar Foam::icoPolynomial<PolySize>::Z(scalar, scalar) const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
+template<int PolySize>
+inline Foam::icoPolynomial<PolySize>& Foam::icoPolynomial<PolySize>::operator=
+(
+    const icoPolynomial<PolySize>& ip
+)
+{
+    specie::operator=(ip);
+
+    rhoPolynomial_ = ip.rhoPolynomial_;
+
+    return *this;
+}
+
+
 template<int PolySize>
 inline void Foam::icoPolynomial<PolySize>::operator+=
 (
@@ -122,6 +147,7 @@ inline void Foam::icoPolynomial<PolySize>::operator-=
 )
 {
     scalar molr1 = this->nMoles();
+
     specie::operator-=(ip);
 
     molr1 /= this->nMoles();
@@ -147,15 +173,15 @@ Foam::icoPolynomial<PolySize> Foam::operator+
     const icoPolynomial<PolySize>& ip2
 )
 {
-    scalar mol1 = ip1.nMoles();
-    scalar mol2 = ip2.nMoles();
-    scalar nMoles = mol1 + mol2;
+    scalar nMoles = ip1.nMoles() + ip2.nMoles();
+    scalar molr1 = ip1.nMoles()/nMoles;
+    scalar molr2 = ip2.nMoles()/nMoles;
 
     return icoPolynomial<PolySize>
     (
         static_cast<const specie&>(ip1)
       + static_cast<const specie&>(ip2),
-        (mol1/nMoles)*ip1.rhoPolynomial_ + (mol2/nMoles)*ip2.rhoPolynomial_
+        molr1*ip1.rhoPolynomial_ + molr2*ip2.rhoPolynomial_
     );
 }
 
@@ -167,15 +193,15 @@ Foam::icoPolynomial<PolySize> Foam::operator-
     const icoPolynomial<PolySize>& ip2
 )
 {
-    scalar mol1 = ip1.nMoles();
-    scalar mol2 = ip2.nMoles();
-    scalar nMoles = mol1 + mol2;
+    scalar nMoles = ip1.nMoles() + ip2.nMoles();
+    scalar molr1 = ip1.nMoles()/nMoles;
+    scalar molr2 = ip2.nMoles()/nMoles;
 
     return icoPolynomial<PolySize>
     (
         static_cast<const specie&>(ip1)
       - static_cast<const specie&>(ip2),
-        (mol1/nMoles)*ip1.rhoPolynomial_ - (mol2/nMoles)*ip2.rhoPolynomial_
+        molr1*ip1.rhoPolynomial_ - molr2*ip2.rhoPolynomial_
     );
 }
 
diff --git a/src/thermophysicalModels/specie/specie/specieI.H b/src/thermophysicalModels/specie/specie/specieI.H
index 0965f0f7cd093d5dd6a17bd450ede63799966d48..a2986a67d4f3882bcf5bc600e0a877df47ad717e 100644
--- a/src/thermophysicalModels/specie/specie/specieI.H
+++ b/src/thermophysicalModels/specie/specie/specieI.H
@@ -54,7 +54,6 @@ inline specie::specie
     const scalar molWeight
 )
 :
-    //name_(),
     nMoles_(nMoles),
     molWeight_(molWeight)
 {}
@@ -113,29 +112,29 @@ inline void specie::operator=(const specie& st)
 
 inline void specie::operator+=(const specie& st)
 {
-    scalar sumNmoles_ = max(nMoles_ + st.nMoles_, SMALL);
+    scalar sumNmoles = max(nMoles_ + st.nMoles_, SMALL);
 
     molWeight_ =
-        nMoles_/sumNmoles_*molWeight_
-      + st.nMoles_/sumNmoles_*st.molWeight_;
+        nMoles_/sumNmoles*molWeight_
+      + st.nMoles_/sumNmoles*st.molWeight_;
 
-    nMoles_ = sumNmoles_;
+    nMoles_ = sumNmoles;
 }
 
 
 inline void specie::operator-=(const specie& st)
 {
-    scalar diffnMoles_ = nMoles_ - st.nMoles_;
-    if (mag(diffnMoles_) < SMALL)
+    scalar diffnMoles = nMoles_ - st.nMoles_;
+    if (mag(diffnMoles) < SMALL)
     {
-        diffnMoles_ = SMALL;
+        diffnMoles = SMALL;
     }
 
     molWeight_ =
-        nMoles_/diffnMoles_*molWeight_
-      - st.nMoles_/diffnMoles_*st.molWeight_;
+        nMoles_/diffnMoles*molWeight_
+      - st.nMoles_/diffnMoles*st.molWeight_;
 
-    nMoles_ = diffnMoles_;
+    nMoles_ = diffnMoles;
 }
 
 
@@ -149,30 +148,30 @@ inline void specie::operator*=(const scalar s)
 
 inline specie operator+(const specie& st1, const specie& st2)
 {
-    scalar sumNmoles_ = max(st1.nMoles_ + st2.nMoles_, SMALL);
+    scalar sumNmoles = max(st1.nMoles_ + st2.nMoles_, SMALL);
 
     return specie
     (
-        sumNmoles_,
-        st1.nMoles_/sumNmoles_*st1.molWeight_
-      + st2.nMoles_/sumNmoles_*st2.molWeight_
+        sumNmoles,
+        st1.nMoles_/sumNmoles*st1.molWeight_
+      + st2.nMoles_/sumNmoles*st2.molWeight_
     );
 }
 
 
 inline specie operator-(const specie& st1, const specie& st2)
 {
-    scalar diffNmoles_ = st1.nMoles_ - st2.nMoles_;
-    if (mag(diffNmoles_) < SMALL)
+    scalar diffNmoles = st1.nMoles_ - st2.nMoles_;
+    if (mag(diffNmoles) < SMALL)
     {
-        diffNmoles_ = SMALL;
+        diffNmoles = SMALL;
     }
 
     return specie
     (
-        diffNmoles_,
-        st1.nMoles_/diffNmoles_*st1.molWeight_
-      - st2.nMoles_/diffNmoles_*st2.molWeight_
+        diffNmoles,
+        st1.nMoles_/diffNmoles*st1.molWeight_
+      - st2.nMoles_/diffNmoles*st2.molWeight_
     );
 }
 
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.C b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.C
index 497dc7c757bf83f5f466a6030efc32851bc91e0e..95504e7df41278aa27c0742df8d55b15194adf88 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.C
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.C
@@ -39,11 +39,21 @@ Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
     Hf_(readScalar(is)),
     Sf_(readScalar(is)),
     cpPolynomial_("cpPolynomial", is),
-    dhPolynomial_(cpPolynomial_.integrate()),
-    dsPolynomial_(cpPolynomial_.integrateMinus1())
+    hPolynomial_(),
+    sPolynomial_()
 {
-    // Offset dh poly so that it is relative to the enthalpy at Tstd
-    dhPolynomial_[0] -= dhPolynomial_.evaluate(specie::Tstd);
+    Hf_ *= this->W();
+    Sf_ *= this->W();
+    cpPolynomial_ *= this->W();
+
+    hPolynomial_ = cpPolynomial_.integrate();
+    sPolynomial_ = cpPolynomial_.integrateMinus1();
+
+    // Offset h poly so that it is relative to the enthalpy at Tstd
+    hPolynomial_[0] += Hf_ - hPolynomial_.evaluate(specie::Tstd);
+
+    // Offset s poly so that it is relative to the entropy at Tstd
+    sPolynomial_[0] += Sf_ - sPolynomial_.evaluate(specie::Tstd);
 }
 
 
@@ -57,15 +67,17 @@ Foam::Ostream& Foam::operator<<
 )
 {
     os  << static_cast<const EquationOfState&>(pt) << tab
-        << pt.Hf_ << tab
+        << pt.Hf_/pt.W() << tab
         << pt.Sf_ << tab
-        << pt.cpPolynomial_ << tab
-        << pt.dhPolynomial_ << tab
-        << pt.dsPolynomial;
+        << "cpPolynomial" << tab << pt.cpPolynomial_/pt.W();
 
     os.check
     (
-        "operator<<(Ostream& os, const hPolynomialThermo<EquationOfState>& pt)"
+        "operator<<"
+        "("
+            "Ostream&, "
+            "const hPolynomialThermo<EquationOfState, PolySize>&"
+        ")"
     );
 
     return os;
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
index a15d05adfe11822dae04268d061233c0cb96f28f..52b4482fa95d42c2986ede646224f2f5a1160f56 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
@@ -100,20 +100,22 @@ class hPolynomialThermo
 {
     // Private data
 
-        //- Heat of formation [J/kg]
+        //- Heat of formation
+        //  Note: input in [J/kg], but internally uses [J/kmol]
         scalar Hf_;
 
-        //- Standard entropy [J/(kg.K)]
+        //- Standard entropy
+        //  Note: input in [J/kg/K], but internally uses [J/kmol/K]
         scalar Sf_;
 
         //- Specific heat at constant pressure [J/(kg.K)]
         Polynomial<PolySize> cpPolynomial_;
 
         //- Enthalpy - derived from cp [J/kg] - relative to Tstd
-        typename Polynomial<PolySize>::intPolyType dhPolynomial_;
+        typename Polynomial<PolySize>::intPolyType hPolynomial_;
 
-        //- Entropy - derived from cp [J/(kg.K)]
-        Polynomial<PolySize> dsPolynomial_;
+        //- Entropy - derived from cp [J/(kg.K)] - relative to Tstd
+        Polynomial<PolySize> sPolynomial_;
 
 
     // Private member functions
@@ -125,8 +127,8 @@ class hPolynomialThermo
             const scalar Hf,
             const scalar Sf,
             const Polynomial<PolySize>& cpPoly,
-            const typename Polynomial<PolySize>::intPolyType& dhPoly,
-            const Polynomial<PolySize>& dsPoly
+            const typename Polynomial<PolySize>::intPolyType& hPoly,
+            const Polynomial<PolySize>& sPoly
         );
 
 
@@ -137,6 +139,9 @@ public:
         //- Construct from dictionary
         hPolynomialThermo(Istream& is);
 
+        //- Construct as copy
+        inline hPolynomialThermo(const hPolynomialThermo&);
+
         //- Construct as a named copy
         inline hPolynomialThermo(const word&, const hPolynomialThermo&);
 
@@ -161,8 +166,10 @@ public:
 
     // Member operators
 
+        inline hPolynomialThermo& operator=(const hPolynomialThermo&);
         inline void operator+=(const hPolynomialThermo&);
         inline void operator-=(const hPolynomialThermo&);
+        inline void operator*=(const scalar);
 
 
     // Friend operators
diff --git a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
index 21e2171ec060931aa865e224ad69e953b835d473..9714724955f0ecb8583548d6fcf1a4f013d508bb 100644
--- a/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
@@ -35,21 +35,36 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
     const scalar Hf,
     const scalar Sf,
     const Polynomial<PolySize>& cpPoly,
-    const typename Polynomial<PolySize>::intPolyType& dhPoly,
-    const Polynomial<PolySize>& dsPoly
+    const typename Polynomial<PolySize>::intPolyType& hPoly,
+    const Polynomial<PolySize>& sPoly
 )
 :
     EquationOfState(pt),
     Hf_(Hf),
     Sf_(Sf),
     cpPolynomial_(cpPoly),
-    dhPolynomial_(dhPoly),
-    dsPolynomial_(dsPoly)
+    hPolynomial_(hPoly),
+    sPolynomial_(sPoly)
 {}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+template<class EquationOfState, int PolySize>
+inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
+(
+    const hPolynomialThermo& pt
+)
+:
+    EquationOfState(pt),
+    Hf_(pt.Hf_),
+    Sf_(pt.Sf_),
+    cpPolynomial_(pt.cpPolynomial_),
+    hPolynomial_(pt.hPolynomial_),
+    sPolynomial_(pt.sPolynomial_)
+{}
+
+
 template<class EquationOfState, int PolySize>
 inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
 (
@@ -61,8 +76,8 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
     Hf_(pt.Hf_),
     Sf_(pt.Sf_),
     cpPolynomial_(pt.cpPolynomial_),
-    dhPolynomial_(pt.dhPolynomial_),
-    dsPolynomial_(pt.dsPolynomial_)
+    hPolynomial_(pt.hPolynomial_),
+    sPolynomial_(pt.sPolynomial_)
 {}
 
 
@@ -74,7 +89,7 @@ inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::cp
     const scalar T
 ) const
 {
-    return cpPolynomial_.evaluate(T)*this->W();
+    return cpPolynomial_.evaluate(T);
 }
 
 
@@ -84,7 +99,7 @@ inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::h
     const scalar T
 ) const
 {
-    return (dhPolynomial_.evaluate(T) + Hf_)*this->W();
+    return hPolynomial_.evaluate(T);
 }
 
 
@@ -94,7 +109,7 @@ inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::hs
     const scalar T
 ) const
 {
-    return dhPolynomial_.evaluate(T)*this->W();
+    return h(T) - hc();
 }
 
 
@@ -102,7 +117,7 @@ template<class EquationOfState, int PolySize>
 inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::hc()
 const
 {
-    return Hf_*this->W();
+    return Hf_;
 }
 
 
@@ -112,12 +127,31 @@ inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::s
     const scalar T
 ) const
 {
-    return (dsPolynomial_.evaluate(T) + Sf_)*this->W();
+    return sPolynomial_.evaluate(T);
 }
 
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
+template<class EquationOfState, int PolySize>
+inline Foam::hPolynomialThermo<EquationOfState, PolySize>&
+Foam::hPolynomialThermo<EquationOfState, PolySize>::operator=
+(
+    const hPolynomialThermo<EquationOfState, PolySize>& pt
+)
+{
+    EquationOfState::operator=(pt);
+
+    Hf_ = pt.Hf_;
+    Sf_ = pt.Sf_;
+    cpPolynomial_ = pt.cpPolynomial_;
+    hPolynomial_ = pt.hPolynomial_;
+    sPolynomial_ = pt.sPolynomial_;
+
+    return *this;
+}
+
+
 template<class EquationOfState, int PolySize>
 inline void Foam::hPolynomialThermo<EquationOfState, PolySize>::operator+=
 (
@@ -134,8 +168,8 @@ inline void Foam::hPolynomialThermo<EquationOfState, PolySize>::operator+=
     Hf_ = molr1*Hf_ + molr2*pt.Hf_;
     Sf_ = molr1*Sf_ + molr2*pt.Sf_;
     cpPolynomial_ = molr1*cpPolynomial_ + molr2*pt.cpPolynomial_;
-    dhPolynomial_ = molr1*dhPolynomial_ + molr2*pt.dhPolynomial_;
-    dsPolynomial_ = molr1*dsPolynomial_ + molr2*pt.dsPolynomial_;
+    hPolynomial_ = molr1*hPolynomial_ + molr2*pt.hPolynomial_;
+    sPolynomial_ = molr1*sPolynomial_ + molr2*pt.sPolynomial_;
 }
 
 
@@ -155,8 +189,18 @@ inline void Foam::hPolynomialThermo<EquationOfState, PolySize>::operator-=
     Hf_ = molr1*Hf_ - molr2*pt.Hf_;
     Sf_ = molr1*Sf_ - molr2*pt.Sf_;
     cpPolynomial_ = molr1*cpPolynomial_ - molr2*pt.cpPolynomial_;
-    dhPolynomial_ = molr1*dhPolynomial_ - molr2*pt.dhPolynomial_;
-    dsPolynomial_ = molr1*dsPolynomial_ - molr2*pt.dsPolynomial_;
+    hPolynomial_ = molr1*hPolynomial_ - molr2*pt.hPolynomial_;
+    sPolynomial_ = molr1*sPolynomial_ - molr2*pt.sPolynomial_;
+}
+
+
+template<class EquationOfState, int PolySize>
+inline void Foam::hPolynomialThermo<EquationOfState, PolySize>::operator*=
+(
+    const scalar s
+)
+{
+    EquationOfState::operator*=(s);
 }
 
 
@@ -169,22 +213,20 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize> Foam::operator+
     const hPolynomialThermo<EquationOfState, PolySize>& pt2
 )
 {
-    EquationOfState eofs
-    (
-        static_cast<const EquationOfState&>(pt1)
-      + static_cast<const EquationOfState&>(pt2)
-    );
+    EquationOfState eofs = pt1;
+    eofs += pt2;
 
     scalar molr1 = pt1.nMoles()/eofs.nMoles();
     scalar molr2 = pt2.nMoles()/eofs.nMoles();
+
     return hPolynomialThermo<EquationOfState, PolySize>
     (
         eofs,
         molr1*pt1.Hf_ + molr2*pt2.Hf_,
         molr1*pt1.Sf_ + molr2*pt2.Sf_,
         molr1*pt1.cpPolynomial_ + molr2*pt2.cpPolynomial_,
-        molr1*pt1.dhPolynomial_ + molr2*pt2.dhPolynomial_,
-        molr1*pt1.dsPolynomial_ + molr2*pt2.dsPolynomial_
+        molr1*pt1.hPolynomial_ + molr2*pt2.hPolynomial_,
+        molr1*pt1.sPolynomial_ + molr2*pt2.sPolynomial_
     );
 }
 
@@ -196,22 +238,20 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize> Foam::operator-
     const hPolynomialThermo<EquationOfState, PolySize>& pt2
 )
 {
-    EquationOfState eofs
-    (
-        static_cast<const EquationOfState&>(pt1)
-      - static_cast<const EquationOfState&>(pt2)
-    );
+    EquationOfState eofs = pt1;
+    eofs -= pt2;
 
     scalar molr1 = pt1.nMoles()/eofs.nMoles();
     scalar molr2 = pt2.nMoles()/eofs.nMoles();
+
     return hPolynomialThermo<EquationOfState, PolySize>
     (
         eofs,
         molr1*pt1.Hf_ - molr2*pt2.Hf_,
         molr1*pt1.Sf_ - molr2*pt2.Sf_,
         molr1*pt1.cpPolynomial_ - molr2*pt2.cpPolynomial_,
-        molr1*pt1.dhPolynomial_ - molr2*pt2.dhPolynomial_,
-        molr1*pt1.dsPolynomial_ - molr2*pt2.dsPolynomial_
+        molr1*pt1.hPolynomial_ - molr2*pt2.hPolynomial_,
+        molr1*pt1.sPolynomial_ - molr2*pt2.sPolynomial_
     );
 }
 
@@ -229,8 +269,8 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize> Foam::operator*
         pt.Hf_,
         pt.Sf_,
         pt.cpPolynomial_,
-        pt.dhPolynomial_,
-        pt.dsPolynomial_
+        pt.hPolynomial_,
+        pt.sPolynomial_
     );
 }
 
diff --git a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
index ea05b4c1234b0a55b154d8c38dd0bf0be7d19ea4..215ea7b2032eff381cd5512f7f239bb6e82484dc 100644
--- a/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
@@ -303,6 +303,7 @@ inline void Foam::specieThermo<thermo>::operator+=
     thermo::operator+=(st);
 }
 
+
 template<class thermo>
 inline void Foam::specieThermo<thermo>::operator-=
 (
@@ -312,6 +313,7 @@ inline void Foam::specieThermo<thermo>::operator-=
     thermo::operator-=(st);
 }
 
+
 template<class thermo>
 inline void Foam::specieThermo<thermo>::operator*=(const scalar s)
 {
diff --git a/src/thermophysicalModels/specie/transport/polynomial/polynomialTransport.C b/src/thermophysicalModels/specie/transport/polynomial/polynomialTransport.C
index 2744d4a97990960136447bbbda36ab1408c92def..0b6e912b0a09f738f683b22ca670c594f5cbb773 100644
--- a/src/thermophysicalModels/specie/transport/polynomial/polynomialTransport.C
+++ b/src/thermophysicalModels/specie/transport/polynomial/polynomialTransport.C
@@ -35,7 +35,10 @@ Foam::polynomialTransport<Thermo, PolySize>::polynomialTransport(Istream& is)
     Thermo(is),
     muPolynomial_("muPolynomial", is),
     kappaPolynomial_("kappaPolynomial", is)
-{}
+{
+    muPolynomial_ *= this->W();
+    kappaPolynomial_ *= this->W();
+}
 
 
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
@@ -47,7 +50,9 @@ Foam::Ostream& Foam::operator<<
     const polynomialTransport<Thermo, PolySize>& pt
 )
 {
-    os << static_cast<const Thermo&>(pt);
+    os  << static_cast<const Thermo&>(pt) << tab
+        << "muPolynomial" << tab << pt.muPolynomial_/pt.W() << tab
+        << "kappaPolynomial" << tab << pt.kappaPolynomial_/pt.W();
 
     os.check
     (
diff --git a/src/thermophysicalModels/specie/transport/polynomial/polynomialTransport.H b/src/thermophysicalModels/specie/transport/polynomial/polynomialTransport.H
index 3ab3879e64a203ac3915828d4bb667a1172b7443..0e8066f9d40033bf362a10f3475567bd3ca55153 100644
--- a/src/thermophysicalModels/specie/transport/polynomial/polynomialTransport.H
+++ b/src/thermophysicalModels/specie/transport/polynomial/polynomialTransport.H
@@ -96,9 +96,11 @@ class polynomialTransport
     // Private data
 
         //- Dynamic viscosity
+        //  Note: input in [Pa.s], but internally uses [Pa.s/kmol]
         Polynomial<PolySize> muPolynomial_;
 
         //- Thermal conductivity
+        //  Note: input in [W/m/K], but internally uses [W/m/K/kmol]
         Polynomial<PolySize> kappaPolynomial_;
 
 
@@ -117,6 +119,9 @@ public:
 
     // Constructors
 
+        //- Construct copy
+        inline polynomialTransport(const polynomialTransport&);
+
         //- Construct as named copy
         inline polynomialTransport(const word&, const polynomialTransport&);
 
@@ -148,6 +153,9 @@ public:
     // Member operators
 
         inline polynomialTransport& operator=(const polynomialTransport&);
+        inline void operator+=(const polynomialTransport&);
+        inline void operator-=(const polynomialTransport&);
+        inline void operator*=(const scalar);
 
 
     // Friend operators
diff --git a/src/thermophysicalModels/specie/transport/polynomial/polynomialTransportI.H b/src/thermophysicalModels/specie/transport/polynomial/polynomialTransportI.H
index 9e6e00e6e3c603146007caa7866e2e9aaf6751b1..ba891b0062d1516d6fed5bf9ca0995755be574a6 100644
--- a/src/thermophysicalModels/specie/transport/polynomial/polynomialTransportI.H
+++ b/src/thermophysicalModels/specie/transport/polynomial/polynomialTransportI.H
@@ -28,6 +28,18 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+template<class Thermo, int PolySize>
+inline Foam::polynomialTransport<Thermo, PolySize>::polynomialTransport
+(
+    const polynomialTransport& pt
+)
+:
+    Thermo(pt),
+    muPolynomial_(pt.muPolynomial_),
+    kappaPolynomial_(pt.kappaPolynomial_)
+{}
+
+
 template<class Thermo, int PolySize>
 inline Foam::polynomialTransport<Thermo, PolySize>::polynomialTransport
 (
@@ -85,7 +97,7 @@ inline Foam::scalar Foam::polynomialTransport<Thermo, PolySize>::mu
     const scalar T
 ) const
 {
-    return muPolynomial_.evaluate(T);
+    return muPolynomial_.evaluate(T)/this->W();
 }
 
 
@@ -95,7 +107,7 @@ inline Foam::scalar Foam::polynomialTransport<Thermo, PolySize>::kappa
     const scalar T
 ) const
 {
-    return kappaPolynomial_.evaluate(T);
+    return kappaPolynomial_.evaluate(T)/this->W();
 }
 
 
@@ -132,6 +144,52 @@ Foam::polynomialTransport<Thermo, PolySize>::operator=
 }
 
 
+template<class Thermo, int PolySize>
+inline void Foam::polynomialTransport<Thermo, PolySize>::operator+=
+(
+    const polynomialTransport<Thermo, PolySize>& pt
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Thermo::operator+=(pt);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = pt.nMoles()/this->nMoles();
+
+    muPolynomial_ = molr1*muPolynomial_ + molr2*pt.muPolynomial_;
+    kappaPolynomial_ = molr1*kappaPolynomial_ + molr2*pt.kappaPolynomial_;
+}
+
+
+template<class Thermo, int PolySize>
+inline void Foam::polynomialTransport<Thermo, PolySize>::operator-=
+(
+    const polynomialTransport<Thermo, PolySize>& pt
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Thermo::operator-=(pt);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = pt.nMoles()/this->nMoles();
+
+    muPolynomial_ = molr1*muPolynomial_ - molr2*pt.muPolynomial_;
+    kappaPolynomial_ = molr1*kappaPolynomial_ - molr2*pt.kappaPolynomial_;
+}
+
+
+template<class Thermo, int PolySize>
+inline void Foam::polynomialTransport<Thermo, PolySize>::operator*=
+(
+    const scalar s
+)
+{
+    Thermo::operator*=(s);
+}
+
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 template<class Thermo, int PolySize>
@@ -146,7 +204,6 @@ inline Foam::polynomialTransport<Thermo, PolySize> Foam::operator+
         static_cast<const Thermo&>(pt1) + static_cast<const Thermo&>(pt2)
     );
 
-
     scalar molr1 = pt1.nMoles()/t.nMoles();
     scalar molr2 = pt2.nMoles()/t.nMoles();
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
index c60d6628e5efea62d1e5c94cb821222104a65586..ed2cc7bf1385b6dd58282bf6f80a35eaf52bfd76 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
@@ -79,6 +79,7 @@ void omegaWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
     os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
     os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
     os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+    os.writeKeyword("beta1") << beta1_ << token::END_STATEMENT << nl;
 }
 
 
@@ -95,6 +96,7 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     Cmu_(0.09),
     kappa_(0.41),
     E_(9.8),
+    beta1_(0.075),
     yPlusLam_(calcYPlusLam(kappa_, E_))
 
 {
@@ -115,6 +117,7 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     Cmu_(ptf.Cmu_),
     kappa_(ptf.kappa_),
     E_(ptf.E_),
+    beta1_(ptf.beta1_),
     yPlusLam_(ptf.yPlusLam_)
 {
     checkType();
@@ -133,6 +136,7 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
     E_(dict.lookupOrDefault<scalar>("E", 9.8)),
+    beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
     yPlusLam_(calcYPlusLam(kappa_, E_))
 {
     checkType();
@@ -149,6 +153,7 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     Cmu_(owfpsf.Cmu_),
     kappa_(owfpsf.kappa_),
     E_(owfpsf.E_),
+    beta1_(owfpsf.beta1_),
     yPlusLam_(owfpsf.yPlusLam_)
 {
     checkType();
@@ -166,6 +171,7 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     Cmu_(owfpsf.Cmu_),
     kappa_(owfpsf.kappa_),
     E_(owfpsf.E_),
+    beta1_(owfpsf.beta1_),
     yPlusLam_(owfpsf.yPlusLam_)
 
 {
@@ -222,7 +228,11 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
             Cmu25*y[faceI]*sqrt(k[faceCellI])
            /(muw[faceI]/rhow[faceI]);
 
-        omega[faceCellI] = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
+        scalar omegaVis = 6.0*muw[faceI]/(rhow[faceI]*beta1_*sqr(y[faceI]));
+
+        scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
+
+        omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog));
 
         if (yPlus > yPlusLam_)
         {
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
index 2efcc1e191f977c0a6ac058e66d4f775e0eca42b..f57a11074e4e713b4df0cca6753db7ea6e0d4fb5 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
@@ -26,7 +26,23 @@ Class
     Foam::compressible::RASModels::omegaWallFunctionFvPatchScalarField
 
 Description
-    Replaces functionality in wallFunctionsI.H
+    Provides a wall function boundary condition/constraint on omega
+
+    Computed value is:
+
+        omega = sqrt(omega_vis^2 + omega_log^2)
+
+    where
+        omega_vis = omega in viscous region
+        omega_log = omega in logarithmic region
+
+    Model described by Eq.(15) of:
+    @verbatim
+        Menter, F., Esch, T.
+        "Elements of Industrial Heat Transfer Prediction"
+        16th Brazilian Congress of Mechanical Engineering (COBEM),
+        Nov. 2001
+    @endverbatim
 
 SourceFiles
     omegaWallFunctionFvPatchScalarField.C
@@ -71,6 +87,9 @@ protected:
         //- E coefficient
         scalar E_;
 
+        //- beta1 coefficient
+        scalar beta1_;
+
         //- Y+ at the edge of the laminar sublayer
         scalar yPlusLam_;
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
index b49454184672c7a8343059053502e4d5b743305c..084bbc2658a053efb8e84dcc385722a948b431ad 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
@@ -79,6 +79,7 @@ void omegaWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
     os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
     os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
     os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+    os.writeKeyword("beta1") << beta1_ << token::END_STATEMENT << nl;
 }
 
 
@@ -95,6 +96,7 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     Cmu_(0.09),
     kappa_(0.41),
     E_(9.8),
+    beta1_(0.075),
     yPlusLam_(calcYPlusLam(kappa_, E_))
 {
     checkType();
@@ -114,6 +116,7 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     Cmu_(ptf.Cmu_),
     kappa_(ptf.kappa_),
     E_(ptf.E_),
+    beta1_(ptf.beta1_),
     yPlusLam_(ptf.yPlusLam_)
 {
     checkType();
@@ -132,6 +135,7 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
     E_(dict.lookupOrDefault<scalar>("E", 9.8)),
+    beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
     yPlusLam_(calcYPlusLam(kappa_, E_))
 {
     checkType();
@@ -148,6 +152,7 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     Cmu_(owfpsf.Cmu_),
     kappa_(owfpsf.kappa_),
     E_(owfpsf.E_),
+    beta1_(owfpsf.beta1_),
     yPlusLam_(owfpsf.yPlusLam_)
 {
     checkType();
@@ -165,8 +170,8 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     Cmu_(owfpsf.Cmu_),
     kappa_(owfpsf.kappa_),
     E_(owfpsf.E_),
+    beta1_(owfpsf.beta1_),
     yPlusLam_(owfpsf.yPlusLam_)
-
 {
     checkType();
 }
@@ -217,7 +222,11 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
 
         scalar yPlus = Cmu25*y[faceI]*sqrt(k[faceCellI])/nuw[faceI];
 
-        omega[faceCellI] = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
+        scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI]));
+
+        scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
+
+        omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog));
 
         if (yPlus > yPlusLam_)
         {
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
index 1b3915b57594ef8ecbc3e1618488ef314483ab67..005f0331f26d5499b4da13ec2e93dfe2955ff89b 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
@@ -26,7 +26,23 @@ Class
     Foam::incompressible::RASModels::omegaWallFunctionFvPatchScalarField
 
 Description
-    Replaces functionality in wallFunctionsI.H
+    Provides a wall function boundary condition/constraint on omega
+
+    Computed value is:
+
+        omega = sqrt(omega_vis^2 + omega_log^2)
+
+    where
+        omega_vis = omega in viscous region
+        omega_log = omega in logarithmic region
+
+    Model described by Eq.(15) of:
+    @verbatim
+        Menter, F., Esch, T.
+        "Elements of Industrial Heat Transfer Prediction"
+        16th Brazilian Congress of Mechanical Engineering (COBEM),
+        Nov. 2001
+    @endverbatim
 
 SourceFiles
     omegaWallFunctionFvPatchScalarField.C
@@ -71,6 +87,9 @@ protected:
         //- E coefficient
         scalar E_;
 
+        //- beta1 coefficient
+        scalar beta1_;
+
         //- Y+ at the edge of the laminar sublayer
         scalar yPlusLam_;
 
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
index ef1151639172571123055850ae91da79c35fc4a9..f1b77471c9936684be2ea56083dd7954ccdeed3f 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
@@ -144,6 +144,8 @@ SingleMixtureFractionCoeffs
 
 LiquidEvaporationCoeffs
 {
+    enthalpyTransfer enthalpyDifference;
+
     activeLiquids
     (
         H2O
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties
index b1b89e17eb08dd6fd511ab06326e1c3126dfe826..eb71a37830768a2538b39d21f3676c22359698ca 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties
@@ -146,6 +146,8 @@ SinglePhaseMixtureCoeffs
 
 LiquidEvaporationCoeffs
 {
+    enthalpyTransfer enthalpyDifference;
+
     activeLiquids
     (
         H2O
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/thermo.incompressiblePoly b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/thermo.incompressiblePoly
index 6b8aeb6dcfafe9c616d27d0b6a37ff52021f3aa7..16729ac9d90b5f20abac77ae3668fd153cb33959 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/thermo.incompressiblePoly
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/thermo.incompressiblePoly
@@ -1,35 +1,34 @@
 (
-    // Nitrogen
-    N2 N2   1   28.0134
+N2 N2   1   28.0134
     rhoPolynomial
     (
-        3.88E+000
-       -1.64E-002
-        3.18E-005
-       -2.89E-008
-        9.90E-012
+        3.8936E+00
+       -1.6463E-02
+        3.2101E-05
+       -2.9174E-08
+        9.9889E-12
         0
         0
         0
     )
-    0.0                                      // Heat of formation [J/kg]
-    6840                                     // Standard entropy [J/kg/K]
+    0.0                                 // Heat of formation
+    0.0                                 // Standard entropy
     cpPolynomial
     (
-        2.75E+004
-        1.12E+001
-       -3.12E-002
-        4.45E-005
-       -1.92E-008
+        9.7908E+02
+        4.1787E-01
+       -1.1761E-03
+        1.6742E-06
+       -7.2559E-10
         0
         0
         0
     )
     muPolynomial
     (
-        1.54E-006
-        6.15E-008
-       -1.81E-011
+        1.5068E-06
+        6.1598E-08
+       -1.8188E-11
         0
         0
         0
@@ -38,47 +37,45 @@
     )
     kappaPolynomial
     (
-        3.15E-003
-        8.49E-005
-       -1.25E-008
+        3.1494E-03
+        8.4997E-05
+       -1.2621E-08
         0
         0
         0
         0
         0
     )
-
-    // Oxygen
-    O2 O2   1   31.9988
+O2 O2   1   31.9988
     rhoPolynomial
     (
-        4.43E+000
-       -1.87E-002
-        3.64E-005
-       -3.30E-008
-        1.13E-011
+        4.4475E+00
+       -1.8805E-02
+        3.6667E-05
+       -3.3323E-08
+        1.1410E-11
         0
         0
         0
     )
-    0.0                                  // Heat of formation [J/kg]
-    6408                                 // Standard entropy [J/kg/K]
+    0.0                                 // Heat of formation
+    0.0                                 // Standard entropy
     cpPolynomial
     (
-        2.67E+004
-        9.93E+000
-       -7.09E-003
-        1.45E-005
-       -9.13E-009
+        8.3484E+02
+        2.9297E-01
+       -1.4959E-04
+        3.4143E-07
+       -2.2786E-10
         0
         0
         0
     )
     muPolynomial
     (
-        1.54E-006
-        6.15E-008
-       -1.81E-011
+        1.5068E-06
+        6.1598E-08
+       -1.8188E-11
         0
         0
         0
@@ -87,25 +84,23 @@
     )
     kappaPolynomial
     (
-        2.09E-004
-        8.52E-005
-       -1.49E-008
+        1.6082E-04
+        8.5301E-05
+       -1.4998E-08
         0
         0
         0
         0
         0
     )
-
-    // Water vapour
-    H2O H2O 1   18.0153
+H2O H2O 1   18.0153
     rhoPolynomial
     (
-        2.49E+000
-       -1.05E-002
-        2.03E-005
-       -1.84E-008
-        6.28E-012
+        2.5039E+00
+       -1.0587E-02
+        2.0643E-05
+       -1.8761E-08
+        6.4237E-12
         0
         0
         0
@@ -114,20 +109,20 @@
     1.0482e04                                // Standard entropy [J/kg/K]
     cpPolynomial
     (
-        2.85E+004
-        2.63E+001
-       -4.63E-002
-        5.11E-005
-       -1.83E-008
+        1.5631E+03
+        1.6040E+00
+       -2.9334E-03
+        3.2168E-06
+       -1.1571E-09
         0
         0
         0
     )
     muPolynomial
     (
-        1.54E-006
-        6.15E-008
-       -1.81E-011
+        1.5068E-06
+        6.1598E-08
+       -1.8188E-11
         0
         0
         0
@@ -136,9 +131,56 @@
     )
     kappaPolynomial
     (
-        3.79E-003
-        1.53E-004
-       -1.22E-008
+        3.7972E-03
+        1.5336E-04
+       -1.1859E-08
+        0
+        0
+        0
+        0
+        0
+    )
+air air 1   28.85
+    rhoPolynomial
+    (
+        4.0097E+00
+       -1.6954E-02
+        3.3057E-05
+       -3.0042E-08
+        1.0286E-11
+        0
+        0
+        0
+    )
+    0.0
+    0.0
+    cpPolynomial
+    (
+        9.4876E+02
+        3.9171E-01
+       -9.5999E-04
+        1.3930E-06
+       -6.2029E-10
+        0
+        0
+        0
+    )
+    muPolynomial
+    (
+        1.5061E-06
+        6.1600E-08
+       -1.8190E-11
+        0
+        0
+        0
+        0
+        0
+    )
+    kappaPolynomial
+    (
+        2.5219E-03
+        8.5060E-05
+       -1.3120E-08
         0
         0
         0
@@ -146,3 +188,4 @@
         0
     )
 )
+
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/G b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/G
new file mode 100644
index 0000000000000000000000000000000000000000..6e8fb5d0e6c24f35752bcc2dd77f41aa6194dc53
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/G
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      G;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 0 -3 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    walls
+    {
+        type            MarshakRadiation;
+        T               T;
+        emissivity      1;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        type            MarshakRadiation;
+        T               T;
+        emissivity      1;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/H2O b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/H2O
new file mode 100644
index 0000000000000000000000000000000000000000..f49b0231f03bef98edc8f9eafcd17549f50a7a57
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/H2O
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      H2O;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+//        type            zeroGradient;
+        type            inletOutlet;
+        inletValue      uniform 0.01;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.01;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/T b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/T
new file mode 100644
index 0000000000000000000000000000000000000000..5d813907a8063eb389d767d6fe4b05ef0a949f98
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/T
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 473.0;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 473.0;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 473.0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/U b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..2ef228419242135ec9d36dded38f6dd5fedc9001
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/U
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    walls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform (0 1e-8 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/air b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/air
new file mode 100644
index 0000000000000000000000000000000000000000..8a461bd14756023f4e10ab13ffe7e942f541e10a
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/air
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.99;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+//        type            zeroGradient;
+        type            inletOutlet;
+        inletValue      uniform 0.99;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.99;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/p b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..9523758e82c0d8e662a5700317c4db542bd22a19
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/0/p
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 100000;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/README b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/README
new file mode 100644
index 0000000000000000000000000000000000000000..96640d8a25246f83a60db6f9eb7ff6a4b1dcd6f9
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/README
@@ -0,0 +1,5 @@
+Single water droplet
+
+One-way coupling with the Eulerian phase
+
+Droplet should approach the wet-bulb temperature
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/chemistryProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/chemistryProperties
new file mode 100644
index 0000000000000000000000000000000000000000..52f3cad99b4fd5c2d6d2a893f01868505005a65e
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/chemistryProperties
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      chemistryProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+rhoChemistryModel  ODEChemistryModel<icoPoly8ThermoPhysics>;
+
+chemistry       off;
+
+turbulentReaction off;
+
+chemistrySolver ode;
+
+initialChemicalTimeStep 1e-07;
+
+sequentialCoeffs
+{
+    cTauChem        0.001;
+    equilibriumRateLimiter off;
+}
+
+EulerImplicitCoeffs
+{
+    cTauChem        0.05;
+    equilibriumRateLimiter off;
+}
+
+odeCoeffs
+{
+    ODESolver       RK;
+    eps             0.05;
+    scale           1;
+}
+
+Cmix            Cmix [ 0 0 0 0 0 0 0 ] 1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/g b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..51944e7abd1f29f379b1643c1593cad94f7a6bfc
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 0 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/pointMassSourcesProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/pointMassSourcesProperties
new file mode 100644
index 0000000000000000000000000000000000000000..4a838a9e79178010ed63ef126e02e73fb6f4f6c7
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/pointMassSourcesProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      pointMassSourcesProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+active          false;
+
+pointSources
+(
+//none
+);
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/polyMesh/blockMeshDict b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..68bb6d93691a873021761f1c0b6fd8cc63431b67
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/polyMesh/blockMeshDict
@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.1;
+
+vertices        
+(
+    (0 0 0)
+    (1 0 0)
+    (1 1 0)
+    (0 1 0)
+    (0 0 0.1)
+    (1 0 0.1)
+    (1 1 0.1)
+    (0 1 0.1)
+);
+
+blocks          
+(
+    hex (0 1 2 3 4 5 6 7) (5 5 5) simpleGrading (1 1 1)
+);
+
+edges           
+(
+);
+
+patches         
+(
+    wall outlet 
+    (
+        (3 7 6 2)
+    )
+    wall inlet
+    (
+        (1 5 4 0)
+    )
+    wall walls
+    (
+        (0 4 7 3)
+        (2 6 5 1)
+    )
+    symmetryPlane back
+    (
+        (0 3 2 1)
+    )
+    symmetryPlane front
+    (
+        (4 5 6 7)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/polyMesh/boundary b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/polyMesh/boundary
new file mode 100755
index 0000000000000000000000000000000000000000..ad06bb3d15d8a42918e9e1ec3c2eb47177bac3d5
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/polyMesh/boundary
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+5
+(
+    outlet
+    {
+        type            patch;
+        nFaces          25;
+        startFace       300;
+    }
+    inlet
+    {
+        type            patch;
+        nFaces          25;
+        startFace       325;
+    }
+    walls
+    {
+        type            wall;
+        nFaces          50;
+        startFace       350;
+    }
+    back
+    {
+        type            symmetryPlane;
+        nFaces          25;
+        startFace       400;
+    }
+    front
+    {
+        type            symmetryPlane;
+        nFaces          25;
+        startFace       425;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/porousZones b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/porousZones
new file mode 100644
index 0000000000000000000000000000000000000000..2cbdbd636d454ff897d4dde2f881c1fda992995a
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/porousZones
@@ -0,0 +1,38 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      porousZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+(
+/*
+    porousRegion // name of cell zone
+    {
+        coordinateSystem
+        {
+            e1  (1 0 0);
+            e2  (0 1 0);
+        }
+
+        Darcy
+        {
+            d   d [0 -2 0 0 0 0 0] (500000 -1000 -1000);
+            f   f [0 -1 0 0 0 0 0] (0 0 0);
+        }
+    }
+*/
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/radiationProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/radiationProperties
new file mode 100644
index 0000000000000000000000000000000000000000..46e495c56ad200d8820b16acdd43ba33e7f7fc57
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/radiationProperties
@@ -0,0 +1,61 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation       off;
+
+radiationModel  none;
+
+solverFreq      10;
+
+P1Coeffs
+{
+    C               C [ 0 0 0 0 0 0 0 ] 0;
+}
+
+absorptionEmissionModel binaryAbsorptionEmission;
+
+binaryAbsorptionEmissionCoeffs
+{
+    model1
+    {
+        absorptionEmissionModel constantAbsorptionEmission;
+        constantAbsorptionEmissionCoeffs
+        {
+            a               a [ 0 -1 0 0 0 0 0 ] 0.5;
+            e               e [ 0 -1 0 0 0 0 0 ] 0.5;
+            E               E [ 1 -1 -3 0 0 0 0 ] 0;
+        }
+    }
+    model2
+    {
+        absorptionEmissionModel cloudAbsorptionEmission;
+        cloudAbsorptionEmissionCoeffs
+        {
+            cloudNames      ( coalCloud1 limestoneCloud1 );
+        }
+    }
+}
+
+scatterModel    cloudScatter;
+
+cloudScatterCoeffs
+{
+    cloudNames      ( coalCloud1 limestoneCloud1 );
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactingCloud1Positions b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactingCloud1Positions
new file mode 100644
index 0000000000000000000000000000000000000000..98646ad46e93bc7571169e63bc3e3a8ec2592b6e
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactingCloud1Positions
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       vectorField;
+    location    "constant";
+    object      reactingCloud1Positions;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+(
+(0.05 0.05 0.005)
+)
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactingCloud1Properties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
new file mode 100644
index 0000000000000000000000000000000000000000..30d317f52ed4fbc96876808a0d71711c81790cd5
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
@@ -0,0 +1,143 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      reactingCloud1Properties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+InjectionModel  ManualInjection;
+
+DragModel       SphereDrag;
+
+DispersionModel none;
+
+PatchInteractionModel StandardWallInteraction;
+
+HeatTransferModel RanzMarshall;
+
+CompositionModel SinglePhaseMixture;
+
+PhaseChangeModel LiquidEvaporation;
+
+PostProcessingModel none;
+
+radiation       off;
+
+coupled         false;
+
+cellValueSourceCorrection off;
+
+parcelTypeId    1;
+
+constantProperties
+{
+    rhoMin          rhoMin [ 1 -3 0 0 0 ] 1e-15;
+    TMin            TMin [ 0 0 0 1 0 ] 200;
+    pMin            pMin [ 1 -1 2 0 0 ] 1000;
+    rho0            rho0 [ 1 -3 0 0 0 ] 1000;
+    minParticleMass minParticleMass [ 1 0 0 0 0 ] 1e-15;
+    T0              T0 [ 0 0 0 1 0 ] 350;
+    cp0             cp0 [ 0 2 -2 -1 0 ] 4100;
+    epsilon0        epsilon0 [ 0 0 0 0 0 ] 1;
+    f0              f0 [ 0 0 0 0 0 ] 0.5;
+    Tvap            Tvap [ 0 0 0 1 0 ] 284;
+    Tbp             Tbp [ 0 0 0 1 0 ] 373;
+    Pr              Pr [ 0 0 0 0 0 ] 0.7;
+    constantVolume  false;
+}
+
+interpolationSchemes
+{
+    rho             cell;
+    U               cellPointFace;
+    mu              cell;
+    T               cell;
+    Cp              cell;
+    p               cell;
+}
+
+integrationSchemes
+{
+    U               Euler;
+    T               Euler;
+}
+
+particleForces
+{
+    gravity         on;
+    virtualMass     off;
+    Cvm             0.5;
+    pressureGradient off;
+    gradU           gradU;
+}
+
+ManualInjectionCoeffs
+{
+    massTotal       massTotal [ 1 0 0 0 0 ] 5.23599e-10; // 1 droplet of density 1000 kg/m3
+    parcelBasisType mass;
+    SOI             0;
+    positionsFile   "reactingCloud1Positions";
+    U0              (0 0 0);
+    parcelPDF
+    {
+        pdfType         uniform;
+        uniformPDF
+        {
+            minValue        100e-06;
+            maxValue        100e-06;
+        }
+    }
+}
+StandardWallInteractionCoeffs
+{
+    type            rebound;
+}
+
+RanzMarshallCoeffs
+{
+    BirdCorrection  off;
+}
+
+SinglePhaseMixtureCoeffs
+{
+    phases
+    (
+        liquid
+        {
+            H2O     1;
+        }
+    );
+}
+
+LiquidEvaporationCoeffs
+{
+    enthalpyTransfer enthalpyDifference; // latentHeat;
+
+    activeLiquids
+    (
+        H2O
+    );
+}
+
+PatchPostProcessingCoeffs
+{
+    maxStoredParcels 100;
+
+    patches
+    (
+        outlet
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactions b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactions
new file mode 100644
index 0000000000000000000000000000000000000000..f1c604b50761ea2ae27bb7e9aac6f8f0bced67bd
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/reactions
@@ -0,0 +1,9 @@
+species
+(
+    air
+    H2O
+);
+
+reactions
+(
+);
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/thermo.incompressiblePoly b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/thermo.incompressiblePoly
new file mode 100644
index 0000000000000000000000000000000000000000..16729ac9d90b5f20abac77ae3668fd153cb33959
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/thermo.incompressiblePoly
@@ -0,0 +1,191 @@
+(
+N2 N2   1   28.0134
+    rhoPolynomial
+    (
+        3.8936E+00
+       -1.6463E-02
+        3.2101E-05
+       -2.9174E-08
+        9.9889E-12
+        0
+        0
+        0
+    )
+    0.0                                 // Heat of formation
+    0.0                                 // Standard entropy
+    cpPolynomial
+    (
+        9.7908E+02
+        4.1787E-01
+       -1.1761E-03
+        1.6742E-06
+       -7.2559E-10
+        0
+        0
+        0
+    )
+    muPolynomial
+    (
+        1.5068E-06
+        6.1598E-08
+       -1.8188E-11
+        0
+        0
+        0
+        0
+        0
+    )
+    kappaPolynomial
+    (
+        3.1494E-03
+        8.4997E-05
+       -1.2621E-08
+        0
+        0
+        0
+        0
+        0
+    )
+O2 O2   1   31.9988
+    rhoPolynomial
+    (
+        4.4475E+00
+       -1.8805E-02
+        3.6667E-05
+       -3.3323E-08
+        1.1410E-11
+        0
+        0
+        0
+    )
+    0.0                                 // Heat of formation
+    0.0                                 // Standard entropy
+    cpPolynomial
+    (
+        8.3484E+02
+        2.9297E-01
+       -1.4959E-04
+        3.4143E-07
+       -2.2786E-10
+        0
+        0
+        0
+    )
+    muPolynomial
+    (
+        1.5068E-06
+        6.1598E-08
+       -1.8188E-11
+        0
+        0
+        0
+        0
+        0
+    )
+    kappaPolynomial
+    (
+        1.6082E-04
+        8.5301E-05
+       -1.4998E-08
+        0
+        0
+        0
+        0
+        0
+    )
+H2O H2O 1   18.0153
+    rhoPolynomial
+    (
+        2.5039E+00
+       -1.0587E-02
+        2.0643E-05
+       -1.8761E-08
+        6.4237E-12
+        0
+        0
+        0
+    )
+   -1.3423e07                                // Heat of formation [J/kg]
+    1.0482e04                                // Standard entropy [J/kg/K]
+    cpPolynomial
+    (
+        1.5631E+03
+        1.6040E+00
+       -2.9334E-03
+        3.2168E-06
+       -1.1571E-09
+        0
+        0
+        0
+    )
+    muPolynomial
+    (
+        1.5068E-06
+        6.1598E-08
+       -1.8188E-11
+        0
+        0
+        0
+        0
+        0
+    )
+    kappaPolynomial
+    (
+        3.7972E-03
+        1.5336E-04
+       -1.1859E-08
+        0
+        0
+        0
+        0
+        0
+    )
+air air 1   28.85
+    rhoPolynomial
+    (
+        4.0097E+00
+       -1.6954E-02
+        3.3057E-05
+       -3.0042E-08
+        1.0286E-11
+        0
+        0
+        0
+    )
+    0.0
+    0.0
+    cpPolynomial
+    (
+        9.4876E+02
+        3.9171E-01
+       -9.5999E-04
+        1.3930E-06
+       -6.2029E-10
+        0
+        0
+        0
+    )
+    muPolynomial
+    (
+        1.5061E-06
+        6.1600E-08
+       -1.8190E-11
+        0
+        0
+        0
+        0
+        0
+    )
+    kappaPolynomial
+    (
+        2.5219E-03
+        8.5060E-05
+       -1.3120E-08
+        0
+        0
+        0
+        0
+        0
+    )
+)
+
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/thermophysicalProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/thermophysicalProperties
new file mode 100644
index 0000000000000000000000000000000000000000..fbb4d94eb246537263e879359134454cf3e1a1ea
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/thermophysicalProperties
@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType      hRhoMixtureThermo<reactingMixture<icoPoly8ThermoPhysics>>;
+
+chemistryReader foamChemistryReader;
+
+foamChemistryFile "$FOAM_CASE/constant/reactions";
+
+foamChemistryThermoFile "$FOAM_CASE/constant/thermo.incompressiblePoly";
+
+liquidComponents
+(
+    H2O
+);
+
+H2O             H2O defaultCoeffs;
+
+solidComponents
+(
+);
+
+inertSpecie      air;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/turbulenceProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..693c0989a7a80636269ffddc6b231e1c21dff283
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/constant/turbulenceProperties
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/controlDict b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..64c9c37c8d8554bb1edaa163be42ed10b228cc8f
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/controlDict
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+startFoam       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.5;
+
+deltaT          1e-3;
+
+writeControl    adjustableRunTime;
+
+writeInterval   1e-1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  10;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           5;
+
+maxDeltaT       1e-3;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..99e07cad3dba5f11e04b385e3b8ad8dc982deacf
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSchemes
@@ -0,0 +1,78 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+//    grad(U)         cellLimited Gauss linear 1;
+//    grad(k)         cellLimited Gauss linear 1;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss upwind;
+    div(phid,p)     Gauss upwind;
+    div(phiU,p)     Gauss linear;
+    div(phi,h)      Gauss upwind;
+    div(phi,k)      Gauss upwind;
+    div(phi,epsilon) Gauss upwind;
+    div(phi,omega) Gauss upwind;
+//    div(U)          Gauss linear;
+    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div(phi,Yi_h)   Gauss upwind;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+/*
+    laplacian(muEff,U) Gauss linear corrected;
+    laplacian(mut,U) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DepsilonEff,omega) Gauss linear corrected;
+    laplacian(DREff,R) Gauss linear corrected;
+    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(alphaEff,h) Gauss linear corrected;
+*/
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..6d53c944858527baf4787780a1cbef7e63335122
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution
@@ -0,0 +1,83 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver           GAMG;
+
+        tolerance        0;
+        relTol           0.01;
+
+        smoother         DICGaussSeidel;
+        nPreSweeps       0;
+        nPostSweeps      2;
+
+        cacheAgglomeration true;
+
+        nCellsInCoarsestLevel 10;
+        agglomerator     faceAreaPair;
+        mergeLevels      1;
+
+        maxIter          100;
+    };
+    pFinal
+    {
+        $p
+        tolerance        1e-6;
+        relTol           0;
+    };
+    "(rho|G)"
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-06;
+        relTol          0;
+    };
+    "(Yi|h|)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-06;
+        relTol          0;
+    };
+    "(U|k|omega)"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-06;
+        relTol          0;
+    };
+}
+
+PISO
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    momentumPredictor yes;
+}
+
+additional
+{
+    dpdt            false;
+    eWork           true;
+    hWork           true;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/probesDict b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/probesDict
new file mode 100644
index 0000000000000000000000000000000000000000..ae04abedd6afb8ba59e3d3da4e67f7ade161f0e9
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/probesDict
@@ -0,0 +1,32 @@
+/*-------------------------------*- C++ -*---------------------------------*\
+|    =========                                                              |
+|    \\      /     OpenFOAM  1.4.1                                          |
+|     \\    /                                                               |
+|      \\  /       The Open Source CFD Toolbox                              |
+|       \\/                                        http://www.OpenFOAM.org  |
+\*-------------------------------------------------------------------------*/
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+    class           dictionary;
+    location        system;
+    object          probesDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// Fields to be probed. runTime modifiable!
+fields
+(
+    T H2O p kT
+);
+
+// Locations to be probed. runTime modifiable!
+probeLocations
+(
+    (0.005 0.0 0.0)
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/H2O b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/H2O
new file mode 100644
index 0000000000000000000000000000000000000000..3e46197b969abc2773cb84a71720294d192a5966
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/H2O
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      H2O;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 0.0;
+    }
+    inletSides
+    {
+        type            fixedValue;
+        value           uniform 0.01;
+    }
+    inletCentral
+    {
+        type            fixedValue;
+        value           uniform 0.01;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/T b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/T
new file mode 100644
index 0000000000000000000000000000000000000000..69f5653ca787665fe7d0ba17a5df972270ca4f80
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/T
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 473.0;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 473.0;
+    }
+    inletSides
+    {
+        type            fixedValue;
+        value           uniform 473.0;
+    }
+    inletCentral
+    {
+        type            fixedValue;
+        value           uniform 573.0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/U b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/U
new file mode 100644
index 0000000000000000000000000000000000000000..d692d79623a73644066733c7189cc454f7fe3491
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/U
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            flowRateInletVelocity;
+        flowRate        0.00379;
+        value           uniform (-0 10.82857143 -0);
+    }
+    inletSides
+    {
+        type            flowRateInletVelocity;
+        flowRate        0.00832;
+        value           uniform (-0 11.88571429 -0);
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    walls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/air b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/air
new file mode 100644
index 0000000000000000000000000000000000000000..5fbc4f5e7f973288ddede8ee4fa61fd762fdec18
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/air
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.99;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1.0;
+    }
+    inletSides
+    {
+        type            fixedValue;
+        value           uniform 0.99;
+    }
+    inletCentral
+    {
+        type            fixedValue;
+        value           uniform 0.99;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/alphat b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/alphat
new file mode 100644
index 0000000000000000000000000000000000000000..e75c946620c588519719d2361479e436daacbaa4
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/alphat
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    inletSides
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/k b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/k
new file mode 100644
index 0000000000000000000000000000000000000000..59479eef3aa4d4cd30dc50805e9d885a48e2e86e
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/k
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 3.75e-11;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.15;
+        value           uniform 3.75e-11;
+    }
+    inletSides
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.16;
+        value           uniform 3.75e-11;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    walls
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 3.75e-11;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/mut b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/mut
new file mode 100644
index 0000000000000000000000000000000000000000..f49efddd3304c690cfddd512da22a77dc5da7c3e
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/mut
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      mut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    inletSides
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            mutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/omega b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/omega
new file mode 100644
index 0000000000000000000000000000000000000000..f7c68ea1c4ff856d52a70ce0e2d6abf76c0b2751
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/omega
@@ -0,0 +1,61 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      omega;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 -1 0 0 0 0];
+
+internalField   uniform 0.0016;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            compressible::turbulentMixingLengthFrequencyInlet;
+        mixingLength    0.007;
+        k               k;
+        value           uniform 0.0016;
+    }
+    inletSides
+    {
+        type            compressible::turbulentMixingLengthFrequencyInlet;
+        mixingLength    0.007;
+        k               k;
+        value           uniform 0.0016;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    walls
+    {
+        type            compressible::omegaWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0.0016;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/p b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/p
new file mode 100644
index 0000000000000000000000000000000000000000..983de9aa2e99af13521e03fb504fd1a3f5d124df
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/p
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            zeroGradient;
+    }
+    inletSides
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 100000;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/H2O b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/H2O
new file mode 100644
index 0000000000000000000000000000000000000000..3e46197b969abc2773cb84a71720294d192a5966
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/H2O
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      H2O;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 0.0;
+    }
+    inletSides
+    {
+        type            fixedValue;
+        value           uniform 0.01;
+    }
+    inletCentral
+    {
+        type            fixedValue;
+        value           uniform 0.01;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/T b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/T
new file mode 100644
index 0000000000000000000000000000000000000000..69f5653ca787665fe7d0ba17a5df972270ca4f80
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/T
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 473.0;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 473.0;
+    }
+    inletSides
+    {
+        type            fixedValue;
+        value           uniform 473.0;
+    }
+    inletCentral
+    {
+        type            fixedValue;
+        value           uniform 573.0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/U b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..d692d79623a73644066733c7189cc454f7fe3491
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/U
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            flowRateInletVelocity;
+        flowRate        0.00379;
+        value           uniform (-0 10.82857143 -0);
+    }
+    inletSides
+    {
+        type            flowRateInletVelocity;
+        flowRate        0.00832;
+        value           uniform (-0 11.88571429 -0);
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    walls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/air b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/air
new file mode 100644
index 0000000000000000000000000000000000000000..5fbc4f5e7f973288ddede8ee4fa61fd762fdec18
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/air
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.99;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1.0;
+    }
+    inletSides
+    {
+        type            fixedValue;
+        value           uniform 0.99;
+    }
+    inletCentral
+    {
+        type            fixedValue;
+        value           uniform 0.99;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/alphat b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/alphat
new file mode 100644
index 0000000000000000000000000000000000000000..e75c946620c588519719d2361479e436daacbaa4
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/alphat
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    inletSides
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/k b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/k
new file mode 100644
index 0000000000000000000000000000000000000000..59479eef3aa4d4cd30dc50805e9d885a48e2e86e
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/k
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 3.75e-11;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.15;
+        value           uniform 3.75e-11;
+    }
+    inletSides
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.16;
+        value           uniform 3.75e-11;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    walls
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 3.75e-11;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/mut b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/mut
new file mode 100644
index 0000000000000000000000000000000000000000..f49efddd3304c690cfddd512da22a77dc5da7c3e
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/mut
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      mut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    inletSides
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            mutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/omega b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/omega
new file mode 100644
index 0000000000000000000000000000000000000000..f7c68ea1c4ff856d52a70ce0e2d6abf76c0b2751
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/omega
@@ -0,0 +1,61 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      omega;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 -1 0 0 0 0];
+
+internalField   uniform 0.0016;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            compressible::turbulentMixingLengthFrequencyInlet;
+        mixingLength    0.007;
+        k               k;
+        value           uniform 0.0016;
+    }
+    inletSides
+    {
+        type            compressible::turbulentMixingLengthFrequencyInlet;
+        mixingLength    0.007;
+        k               k;
+        value           uniform 0.0016;
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    walls
+    {
+        type            compressible::omegaWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0.0016;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/p b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..983de9aa2e99af13521e03fb504fd1a3f5d124df
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/p
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            zeroGradient;
+    }
+    inletSides
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 100000;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allclean b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..4843b685f3beeb2e4bdf69fb14eb483efd5fd773
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allclean
@@ -0,0 +1,12 @@
+#!/bin/sh
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+# remove old time folders
+rm -rf 0 *[1-9]* processor*
+
+# copy 0.org to 0
+cp -r 0.org 0
+
+cleanCase
+
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allrun b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..7aca4e2d349443694520edde716888a6059a92d2
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allrun
@@ -0,0 +1,9 @@
+#!/bin/sh
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# create mesh
+runApplication blockMesh
+
+# run the solver
+runApplication porousExplicitSourceReactingParcelFoam
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/RASProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/RASProperties
new file mode 100644
index 0000000000000000000000000000000000000000..02c52e85d0070b4258c314a2650e2fff5d442929
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/RASProperties
@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                     |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel            kOmegaSST; // kEpsilon;
+
+turbulence          on;
+
+printCoeffs         on;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/chemistryProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/chemistryProperties
new file mode 100644
index 0000000000000000000000000000000000000000..52f3cad99b4fd5c2d6d2a893f01868505005a65e
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/chemistryProperties
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      chemistryProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+rhoChemistryModel  ODEChemistryModel<icoPoly8ThermoPhysics>;
+
+chemistry       off;
+
+turbulentReaction off;
+
+chemistrySolver ode;
+
+initialChemicalTimeStep 1e-07;
+
+sequentialCoeffs
+{
+    cTauChem        0.001;
+    equilibriumRateLimiter off;
+}
+
+EulerImplicitCoeffs
+{
+    cTauChem        0.05;
+    equilibriumRateLimiter off;
+}
+
+odeCoeffs
+{
+    ODESolver       RK;
+    eps             0.05;
+    scale           1;
+}
+
+Cmix            Cmix [ 0 0 0 0 0 0 0 ] 1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/g b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..51944e7abd1f29f379b1643c1593cad94f7a6bfc
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 0 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/particleTrackProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/particleTrackProperties
new file mode 100644
index 0000000000000000000000000000000000000000..7c073e899cfc2f749e7af50c50dc4c4a6c3d0212
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/particleTrackProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      particleTrackProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+cloudName       reactingCloud1;
+
+sampleFrequency 1;
+
+maxPositions    1000000;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/pointMassSourcesProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/pointMassSourcesProperties
new file mode 100644
index 0000000000000000000000000000000000000000..2698204f7ce06b1851c2c6b3a49c4be1e0833ee1
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/pointMassSourcesProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      pointMassSourcesProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+active          false;
+
+pointSources
+(
+    //none
+);
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/blockMeshDict b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..8f1fd10943be05453af59a3795fd180c4a222904
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/blockMeshDict
@@ -0,0 +1,203 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                     |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+    location        "";
+    note            "Created Wed Jul  1 19:20:21 2009. Blocks = 8, cells = 9340, vertices = 36";
+    class           dictionary;
+    object          blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.001;
+
+vertices
+(
+    // front vertices
+    ( 0.00000e+00 -2.30000e+02  2.50000e+01) // v0 0
+    ( 0.00000e+00 -3.00000e+01  2.50000e+01) // v1 1
+    ( 0.00000e+00  0.00000e+00  2.50000e+01) // v2 2
+    ( 0.00000e+00  1.05000e+03  2.50000e+01) // v3 3
+    ( 9.00000e+00  1.05000e+03  2.50000e+01) // v4 4
+    ( 1.60000e+01  1.05000e+03  2.50000e+01) // v5 5
+    ( 2.50000e+01  1.05000e+03  2.50000e+01) // v6 6
+    ( 2.50000e+01  0.00000e+00  2.50000e+01) // v7 7
+    ( 2.50000e+01 -3.00000e+01  2.50000e+01) // v8 8
+    ( 2.50000e+01 -2.30000e+02  2.50000e+01) // v9 9
+    ( 1.80000e+01 -2.30000e+02  2.50000e+01) // v10 10
+    ( 1.80000e+01 -3.00000e+01  2.50000e+01) // v11 11
+    ( 1.60000e+01  0.00000e+00  2.50000e+01) // v12 12
+    ( 1.60000e+01 -2.30000e+02  2.50000e+01) // v13 13
+    ( 9.00000e+00 -2.30000e+02  2.50000e+01) // v14 14
+    ( 9.00000e+00  0.00000e+00  2.50000e+01) // v15 15
+    ( 7.00000e+00 -3.00000e+01  2.50000e+01) // v16 16
+    ( 7.00000e+00 -2.30000e+02  2.50000e+01) // v17 17
+
+    // back vertices
+    ( 0.00000e+00 -2.30000e+02 -2.50000e+01) // v0 18
+    ( 0.00000e+00 -3.00000e+01 -2.50000e+01) // v1 19
+    ( 0.00000e+00  0.00000e+00 -2.50000e+01) // v2 20
+    ( 0.00000e+00  1.05000e+03 -2.50000e+01) // v3 21
+    ( 9.00000e+00  1.05000e+03 -2.50000e+01) // v4 22
+    ( 1.60000e+01  1.05000e+03 -2.50000e+01) // v5 23
+    ( 2.50000e+01  1.05000e+03 -2.50000e+01) // v6 24
+    ( 2.50000e+01  0.00000e+00 -2.50000e+01) // v7 25
+    ( 2.50000e+01 -3.00000e+01 -2.50000e+01) // v8 26
+    ( 2.50000e+01 -2.30000e+02 -2.50000e+01) // v9 27
+    ( 1.80000e+01 -2.30000e+02 -2.50000e+01) // v10 28
+    ( 1.80000e+01 -3.00000e+01 -2.50000e+01) // v11 29
+    ( 1.60000e+01  0.00000e+00 -2.50000e+01) // v12 30
+    ( 1.60000e+01 -2.30000e+02 -2.50000e+01) // v13 31
+    ( 9.00000e+00 -2.30000e+02 -2.50000e+01) // v14 32
+    ( 9.00000e+00  0.00000e+00 -2.50000e+01) // v15 33
+    ( 7.00000e+00 -3.00000e+01 -2.50000e+01) // v16 34
+    ( 7.00000e+00 -2.30000e+02 -2.50000e+01) // v17 35
+);
+
+edges
+(
+);
+
+blocks
+(
+    // block 0
+    hex (0 1 16 17 18 19 34 35)
+    (67 10 10)
+    edgeGrading
+    (
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+    )
+
+    // block 1
+    hex (1 2 15 16 19 20 33 34)
+    (10 10 10)
+    edgeGrading
+    (
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+    )
+
+    // block 2
+    hex (2 3 4 15 20 21 22 33)
+    (234 10 10)
+    edgeGrading
+    (
+         4.00000e+00  4.00000e+00  4.00000e+00  4.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+    )
+
+    // block 3
+    hex (14 15 12 13 32 33 30 31)
+    (77 10 10)
+    edgeGrading
+    (
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+    )
+
+    // block 4
+    hex (15 4 5 12 33 22 23 30)
+    (234 10 10)
+    edgeGrading
+    (
+         4.00000e+00  4.00000e+00  4.00000e+00  4.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+    )
+
+    // block 5
+    hex (10 11 8 9 28 29 26 27)
+    (67 10 10)
+    edgeGrading
+    (
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+    )
+
+    // block 6
+    hex (11 12 7 8 29 30 25 26)
+    (11 10 10)
+    edgeGrading
+    (
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+    )
+
+    // block 7
+    hex (12 5 6 7 30 23 24 25)
+    (234 10 10)
+    edgeGrading
+    (
+         4.00000e+00  4.00000e+00  4.00000e+00  4.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+         1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
+    )
+
+);
+
+defaultPatch
+{
+    name walls;
+    type wall;
+}
+
+patches
+(
+    symmetryPlane back
+    (
+        (0 1 16 17)
+        (1 2 15 16)
+        (2 3 4 15)
+        (14 15 12 13)
+        (15 4 5 12)
+        (10 11 8 9)
+        (11 12 7 8)
+        (12 5 6 7)
+    )
+
+    symmetryPlane front
+    (
+        (18 19 34 35)
+        (19 20 33 34)
+        (20 21 22 33)
+        (32 33 30 31)
+        (33 22 23 30)
+        (28 29 26 27)
+        (29 30 25 26)
+        (30 23 24 25)
+    )
+
+    patch inletCentral
+    (
+        (13 14 32 31)
+    )
+
+    patch inletSides
+    (
+        (17 0 18 35)
+        (9 10 28 27)
+    )
+
+    patch outlet
+    (
+        (3 4 22 21)
+        (4 5 23 22)
+        (5 6 24 23)
+    )
+);
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/boundary b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/boundary
new file mode 100644
index 0000000000000000000000000000000000000000..f62411e7874ded628c4386d9b1d8c485b2fa36a8
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/boundary
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+6
+(
+    back
+    {
+        type            symmetryPlane;
+        nFaces          9340;
+        startFace       265900;
+    }
+    front
+    {
+        type            symmetryPlane;
+        nFaces          9340;
+        startFace       275240;
+    }
+    inletCentral
+    {
+        type            patch;
+        nFaces          100;
+        startFace       284580;
+    }
+    inletSides
+    {
+        type            patch;
+        nFaces          200;
+        startFace       284680;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          300;
+        startFace       284880;
+    }
+    walls
+    {
+        type            wall;
+        nFaces          9320;
+        startFace       285180;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/porousZones b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/porousZones
new file mode 100644
index 0000000000000000000000000000000000000000..2cbdbd636d454ff897d4dde2f881c1fda992995a
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/porousZones
@@ -0,0 +1,38 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      porousZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+(
+/*
+    porousRegion // name of cell zone
+    {
+        coordinateSystem
+        {
+            e1  (1 0 0);
+            e2  (0 1 0);
+        }
+
+        Darcy
+        {
+            d   d [0 -2 0 0 0 0 0] (500000 -1000 -1000);
+            f   f [0 -1 0 0 0 0 0] (0 0 0);
+        }
+    }
+*/
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/radiationProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/radiationProperties
new file mode 100644
index 0000000000000000000000000000000000000000..b05f5929254b27457f8d501abb44f7ebb6ce318e
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/radiationProperties
@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation       off;
+
+radiationModel  none;
+
+solverFreq      10;
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
new file mode 100644
index 0000000000000000000000000000000000000000..dcf5897492dfd0b752f6fe3842c40ce045c16fc6
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
@@ -0,0 +1,165 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      reactingCloud1Properties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+InjectionModel  PatchInjection;
+
+DragModel       SphereDrag;
+
+DispersionModel StochasticDispersionRAS;
+
+PatchInteractionModel StandardWallInteraction;
+
+HeatTransferModel RanzMarshall;
+
+CompositionModel SinglePhaseMixture;
+
+PhaseChangeModel LiquidEvaporation;
+
+PostProcessingModel PatchPostProcessing;
+
+radiation       off;
+
+coupled         true;
+
+cellValueSourceCorrection on;
+
+parcelTypeId    1;
+
+constantProperties
+{
+    rhoMin          rhoMin [ 1 -3 0 0 0 ] 1e-15;
+    TMin            TMin [ 0 0 0 1 0 ] 200;
+    pMin            pMin [ 1 -1 2 0 0 ] 1000;
+    rho0            rho0 [ 1 -3 0 0 0 ] 1000;
+    minParticleMass minParticleMass [ 1 0 0 0 0 ] 1e-15;
+    T0              T0 [ 0 0 0 1 0 ] 350;
+    cp0             cp0 [ 0 2 -2 -1 0 ] 4100;
+    epsilon0        epsilon0 [ 0 0 0 0 0 ] 1;
+    f0              f0 [ 0 0 0 0 0 ] 0.5;
+    Tvap            Tvap [ 0 0 0 1 0 ] 273;
+    Tbp             Tbp [ 0 0 0 1 0 ] 373;
+    Pr              Pr [ 0 0 0 0 0 ] 0.7;
+    constantVolume  false;
+}
+
+interpolationSchemes
+{
+    rho             cell;
+    U               cellPointFace;
+    mu              cell;
+    T               cell;
+    Cp              cell;
+    p               cell;
+}
+
+integrationSchemes
+{
+    U               Euler;
+    T               Analytical;
+}
+
+particleForces
+{
+    gravity         on;
+    virtualMass     off;
+    Cvm             0.5;
+    pressureGradient off;
+    gradU           gradU;
+}
+
+PatchInjectionCoeffs
+{
+    SOI             0.01;
+    massTotal       massTotal [ 1 0 0 0 0 ] 8;
+    parcelBasisType mass;
+    patchName       inletCentral;
+    duration        10000;
+    parcelsPerSecond 1e5;
+    U0              (40 0 0);
+    volumeFlowRate  constant 1;
+    parcelPDF
+    {
+        pdfType         general;
+        generalPDF
+        {
+            distribution
+            (
+                (10e-06      0.0025)
+                (15e-06      0.0528)
+                (20e-06      0.2795)
+                (25e-06      1.0918)
+                (30e-06      2.3988)
+                (35e-06      4.4227)
+                (40e-06      6.3888)
+                (45e-06      8.6721)
+                (50e-06      10.3153)
+                (55e-06      11.6259)
+                (60e-06      12.0030)
+                (65e-06      10.4175)
+                (70e-06      10.8427)
+                (75e-06      8.0016)
+                (80e-06      6.1333)
+                (85e-06      3.8827)
+                (90e-06      3.4688)
+            );
+        }
+    }
+}
+
+StandardWallInteractionCoeffs
+{
+    type            rebound;
+}
+
+RanzMarshallCoeffs
+{
+    BirdCorrection  on; // off;
+}
+
+SinglePhaseMixtureCoeffs
+{
+    phases
+    (
+        liquid
+        {
+            H2O     1;
+        }
+    );
+}
+
+LiquidEvaporationCoeffs
+{
+    enthalpyTransfer enthalpyDifference;
+
+    activeLiquids
+    (
+        H2O
+    );
+}
+
+PatchPostProcessingCoeffs
+{
+    maxStoredParcels 100;
+
+    patches
+    (
+        outlet
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactions b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactions
new file mode 100644
index 0000000000000000000000000000000000000000..f1c604b50761ea2ae27bb7e9aac6f8f0bced67bd
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactions
@@ -0,0 +1,9 @@
+species
+(
+    air
+    H2O
+);
+
+reactions
+(
+);
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/thermo.incompressiblePoly b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/thermo.incompressiblePoly
new file mode 100644
index 0000000000000000000000000000000000000000..16729ac9d90b5f20abac77ae3668fd153cb33959
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/thermo.incompressiblePoly
@@ -0,0 +1,191 @@
+(
+N2 N2   1   28.0134
+    rhoPolynomial
+    (
+        3.8936E+00
+       -1.6463E-02
+        3.2101E-05
+       -2.9174E-08
+        9.9889E-12
+        0
+        0
+        0
+    )
+    0.0                                 // Heat of formation
+    0.0                                 // Standard entropy
+    cpPolynomial
+    (
+        9.7908E+02
+        4.1787E-01
+       -1.1761E-03
+        1.6742E-06
+       -7.2559E-10
+        0
+        0
+        0
+    )
+    muPolynomial
+    (
+        1.5068E-06
+        6.1598E-08
+       -1.8188E-11
+        0
+        0
+        0
+        0
+        0
+    )
+    kappaPolynomial
+    (
+        3.1494E-03
+        8.4997E-05
+       -1.2621E-08
+        0
+        0
+        0
+        0
+        0
+    )
+O2 O2   1   31.9988
+    rhoPolynomial
+    (
+        4.4475E+00
+       -1.8805E-02
+        3.6667E-05
+       -3.3323E-08
+        1.1410E-11
+        0
+        0
+        0
+    )
+    0.0                                 // Heat of formation
+    0.0                                 // Standard entropy
+    cpPolynomial
+    (
+        8.3484E+02
+        2.9297E-01
+       -1.4959E-04
+        3.4143E-07
+       -2.2786E-10
+        0
+        0
+        0
+    )
+    muPolynomial
+    (
+        1.5068E-06
+        6.1598E-08
+       -1.8188E-11
+        0
+        0
+        0
+        0
+        0
+    )
+    kappaPolynomial
+    (
+        1.6082E-04
+        8.5301E-05
+       -1.4998E-08
+        0
+        0
+        0
+        0
+        0
+    )
+H2O H2O 1   18.0153
+    rhoPolynomial
+    (
+        2.5039E+00
+       -1.0587E-02
+        2.0643E-05
+       -1.8761E-08
+        6.4237E-12
+        0
+        0
+        0
+    )
+   -1.3423e07                                // Heat of formation [J/kg]
+    1.0482e04                                // Standard entropy [J/kg/K]
+    cpPolynomial
+    (
+        1.5631E+03
+        1.6040E+00
+       -2.9334E-03
+        3.2168E-06
+       -1.1571E-09
+        0
+        0
+        0
+    )
+    muPolynomial
+    (
+        1.5068E-06
+        6.1598E-08
+       -1.8188E-11
+        0
+        0
+        0
+        0
+        0
+    )
+    kappaPolynomial
+    (
+        3.7972E-03
+        1.5336E-04
+       -1.1859E-08
+        0
+        0
+        0
+        0
+        0
+    )
+air air 1   28.85
+    rhoPolynomial
+    (
+        4.0097E+00
+       -1.6954E-02
+        3.3057E-05
+       -3.0042E-08
+        1.0286E-11
+        0
+        0
+        0
+    )
+    0.0
+    0.0
+    cpPolynomial
+    (
+        9.4876E+02
+        3.9171E-01
+       -9.5999E-04
+        1.3930E-06
+       -6.2029E-10
+        0
+        0
+        0
+    )
+    muPolynomial
+    (
+        1.5061E-06
+        6.1600E-08
+       -1.8190E-11
+        0
+        0
+        0
+        0
+        0
+    )
+    kappaPolynomial
+    (
+        2.5219E-03
+        8.5060E-05
+       -1.3120E-08
+        0
+        0
+        0
+        0
+        0
+    )
+)
+
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/thermophysicalProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/thermophysicalProperties
new file mode 100644
index 0000000000000000000000000000000000000000..fbb4d94eb246537263e879359134454cf3e1a1ea
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/thermophysicalProperties
@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType      hRhoMixtureThermo<reactingMixture<icoPoly8ThermoPhysics>>;
+
+chemistryReader foamChemistryReader;
+
+foamChemistryFile "$FOAM_CASE/constant/reactions";
+
+foamChemistryThermoFile "$FOAM_CASE/constant/thermo.incompressiblePoly";
+
+liquidComponents
+(
+    H2O
+);
+
+H2O             H2O defaultCoeffs;
+
+solidComponents
+(
+);
+
+inertSpecie      air;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/turbulenceProperties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..d49f3cf71d02a4f2922cbece9fb87ac9508635b8
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/turbulenceProperties
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RASModel;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/controlDict b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..d31fcef49d8971ddc8034103da5d394d5af7e644
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/controlDict
@@ -0,0 +1,75 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+startFoam       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.5;
+
+deltaT          1e-05;
+
+writeControl    adjustableRunTime;
+
+writeInterval   1e-02;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  10;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           5;
+
+maxDeltaT       1e-03;
+
+functions
+{
+    faceSource1
+    {
+        type            faceSource;
+        functionObjectLibs ("libfieldFunctionObjects.so");
+        enabled         true;
+        outputControl   outputTime;
+        log             true;
+        valueOutput     true;
+        source          patch;
+        sourceName      outlet;
+        operation       weightedAverage;
+        weightField     phi;
+        fields
+        (
+            H2O
+            T
+        );
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..12543c5fd7c97d1402c4a9adef8a8f84fb0937eb
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSchemes
@@ -0,0 +1,67 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+    grad(U)         cellLimited Gauss linear 1;
+    grad(k)         cellLimited Gauss linear 1;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss upwind;
+    div(phid,p)     Gauss upwind;
+    div(phiU,p)     Gauss linear;
+    div(phi,h)      Gauss upwind;
+    div(phi,k)      Gauss upwind;
+    div(phi,epsilon) Gauss upwind;
+    div(phi,omega) Gauss upwind;
+    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div(phi,Yi_h)   Gauss upwind;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..1ccb5e9656f19e00c4eb61bba4b3451c0bcd7d10
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution
@@ -0,0 +1,91 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver           GAMG;
+        tolerance        0;
+        relTol           0.01;
+
+        smoother         DICGaussSeidel;
+        nPreSweeps       0;
+        nPostSweeps      2;
+
+        cacheAgglomeration true;
+
+        nCellsInCoarsestLevel 10;
+        agglomerator     faceAreaPair;
+        mergeLevels      1;
+
+        maxIter          50;
+    };
+    pFinal
+    {
+        $p
+        tolerance        1e-6;
+        relTol           0;
+    };
+    "(rho|G)"
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        maxIter         200;
+        tolerance       1e-06;
+        relTol          0;
+    };
+    "(Yi|h|)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        maxIter         200;
+        tolerance       1e-06;
+        relTol          0;
+    };
+    "(U|k|omega)"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        maxIter         200;
+        tolerance       1e-06;
+        relTol          0;
+    };
+}
+
+PISO
+{
+    nCorrectors     3;
+    nNonOrthogonalCorrectors 0;
+    momentumPredictor yes;
+}
+
+SIMPLE
+{
+    // used for potentialFoam initialisation
+    nNonOrthogonalCorrectors 10;
+}
+
+additional
+{
+    dpdt            false;
+    eWork           true;
+    hWork           true;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties
index 97619b95a78721d476e1ff9c26bfef1f8d3ce3bc..5a60bcc7155d15e5b09370185754f6bcb117a2a0 100644
--- a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties
@@ -120,6 +120,8 @@ SinglePhaseMixtureCoeffs
 
 LiquidEvaporationCoeffs
 {
+    enthalpyTransfer enthalpyDifference;
+
     activeLiquids
     (
         H2O