diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
index 0ddfa64588..000948514e 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
@@ -30,7 +30,7 @@ Description
 
     Sub-models include:
     - turbulence modelling, i.e. laminar, RAS or LES
-    - run-time selectable fvOptions, e.g. MRF, explicit porosity
+    - run-time selectable finitie volume options, e.g. MRF, explicit porosity
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/Allwmake b/applications/solvers/lagrangian/reactingParcelFoam/Allwmake
index ecd3f26015..f615df216c 100755
--- a/applications/solvers/lagrangian/reactingParcelFoam/Allwmake
+++ b/applications/solvers/lagrangian/reactingParcelFoam/Allwmake
@@ -3,6 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake
+wmake simpleReactingParcelFoam
 wmake LTSReactingParcelFoam
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C
index 991de15ad6..34aecb65b4 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C
@@ -28,7 +28,7 @@ Description
     Local time stepping (LTS) solver for steady, compressible, laminar or
     turbulent reacting and non-reacting flow with multiphase Lagrangian
     parcels and porous media, including run-time selectable finitite volume
-    options, e.g. fvOptions, constraints
+    options, e.g. sources, constraints
 
     Note: ddtPhiCorr not used here when porous zones are active
     - not well defined for porous calculations
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index 424c23b8b9..c4503cc101 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -22,12 +22,12 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    porousExplicitSourceReactingParcelFoam
+    reactingParcelFoam
 
 Description
     Transient PIMPLE solver for compressible, laminar or turbulent flow with
     reacting multiphase Lagrangian parcels, including run-time selectable
-    finite volume options, e.g. fvOptions, constraints
+    finite volume options, e.g. sources, constraints
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/EEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/EEqn.H
new file mode 100644
index 0000000000..d4d686e355
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/EEqn.H
@@ -0,0 +1,32 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        mvConvection->fvmDiv(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        parcels.Sh(he)
+      + radiation->Sh(thermo)
+      + combustion->Sh()
+      + fvOptions(rho, he)
+    );
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+    thermo.correct();
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/files b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/files
new file mode 100644
index 0000000000..4a202fcd4d
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/files
@@ -0,0 +1,3 @@
+simpleReactingParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/simpleReactingParcelFoam
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
new file mode 100644
index 0000000000..705cfc9433
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
@@ -0,0 +1,53 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lspecie \
+    -lfluidThermophysicalModels \
+    -lliquidProperties \
+    -lliquidMixtureProperties \
+    -lsolidProperties \
+    -lsolidMixtureProperties \
+    -lthermophysicalFunctions \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lradiationModels \
+    -lODE \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lcombustionModels \
+    -lfvOptions \
+    -lsampling
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H
new file mode 100644
index 0000000000..eabbb40454
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H
@@ -0,0 +1,17 @@
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho.dimensionedInternalField()*g
+      + parcels.SU(U)
+      + fvOptions(rho, U)
+    );
+
+    UEqn().relax();
+
+    fvOptions.constrain(UEqn());
+
+    solve(UEqn() == -fvc::grad(p));
+
+    fvOptions.correct(U);
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H
new file mode 100644
index 0000000000..cd0a45f0f0
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H
@@ -0,0 +1,53 @@
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+{
+    combustion->correct();
+    dQ = combustion->dQ();
+    label inertIndex = -1;
+    volScalarField Yt(0.0*Y[0]);
+
+    forAll(Y, i)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YEqn
+            (
+                mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+             ==
+                parcels.SYi(i, Yi)
+              + combustion->R(Yi)
+              + fvOptions(rho, Yi)
+            );
+
+            YEqn.relax();
+
+            fvOptions.constrain(YEqn);
+
+            YEqn.solve(mesh.solver("Yi"));
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H
new file mode 100644
index 0000000000..954b74e069
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createClouds.H
@@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicReactingMultiphaseCloud parcels
+(
+    "reactingCloud1",
+    rho,
+    U,
+    g,
+    slgThermo
+);
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H
new file mode 100644
index 0000000000..d6df24cb48
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H
@@ -0,0 +1,98 @@
+    Info<< "Creating combustion model\n" << endl;
+
+    autoPtr<combustionModels::rhoCombustionModel> combustion
+    (
+        combustionModels::rhoCombustionModel::New(mesh)
+    );
+
+    rhoReactionThermo& thermo = combustion->thermo();
+    thermo.validate(args.executable(), "h", "e");
+
+    SLGThermo slgThermo(mesh, thermo);
+
+    basicMultiComponentMixture& composition = thermo.composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
+    const word inertSpecie(thermo.lookup("inertSpecie"));
+
+    if (!composition.contains(inertSpecie))
+    {
+        FatalErrorIn(args.executable())
+            << "Specified inert specie '" << inertSpecie << "' not found in "
+            << "species list. Available species:" << composition.species()
+            << exit(FatalError);
+    }
+
+    volScalarField& p = thermo.p();
+    const volScalarField& T = thermo.T();
+    const volScalarField& psi = thermo.psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+    dimensionedScalar rhoMax(simple.dict().lookup("rhoMax"));
+    dimensionedScalar rhoMin(simple.dict().lookup("rhoMin"));
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    // Set the turbulence into the combustion model
+    combustion->setTurbulence(turbulence());
+
+    Info<< "Creating multi-variate interpolation scheme\n" << endl;
+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+    forAll(Y, i)
+    {
+        fields.add(Y[i]);
+    }
+    fields.add(thermo.he());
+
+    volScalarField dQ
+    (
+        IOobject
+        (
+            "dQ",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    );
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
new file mode 100644
index 0000000000..0d1aa2e238
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
@@ -0,0 +1,59 @@
+{
+    rho = thermo.rho();
+
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p;
+
+    volScalarField rAU(1.0/UEqn().A());
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
+
+    UEqn.clear();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    fvOptions.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rho*rAU, p)
+         ==
+            parcels.Srho()
+          + fvOptions(psi, p, rho.name())
+        );
+
+        fvOptions.constrain(pEqn);
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+
+    p.relax();
+
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
+
+    #include "compressibleContinuityErrs.H"
+
+    U = HbyA - rAU*fvc::grad(p);
+    U.correctBoundaryConditions();
+    fvOptions.correct(U);
+
+    rho = thermo.rho();
+    rho = max(rho, rhoMin);
+    rho = min(rho, rhoMax);
+
+    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
+}
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C
new file mode 100644
index 0000000000..6620d2af52
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C
@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    simpleReactingParcelFoam
+
+Description
+    Steady state SIMPLE solver for compressible, laminar or turbulent flow with
+    reacting multiphase Lagrangian parcels, including run-time selectable
+    finite volume options, e.g. sources, constraints
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "turbulenceModel.H"
+#include "basicReactingMultiphaseCloud.H"
+#include "rhoCombustionModel.H"
+#include "radiationModel.H"
+#include "IOporosityModelList.H"
+#include "fvIOoptionList.H"
+#include "SLGThermo.H"
+#include "simpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+
+    simpleControl simple(mesh);
+
+    #include "createFields.H"
+    #include "createRadiationModel.H"
+    #include "createClouds.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        parcels.evolve();
+
+        // --- Pressure-velocity SIMPLE corrector loop
+        {
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "EEqn.H"
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/Circulator/Test-Circulator.C b/applications/test/Circulator/Test-Circulator.C
index 19a21f6b79..2d8ecd3765 100644
--- a/applications/test/Circulator/Test-Circulator.C
+++ b/applications/test/Circulator/Test-Circulator.C
@@ -38,94 +38,6 @@ Description
 using namespace Foam;
 
 
-// return
-//   0: no match
-//  +1: identical
-//  -1: same face, but different orientation
-label compare(const face& a, const face& b)
-{
-    // Basic rule: we assume that the sequence of labels in each list
-    // will be circular in the same order (but not necessarily in the
-    // same direction or from the same starting point).
-
-    // Trivial reject: faces are different size
-    label sizeA = a.size();
-    label sizeB = b.size();
-
-    if (sizeA != sizeB || sizeA == 0)
-    {
-        return 0;
-    }
-
-    const_circulator<face> aCirc(a);
-    const_circulator<face> bCirc(b);
-
-    // Rotate face b until its element matches the starting element of face a.
-    do
-    {
-        if (aCirc() == bCirc())
-        {
-            // Set bCirc fulcrum to its iterator and increment the iterators
-            bCirc.setFulcrumToIterator();
-            ++aCirc;
-            ++bCirc;
-
-            break;
-        }
-    } while (bCirc.circulate(CirculatorBase::CLOCKWISE));
-
-    // Look forwards around the faces for a match
-    do
-    {
-        if (aCirc() != bCirc())
-        {
-            break;
-        }
-    }
-    while
-    (
-        aCirc.circulate(CirculatorBase::CLOCKWISE),
-        bCirc.circulate(CirculatorBase::CLOCKWISE)
-    );
-
-    // If the circulator has stopped then faces a and b matched.
-    if (!aCirc.circulate())
-    {
-        return 1;
-    }
-    else
-    {
-        // Reset the circulators back to their fulcrum
-        aCirc.setIteratorToFulcrum();
-        bCirc.setIteratorToFulcrum();
-        ++aCirc;
-        --bCirc;
-    }
-
-    // Look backwards around the faces for a match
-    do
-    {
-        if (aCirc() != bCirc())
-        {
-            break;
-        }
-    }
-    while
-    (
-        aCirc.circulate(CirculatorBase::CLOCKWISE),
-        bCirc.circulate(CirculatorBase::ANTICLOCKWISE)
-    );
-
-    // If the circulator has stopped then faces a and b matched.
-    if (!aCirc.circulate())
-    {
-        return -1;
-    }
-
-    return 0;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:
 
@@ -184,40 +96,44 @@ int main(int argc, char *argv[])
 
     Info<< nl << nl << "Compare two faces: " << endl;
     face a(identity(5));
-    Info<< "Compare " << a << " and " << a << " Match = " << compare(a, a)
+    Info<< "Compare " << a << " and " << a << " Match = " << face::compare(a, a)
         << endl;
 
     face b(reverseList(a));
-    Info<< "Compare " << a << " and " << b << " Match = " << compare(a, b)
+    Info<< "Compare " << a << " and " << b << " Match = " << face::compare(a, b)
         << endl;
 
     face c(a);
     c[4] = 3;
-    Info<< "Compare " << a << " and " << c << " Match = " << compare(a, c)
+    Info<< "Compare " << a << " and " << c << " Match = " << face::compare(a, c)
         << endl;
 
     face d(rotateList(a, 2));
-    Info<< "Compare " << a << " and " << d << " Match = " << compare(a, d)
+    Info<< "Compare " << a << " and " << d << " Match = " << face::compare(a, d)
         << endl;
 
     face g(labelList(5, 1));
     face h(g);
-    Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+    Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h)
         << endl;
 
     g[0] = 2;
     h[3] = 2;
-    Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+    Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h)
         << endl;
 
     g[4] = 3;
     h[4] = 3;
-    Info<< "Compare " << g << " and " << h << " Match = " << compare(g, h)
+    Info<< "Compare " << g << " and " << h << " Match = " << face::compare(g, h)
         << endl;
 
     face face1(identity(1));
     Info<< "Compare " << face1 << " and " << face1
-        << " Match = " << compare(face1, face1) << endl;
+        << " Match = " << face::compare(face1, face1) << endl;
+
+    face face2(identity(1)+1);
+    Info<< "Compare " << face1 << " and " << face2
+        << " Match = " << face::compare(face1, face2) << endl;
 
     Info<< nl << nl << "Zero face" << nl << endl;
 
diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseDict b/applications/utilities/mesh/advanced/collapseEdges/collapseDict
index 60f4cb2f43..340c33049f 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/collapseDict
+++ b/applications/utilities/mesh/advanced/collapseEdges/collapseDict
@@ -17,7 +17,7 @@ FoamFile
 collapseEdgesCoeffs
 {
     // Edges shorter than this absolute value will be merged
-    minimumEdgeLength   1e-8;
+    minimumEdgeLength   1e-6;
 
     // The maximum angle between two edges that share a point attached to
     // no other edges
@@ -25,7 +25,7 @@ collapseEdgesCoeffs
 
     // The amount that minimumEdgeLength will be reduced by for each
     // edge if that edge's collapse generates a poor quality face
-    reductionFactor                         0.5;
+    reductionFactor     0.5;
 }
 
 
@@ -67,14 +67,17 @@ meshQualityCoeffs
 {
     // Name of the dictionary that has the mesh quality coefficients used
     // by motionSmoother::checkMesh
-    meshQualityCoeffDict    meshQualityDict;
+    #include                    "meshQualityDict";
+
+    // Maximum number of smoothing iterations for the reductionFactors
+    maximumSmoothingIterations  2;
 
     // Maximum number of outer iterations is mesh quality checking is enabled
-    maximumIterations       10;
+    maximumIterations           10;
 
     // Maximum number of iterations deletion of a point can cause a bad face
     // to be constructed before it is forced to not be deleted
-    maxPointErrorCount      5;
+    maxPointErrorCount          5;
 }
 
 
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options
index f1db60d53f..fe1283a73a 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options
@@ -1,5 +1,4 @@
 EXE_INC = \
-    /* -DFULLDEBUG -g -O0 */ \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
index 7979fa966a..1b83599b31 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
@@ -26,12 +26,6 @@ License
 #include "patchToPoly2DMesh.H"
 #include "PatchTools.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::patchToPoly2DMesh::flipFaceOrder()
@@ -275,12 +269,9 @@ void Foam::patchToPoly2DMesh::createPolyMeshComponents()
     addPatchFacesToOwner();
 }
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-//- Construct from a primitivePatch
 Foam::patchToPoly2DMesh::patchToPoly2DMesh
 (
     const MeshedSurface<face>& patch,
@@ -321,13 +312,5 @@ void Foam::patchToPoly2DMesh::createMesh()
     }
 }
 
-// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
-
 
 // ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
index ea27dc6f53..b7c81e68a2 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
@@ -51,6 +51,7 @@ class patchToPoly2DMesh
 {
     // Private data
 
+        // Reference to the meshed surface
         const MeshedSurface<face>& patch_;
 
         const wordList& patchNames_;
@@ -62,11 +63,16 @@ class patchToPoly2DMesh
         const EdgeMap<label>& mapEdgesRegion_;
 
         pointField points_;
+
         faceList faces_;
+
         labelList owner_;
+
         labelList neighbour_;
 
-        //- Description of data_
+
+    // Private Member Functions
+
         void flipFaceOrder();
 
         void createNeighbours();
@@ -79,9 +85,6 @@ class patchToPoly2DMesh
 
         void createPolyMeshComponents();
 
-
-    // Private Member Functions
-
         //- Disallow default bitwise copy construct
         patchToPoly2DMesh(const patchToPoly2DMesh&);
 
@@ -110,47 +113,47 @@ public:
     // Member Functions
 
         // Access
-        pointField& points()
-        {
-            return points_;
-        }
 
-        faceList& faces()
-        {
-            return faces_;
-        }
+            pointField& points()
+            {
+                return points_;
+            }
 
-        labelList& owner()
-        {
-            return owner_;
-        }
+            faceList& faces()
+            {
+                return faces_;
+            }
 
-        labelList& neighbour()
-        {
-            return neighbour_;
-        }
+            labelList& owner()
+            {
+                return owner_;
+            }
 
-        const wordList& patchNames() const
-        {
-            return patchNames_;
-        }
+            labelList& neighbour()
+            {
+                return neighbour_;
+            }
 
-        const labelList& patchSizes() const
-        {
-            return patchSizes_;
-        }
+            const wordList& patchNames() const
+            {
+                return patchNames_;
+            }
+
+            const labelList& patchSizes() const
+            {
+                return patchSizes_;
+            }
+
+            const labelList& patchStarts() const
+            {
+                return patchStarts_;
+            }
 
-        const labelList& patchStarts() const
-        {
-            return patchStarts_;
-        }
-        // Check
 
         // Edit
-        void createMesh();
-
-        // Write
 
+            //- Create the mesh
+            void createMesh();
 };
 
 
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
index 088d0594f8..0b04d91b72 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
@@ -317,12 +317,12 @@ int main(int argc, char *argv[])
     }
 
     // Take over refinement levels and write to new time directory.
-    Pout<< "\nWriting extruded mesh to time = " << runTimeExtruded.timeName()
+    Info<< "\nWriting extruded mesh to time = " << runTimeExtruded.timeName()
         << nl << endl;
 
     mesh().write();
 
-    Pout<< "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return 0;
 }
diff --git a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H
index f5186c3a25..7aa4459e89 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H
+++ b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/CGALPolyhedronRings.H
@@ -2,31 +2,31 @@
 #define CGAL_PSURF_RINGS_H_
 
 // This file adapted from
-//     CGAL-3.7/examples/Jet_fitting_3//PolyhedralSurf_rings.h
-// Licensed under CGAL-3.7/LICENSE.FREE_USE
+//     CGAL-4.0/examples/Jet_fitting_3/PolyhedralSurf_rings.h
+// Licensed under CGAL-4.0/LICENSE.FREE_USE
 
 // Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007
-// Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie
-// Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France),
-// Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute
-// Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University
-// (Israel).  All rights reserved.
-
-// Permission is hereby granted, free of charge, to any person obtaining a copy
-// of this software and associated documentation files (the "Software"), to deal
-// in the Software without restriction, including without limitation the rights
-// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the Software is
-// furnished to do so, subject to the following conditions:
-
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
-// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-// SOFTWARE.
-
+// Utrecht University (The Netherlands),
+// ETH Zurich (Switzerland),
+// INRIA Sophia-Antipolis (France),
+// Max-Planck-Institute Saarbruecken (Germany),
+// and Tel-Aviv University (Israel).  All rights reserved.
+//
+// Permission is hereby granted, free of charge, to any person obtaining a
+// copy of this software and associated documentation files (the
+// "Software"), to deal in the Software without restriction, including
+// without limitation the rights to use, copy, modify, merge, publish,
+// distribute, sublicense, and/or sell copies of the Software, and to
+// permit persons to whom the Software is furnished to do so, subject to
+// the following conditions:
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
 #include <cassert>
 
diff --git a/etc/codeTemplates/dynamicCode/functionObjectTemplate.H b/etc/codeTemplates/dynamicCode/functionObjectTemplate.H
index bd07fff034..211c0868ec 100644
--- a/etc/codeTemplates/dynamicCode/functionObjectTemplate.H
+++ b/etc/codeTemplates/dynamicCode/functionObjectTemplate.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,6 +48,7 @@ namespace Foam
 // Forward declaration of classes
 class objectRegistry;
 class dictionary;
+class polyMesh;
 class mapPolyMesh;
 class fvMesh;
 
@@ -131,7 +132,7 @@ public:
         {}
 
         //- Update for changes of mesh
-        virtual void movePoints(const pointField&)
+        virtual void movePoints(const polyMesh&)
         {}
 };
 
diff --git a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H
index a50d5e3995..7309af83e1 100644
--- a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H
+++ b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H
@@ -25,7 +25,8 @@ Class
     Foam::dynamicIndexedOctree
 
 Description
-    Non-pointer based hierarchical recursive searching
+    Non-pointer based hierarchical recursive searching.
+    Storage is dynamic, so elements can be deleted.
 
 SourceFiles
     dynamicIndexedOctree.C
@@ -451,6 +452,7 @@ public:
             );
         }
 
+
     // Member Functions
 
         // Access
@@ -656,6 +658,7 @@ public:
 
             label removeIndex(const label nodIndex, const label index);
 
+
         // Write
 
             //- Print tree. Either print all indices (printContent = true) or
@@ -671,6 +674,7 @@ public:
 
             void writeTreeInfo() const;
 
+
     // IOstream Operators
 
         friend Ostream& operator<< <Type>
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C
index 1d7ee06c2e..21576535e9 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.C
@@ -318,6 +318,17 @@ int Foam::face::compare(const face& a, const face& b)
     {
         return 0;
     }
+    else if (sizeA == 1)
+    {
+        if (a[0] == b[0])
+        {
+            return 1;
+        }
+        else
+        {
+            return 0;
+        }
+    }
 
     const_circulator<face> aCirc(a);
     const_circulator<face> bCirc(b);
@@ -337,7 +348,7 @@ int Foam::face::compare(const face& a, const face& b)
     } while (bCirc.circulate(CirculatorBase::CLOCKWISE));
 
     // If the circulator has stopped then faces a and b do not share a matching
-    // point
+    // point. Doesn't work on matching, single element face.
     if (!bCirc.circulate())
     {
         return 0;
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
index ef697325b6..df95c96818 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -538,6 +538,9 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
     label cI
 )
 {
+    static label nWarnings = 0;
+    static const label maxWarnings = 100;
+
     const faceList& pFaces = mesh.faces();
     const labelList& pOwner = mesh.faceOwner();
 
@@ -553,19 +556,27 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
 
     if (tetBasePtI == -1)
     {
-        WarningIn
-        (
-            "Foam::List<Foam::tetIndices> "
-            "Foam::polyMeshTetDecomposition::faceTetIndices"
-            "("
-                "const polyMesh&, "
-                "label, "
-                "label"
-            ")"
-        )
-            << "No base point for face " << fI << ", " << f
-            << ", produces a valid tet decomposition."
-            << endl;
+        if (nWarnings < maxWarnings)
+        {
+            WarningIn
+            (
+                "Foam::List<Foam::tetIndices> "
+                "Foam::polyMeshTetDecomposition::faceTetIndices"
+                "("
+                    "const polyMesh&, "
+                    "label, "
+                    "label"
+                ")"
+            )   << "No base point for face " << fI << ", " << f
+                << ", produces a valid tet decomposition."
+                << endl;
+            nWarnings++;
+        }
+        if (nWarnings == maxWarnings)
+        {
+            Warning<< "Suppressing any further warnings." << endl;
+            nWarnings++;
+        }
 
         tetBasePtI = 0;
     }
diff --git a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
index c44f646025..9b184f994e 100644
--- a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
+++ b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
@@ -35,6 +35,12 @@ License
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
+namespace Foam
+{
+defineTypeNameAndDebug(polyMeshFilter, 0);
+}
+
+
 Foam::autoPtr<Foam::fvMesh> Foam::polyMeshFilter::copyMesh(const fvMesh& mesh)
 {
     polyTopoChange originalMeshToNewMesh(mesh);
diff --git a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
index 135a52520c..28aaf7ea06 100644
--- a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
+++ b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
@@ -25,6 +25,10 @@ Class
     Foam::polyMeshFilter
 
 Description
+    Filter the edges and faces of a polyMesh whilst satisfying the given mesh
+    quality criteria.
+
+    Works on a copy of the mesh.
 
 SourceFiles
     polyMeshFilter.C
@@ -102,10 +106,10 @@ class polyMeshFilter
         //  faces
         const scalar faceReductionFactor_;
 
-        //-
+        //- Maximum number of times a deleted point can be associated with the
+        //  creation of a bad face it is forced to be kept.
         const label maxPointErrorCount_;
 
-
         //- The minimum edge length for each edge
         scalarField minEdgeLen_;
 
@@ -195,6 +199,10 @@ class polyMeshFilter
 
 public:
 
+    //- Runtime type information
+    ClassName("polyMeshFilter");
+
+
     // Constructors
 
         //- Construct from fvMesh
@@ -225,6 +233,7 @@ public:
             //- Filter edges only.
             label filterEdges(const label nOriginalBadFaces);
 
+            //- Filter all faces that are in the face zone indirectPatchFaces
             label filterIndirectPatchFaces();
 };
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
index af911c7ad7..f6c35df41e 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
@@ -35,6 +35,12 @@ License
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
+namespace Foam
+{
+defineTypeNameAndDebug(edgeCollapser, 0);
+}
+
+
 Foam::label Foam::edgeCollapser::longestEdge
 (
     const face& f,
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
index 02359f39ac..dc6b1578c1 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
@@ -61,7 +61,7 @@ class face;
 class edge;
 
 /*---------------------------------------------------------------------------*\
-                           Class edgeCollapser Declaration
+                        Class edgeCollapser Declaration
 \*---------------------------------------------------------------------------*/
 
 class edgeCollapser
@@ -84,12 +84,18 @@ private:
         //- Reference to mesh
         const polyMesh& mesh_;
 
+        //- Controls collapse of a face to an edge
         const scalar guardFraction_;
 
+        //- Only collapse face to a point if high aspect ratio
         const scalar maxCollapseFaceToPointSideLengthCoeff_;
 
+        //- Allow a face to be collapsed to a point early, before the test
+        //  to collapse to an edge
         const Switch allowEarlyCollapseToPoint_;
 
+        //- Fraction of maxCollapseFaceToPointSideLengthCoeff_ to use when
+        //  allowEarlyCollapseToPoint_ is on
         const scalar allowEarlyCollapseCoeff_;
 
 
@@ -266,8 +272,8 @@ public:
                 const dictionary& meshQualityDict
             );
 
-            // Check mesh and mark points on faces in error
-            // Returns boolList with points in error set
+            //- Check mesh and mark points on faces in error
+            //  Returns boolList with points in error set
             static label checkMeshQuality
             (
                 const polyMesh& mesh,
@@ -300,7 +306,7 @@ public:
                 polyTopoChange& meshMod
             ) const;
 
-            // Mark (in collapseEdge) any edges to collapse
+            //- Mark (in collapseEdge) any edges to collapse
             label markSmallEdges
             (
                 const scalarField& minEdgeLen,
@@ -309,7 +315,7 @@ public:
                 Map<point>& collapsePointToLocation
             ) const;
 
-            // Mark (in collapseEdge) any edges to merge
+            //- Mark (in collapseEdge) any edges to merge
             label markMergeEdges
             (
                 const scalar maxCos,
@@ -332,6 +338,7 @@ public:
                 Map<point>& collapsePointToLocation
             ) const;
 
+            //- Marks edges in the faceZone indirectPatchFaces for collapse
             void markIndirectPatchFaces
             (
                 PackedBoolList& collapseEdge,
diff --git a/src/fvOptions/Make/files b/src/fvOptions/Make/files
index d6c2e970a5..23fdf0282f 100644
--- a/src/fvOptions/Make/files
+++ b/src/fvOptions/Make/files
@@ -28,11 +28,12 @@ $(derivedSources)/rotorDiskSource/trimModel/trimModel/trimModelNew.C
 $(derivedSources)/rotorDiskSource/trimModel/fixed/fixedTrim.C
 $(derivedSources)/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
 
-interRegion = $(derivedSources)/interRegionHeatTransferModel
-$(interRegion)/constantHeatTransfer/constantHeatTransfer.C
-$(interRegion)/interRegionHeatTransferModel/interRegionHeatTransferModel.C
-$(interRegion)/tabulatedHeatTransfer/tabulatedHeatTransfer.C
-$(interRegion)/variableHeatTransfer/variableHeatTransfer.C
+interRegion = sources/interRegion
+$(interRegion)/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
+$(interRegion)/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
+$(interRegion)/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
+$(interRegion)/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
+$(interRegion)/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
 
 
 /* constraints */
diff --git a/src/fvOptions/fvOptions/fvIOoptionList.C b/src/fvOptions/fvOptions/fvIOoptionList.C
index 1ea4be17d2..287e4171be 100644
--- a/src/fvOptions/fvOptions/fvIOoptionList.C
+++ b/src/fvOptions/fvOptions/fvIOoptionList.C
@@ -45,14 +45,15 @@ Foam::IOobject Foam::fv::IOoptionList::createIOobject
 
     if (io.headerOk())
     {
-        Info<< "Creating field source list from " << io.name() << nl << endl;
+        Info<< "Creating fintite volume options from " << io.name() << nl
+            << endl;
 
         io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
         return io;
     }
     else
     {
-        Info<< "No field sources present" << nl << endl;
+        Info<< "No finite volume options present" << nl << endl;
 
         io.readOpt() = IOobject::NO_READ;
         return io;
diff --git a/src/fvOptions/fvOptions/fvOption.C b/src/fvOptions/fvOptions/fvOption.C
index 4a66ded16e..8cb72b8e99 100644
--- a/src/fvOptions/fvOptions/fvOption.C
+++ b/src/fvOptions/fvOptions/fvOption.C
@@ -27,6 +27,7 @@ License
 #include "fvMesh.H"
 #include "fvMatrices.H"
 #include "volFields.H"
+#include "ListOps.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -86,9 +87,7 @@ void Foam::fv::option::setSelection(const dictionary& dict)
         }
         case smMapRegion:
         {
-            dict.lookup("nbrModelName") >> nbrModelName_;
             dict.lookup("nbrRegionName") >> nbrRegionName_;
-            master_ = readBool(dict.lookup("master"));
             break;
         }
         case smAll:
@@ -191,7 +190,10 @@ void Foam::fv::option::setCellSet()
                         (
                             mesh_,
                             nbrMesh,
-                            readBool(dict_.lookup("consistentMeshes"))
+                            meshToMeshNew::interpolationMethodNames_.read
+                            (
+                                dict_.lookup("interpolationMethod")
+                            )
                         )
                     );
                 }
@@ -247,7 +249,8 @@ Foam::fv::option::option
     const word& name,
     const word& modelType,
     const dictionary& dict,
-    const fvMesh& mesh
+    const fvMesh& mesh,
+    const bool master
 )
 :
     name_(name),
@@ -261,9 +264,8 @@ Foam::fv::option::option
     cellSetName_("none"),
     V_(0.0),
     meshInterpPtr_(),
-    nbrModelName_("none"),
     nbrRegionName_("none"),
-    master_(false),
+    master_(master),
     fieldNames_(),
     applied_()
 {
@@ -350,15 +352,7 @@ Foam::label Foam::fv::option::applyToField(const word& fieldName) const
         return 0;
     }
 
-    forAll(fieldNames_, i)
-    {
-        if (fieldNames_[i] == fieldName)
-        {
-            return i;
-        }
-    }
-
-    return -1;
+    return findIndex(fieldNames_, fieldName);
 }
 
 
diff --git a/src/fvOptions/fvOptions/fvOption.H b/src/fvOptions/fvOptions/fvOption.H
index 3501fdfa45..28509ee4c0 100644
--- a/src/fvOptions/fvOptions/fvOption.H
+++ b/src/fvOptions/fvOptions/fvOption.H
@@ -134,9 +134,6 @@ protected:
             //- Mesh to mesh interpolation object
             autoPtr<meshToMeshNew> meshInterpPtr_;
 
-            //- Name of the model in the neighbour mesh
-            word nbrModelName_;
-
             //- Name of the neighbour region to map
             word nbrRegionName_;
 
@@ -194,7 +191,8 @@ public:
             const word& name,
             const word& modelType,
             const dictionary& dict,
-            const fvMesh& mesh
+            const fvMesh& mesh,
+            const bool master = false
         );
 
         //- Return clone
@@ -288,9 +286,6 @@ public:
             //- Return const access to the total cell volume
             inline scalar V() const;
 
-            //- Return const access to the neighbour model name
-            inline const word& nbrModelName() const;
-
             //- Return const access to the neighbour region name
             inline const word& nbrRegionName() const;
 
diff --git a/src/fvOptions/fvOptions/fvOptionI.H b/src/fvOptions/fvOptions/fvOptionI.H
index aed4c14e55..796b2fdbc9 100644
--- a/src/fvOptions/fvOptions/fvOptionI.H
+++ b/src/fvOptions/fvOptions/fvOptionI.H
@@ -126,12 +126,6 @@ inline Foam::scalar& Foam::fv::option::duration()
 }
 
 
-inline const Foam::word& Foam::fv::option::nbrModelName() const
-{
-    return nbrModelName_;
-}
-
-
 inline const Foam::word& Foam::fv::option::nbrRegionName() const
 {
     return nbrRegionName_;
diff --git a/src/fvOptions/fvOptions/fvOptionIO.C b/src/fvOptions/fvOptions/fvOptionIO.C
index c89c97b864..c51b6b9b48 100644
--- a/src/fvOptions/fvOptions/fvOptionIO.C
+++ b/src/fvOptions/fvOptions/fvOptionIO.C
@@ -90,8 +90,11 @@ void Foam::fv::option::writeData(Ostream& os) const
 bool Foam::fv::option::read(const dictionary& dict)
 {
     active_ = readBool(dict.lookup("active"));
-    timeStart_ = readScalar(dict.lookup("timeStart"));
-    duration_  = readScalar(dict.lookup("duration"));
+
+    if (dict.readIfPresent("timeStart", timeStart_))
+    {
+        dict.lookup("duration") >> duration_;
+    }
 
     coeffs_ = dict.subDict(type() + "Coeffs");
 
diff --git a/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C b/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
new file mode 100644
index 0000000000..0c0cf7f284
--- /dev/null
+++ b/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C
@@ -0,0 +1,230 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "interRegionExplicitPorositySource.H"
+#include "fvMesh.H"
+#include "fvMatrices.H"
+#include "porosityModel.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(interRegionExplicitPorositySource, 0);
+    addToRunTimeSelectionTable
+    (
+        option,
+        interRegionExplicitPorositySource,
+        dictionary
+    );
+}
+}
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void Foam::fv::interRegionExplicitPorositySource::initialise()
+{
+    if (!firstIter_)
+    {
+        return;
+    }
+
+    const word zoneName(name_ + ".porous");
+
+    const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
+    const cellZoneMesh& cellZones = nbrMesh.cellZones();
+    label zoneID = cellZones.findZoneID(zoneName);
+
+    if (zoneID == -1)
+    {
+        cellZoneMesh& cz = const_cast<cellZoneMesh&>(cellZones);
+
+        zoneID = cz.size();
+
+        cz.setSize(zoneID + 1);
+
+        cz.set
+        (
+            zoneID,
+            new cellZone
+            (
+                zoneName,
+                nbrMesh.faceNeighbour(), // neighbour internal cells
+                zoneID,
+                cellZones
+            )
+        );
+
+        cz.clearAddressing();
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "void Foam::fv::interRegionExplicitPorositySource::initialise()"
+        )
+            << "Unable to create porous cellZone " << zoneName
+            << ": zone already exists"
+            << abort(FatalError);
+    }
+
+    porosityPtr_.reset
+    (
+        porosityModel::New
+        (
+            name_,
+            nbrMesh,
+            coeffs_,
+            zoneName
+        ).ptr()
+    ),
+
+    firstIter_ = false;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::interRegionExplicitPorositySource::interRegionExplicitPorositySource
+(
+    const word& name,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    option(name, modelType, dict, mesh, true),
+    porosityPtr_(NULL),
+    firstIter_(-1),
+    UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
+    rhoName_(coeffs_.lookupOrDefault<word>("rhoName", "rho")),
+    muName_(coeffs_.lookupOrDefault<word>("muName", "thermo:mu"))
+{
+    if (active_)
+    {
+        fieldNames_.setSize(1, UName_);
+        applied_.setSize(1, false);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fv::interRegionExplicitPorositySource::addSup
+(
+    fvMatrix<vector>& eqn,
+    const label fieldI
+)
+{
+    initialise();
+
+    const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
+
+    const volVectorField& U = eqn.psi();
+
+    volVectorField UNbr
+    (
+        IOobject
+        (
+            name_ + ".UNbr",
+            nbrMesh.time().timeName(),
+            nbrMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        nbrMesh,
+        dimensionedVector("zero", U.dimensions(), vector::zero)
+    );
+
+    // map local velocity onto neighbour region
+    meshInterp().mapSrcToTgt
+    (
+        U.internalField(),
+        plusEqOp<vector>(),
+        UNbr.internalField()
+    );
+
+    fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());
+
+    if (eqn.dimensions() == dimForce)
+    {
+        const volScalarField& rhoNbr =
+            nbrMesh.lookupObject<volScalarField>(rhoName_);
+        const volScalarField& muNbr =
+            nbrMesh.lookupObject<volScalarField>(muName_);
+
+        porosityPtr_->addResistance(nbrEqn, rhoNbr, muNbr);
+    }
+    else
+    {
+        porosityPtr_->addResistance(nbrEqn);
+    }
+
+    // convert source from neighbour to local region
+    fvMatrix<vector> porosityEqn(U, eqn.dimensions());
+    scalarField& Udiag = porosityEqn.diag();
+    vectorField& Usource = porosityEqn.source();
+
+    Udiag.setSize(eqn.diag().size(), 0.0);
+    Usource.setSize(eqn.source().size(), vector::zero);
+
+    meshInterp().mapTgtToSrc(nbrEqn.diag(), plusEqOp<scalar>(), Udiag);
+    meshInterp().mapTgtToSrc(nbrEqn.source(), plusEqOp<vector>(), Usource);
+
+    eqn -= porosityEqn;
+}
+
+
+void Foam::fv::interRegionExplicitPorositySource::writeData(Ostream& os) const
+{
+    os  << indent << name_ << endl;
+    dict_.write(os);
+}
+
+
+bool Foam::fv::interRegionExplicitPorositySource::read(const dictionary& dict)
+{
+    if (option::read(dict))
+    {
+        coeffs_.readIfPresent("UName", UName_);
+        coeffs_.readIfPresent("rhoName", rhoName_);
+        coeffs_.readIfPresent("muName", muName_);
+
+        // reset the porosity model?
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H b/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H
new file mode 100644
index 0000000000..da93314a0a
--- /dev/null
+++ b/src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H
@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::fv::interRegionExplicitPorositySource
+
+Description
+    Inter-region explicit porosity source
+
+    Sources described by, for example using the DarcyForchheimer model:
+
+        interRegionExplicitPorositySourceCoeffs
+        {
+            type            DarcyForchheimer;
+            DarcyForchheimerCoeffs
+            {
+                d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
+                f   f [0 -1 0 0 0 0 0] (0 0 0);
+
+                coordinateSystem
+                {
+                    e1  (0.70710678 0.70710678 0);
+                    e2  (0 0 1);
+                }
+            }
+        }
+
+Note:
+    The porous region must be selected as a cellZone.
+
+SourceFiles
+    interRegionExplicitPorositySource.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interRegionExplicitPorositySource_H
+#define interRegionExplicitPorositySource_H
+
+#include "fvOption.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class porosityModel;
+
+namespace fv
+{
+
+
+/*---------------------------------------------------------------------------*\
+                   Class interRegionExplicitPorositySource Declaration
+\*---------------------------------------------------------------------------*/
+
+class interRegionExplicitPorositySource
+:
+    public option
+{
+
+protected:
+
+    // Protected data
+
+        //- Run-time selectable porosity model
+        autoPtr<porosityModel> porosityPtr_;
+
+        //- First iteration
+        bool firstIter_;
+
+        //- Velocity field name, default = U
+        word UName_;
+
+        //- Density field name (compressible case only), default = rho
+        word rhoName_;
+
+        //- Dynamic viscosity field name (compressible case only)
+        //  default = thermo:mu
+        word muName_;
+
+
+    // Protected Member Functions
+
+        //- Initialise
+        void initialise();
+
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        interRegionExplicitPorositySource
+        (
+            const interRegionExplicitPorositySource&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const interRegionExplicitPorositySource&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("interRegionExplicitPorositySource");
+
+
+    // Constructors
+
+        //- Construct from components
+        interRegionExplicitPorositySource
+        (
+            const word& name,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~interRegionExplicitPorositySource()
+    {}
+
+
+    // Member Functions
+
+        // Add explicit and implicit contributions
+
+            //- Vector
+            virtual void addSup
+            (
+                fvMatrix<vector>& eqn,
+                const label fieldI
+            );
+
+
+        // I-O
+
+            //- Write data
+            virtual void writeData(Ostream&) const;
+
+            //- Read dictionary
+            virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
similarity index 87%
rename from src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
rename to src/fvOptions/sources/interRegion/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
index cc8879877c..1cafec08c4 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.C
@@ -91,14 +91,7 @@ Foam::fv::constantHeatTransfer::constantHeatTransfer
             )
         );
 
-        const DimensionedField<scalar, volMesh>& htcConsti =
-            htcConst_().dimensionedInternalField();
-        const DimensionedField<scalar, volMesh>& AoVi =
-            AoV_().dimensionedInternalField();
-        dimensionedScalar interVol("V", dimVolume, meshInterp().V());
-
-        htc_.dimensionedInternalField() = htcConsti*AoVi*interVol/mesh.V();
-        htc_.correctBoundaryConditions();
+        htc_ = htcConst_()*AoV_();
     }
 }
 
@@ -111,10 +104,9 @@ Foam::fv::constantHeatTransfer::~constantHeatTransfer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::tmp<Foam::volScalarField>
-Foam::fv::constantHeatTransfer::calculateHtc()
+void Foam::fv::constantHeatTransfer::calculateHtc()
 {
-    return htc_;
+    // do nothing
 }
 
 
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H
similarity index 98%
rename from src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H
rename to src/fvOptions/sources/interRegion/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H
index 52facea986..ceb4aaf15a 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/constantHeatTransfer/constantHeatTransfer.H
@@ -87,7 +87,7 @@ public:
     // Public Functions
 
         //- Calculate the heat transfer coefficient
-        virtual const tmp<volScalarField> calculateHtc();
+        virtual void calculateHtc();
 
 
         // I-O
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
similarity index 90%
rename from src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
rename to src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
index 7cfde633e1..710a43cf3e 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C
@@ -41,7 +41,7 @@ namespace fv
 }
 
 
-// * * * * * * * * * * * *  Private member functions * * * * * * * * * * * //
+// * * * * * * * * * * * *  Protected member functions * * * * * * * * * * * //
 
 void Foam::fv::interRegionHeatTransferModel::setNbrModel()
 {
@@ -71,7 +71,7 @@ void Foam::fv::interRegionHeatTransferModel::setNbrModel()
 
     if (!nbrModelFound)
     {
-        FatalErrorIn("interRegionHeatTransferModel::check()")
+        FatalErrorIn("interRegionHeatTransferModel::setNbrModel()")
             << "Neighbour model not found" << nbrModelName_
             << " in region " << nbrMesh.name() << nl
             << exit(FatalError);
@@ -81,6 +81,24 @@ void Foam::fv::interRegionHeatTransferModel::setNbrModel()
 }
 
 
+void Foam::fv::interRegionHeatTransferModel::correct()
+{
+    if (master_)
+    {
+        if (mesh_.time().timeIndex() != timeIndex_)
+        {
+            calculateHtc();
+            timeIndex_ = mesh_.time().timeIndex();
+        }
+    }
+    else
+    {
+        nbrModel().correct();
+        interpolate(nbrModel().htc(), htc_);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fv::interRegionHeatTransferModel::interRegionHeatTransferModel
@@ -91,9 +109,11 @@ Foam::fv::interRegionHeatTransferModel::interRegionHeatTransferModel
     const fvMesh& mesh
 )
 :
-    option(name, modelType, dict, mesh),
+    option(name, modelType, dict, mesh, readBool(dict.lookup("master"))),
     nbrModel_(NULL),
+    nbrModelName_(word::null),
     firstIter_(true),
+    timeIndex_(-1),
     htc_
     (
         IOobject
@@ -119,6 +139,8 @@ Foam::fv::interRegionHeatTransferModel::interRegionHeatTransferModel
 {
     if (active())
     {
+        coeffs_.lookup("nbrModelName") >> nbrModelName_;
+
         coeffs_.lookup("fieldNames") >> fieldNames_;
         applied_.setSize(fieldNames_.size(), false);
 
@@ -143,6 +165,8 @@ void Foam::fv::interRegionHeatTransferModel::addSup
 {
     setNbrModel();
 
+    correct();
+
     const volScalarField& h = eqn.psi();
 
     const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);
@@ -172,11 +196,6 @@ void Foam::fv::interRegionHeatTransferModel::addSup
 
     interpolate(Tnbr, Tmapped.internalField());
 
-    if (!master_)
-    {
-        interpolate(nbrModel().calculateHtc()(), htc_);
-    }
-
     if (debug)
     {
         Info<< "Volumetric integral of htc: "
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H
similarity index 93%
rename from src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H
rename to src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H
index dd7104de13..55aabd45e1 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.H
@@ -65,22 +65,21 @@ private:
         //- Pointer to neighbour interRegionHeatTransferModel
         interRegionHeatTransferModel* nbrModel_;
 
+        //- Name of the model in the neighbour mesh
+        word nbrModelName_;
+
         //- First iteration
         bool firstIter_;
 
-
-    // Private members
-
-        //- Set the neighbour interRegionHeatTransferModel
-        void setNbrModel();
+        //- Time index - used for updating htc
+        label timeIndex_;
 
 
 protected:
 
     // Protected data
 
-        //- Heat transfer coefficient [W/m2/k] by area/volume [1/m]
-        // registered on the master mesh
+        //- Heat transfer coefficient [W/m2/k] times area/volume [1/m]
         volScalarField htc_;
 
         //- Flag to activate semi-implicit coupling
@@ -95,6 +94,12 @@ protected:
 
     // Protected member functions
 
+        //- Set the neighbour interRegionHeatTransferModel
+        void setNbrModel();
+
+        //- Correct to calculate the inter-region heat transfer coefficient
+        void correct();
+
         //- Interpolate field with nbrModel specified
         template<class Type>
         tmp<Field<Type> > interpolate
@@ -168,7 +173,7 @@ public:
             virtual void addSup(fvMatrix<scalar>& eqn, const label fieldI);
 
             //- Calculate heat transfer coefficient
-            virtual const tmp<volScalarField> calculateHtc() = 0;
+            virtual void calculateHtc() = 0;
 
 
         // I-O
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelI.H b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelI.H
similarity index 100%
rename from src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelI.H
rename to src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelI.H
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelTemplates.C b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelTemplates.C
similarity index 100%
rename from src/fvOptions/sources/derived/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelTemplates.C
rename to src/fvOptions/sources/interRegion/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModelTemplates.C
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
similarity index 96%
rename from src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
rename to src/fvOptions/sources/interRegion/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
index cf8ec5a67e..fd49c31963 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.C
@@ -109,8 +109,7 @@ Foam::fv::tabulatedHeatTransfer::~tabulatedHeatTransfer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::tmp<Foam::volScalarField>
-Foam::fv::tabulatedHeatTransfer::calculateHtc()
+void Foam::fv::tabulatedHeatTransfer::calculateHtc()
 {
     const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName());
 
@@ -130,9 +129,7 @@ Foam::fv::tabulatedHeatTransfer::calculateHtc()
         htcc[i] = hTable()(mag(U[i]), UMagNbrMapped[i]);
     }
 
-    htcc = htcc*AoV()*meshInterp().V()/mesh_.V();
-
-    return htc_;
+    htcc = htcc*AoV();
 }
 
 
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H
similarity index 98%
rename from src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H
rename to src/fvOptions/sources/interRegion/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H
index 4656560eac..8feb9d68df 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/tabulatedHeatTransfer/tabulatedHeatTransfer.H
@@ -104,7 +104,7 @@ public:
     // Public Functions
 
         //- Calculate the heat transfer coefficient
-        virtual const tmp<volScalarField> calculateHtc();
+        virtual void calculateHtc();
 
 
         // I-O
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
similarity index 96%
rename from src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
rename to src/fvOptions/sources/interRegion/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
index 4445306bdd..cf96fa0433 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.C
@@ -97,8 +97,7 @@ Foam::fv::variableHeatTransfer::~variableHeatTransfer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::tmp<Foam::volScalarField>
-Foam::fv::variableHeatTransfer::calculateHtc()
+void Foam::fv::variableHeatTransfer::calculateHtc()
 {
     const fvMesh& nbrMesh =
         mesh_.time().lookupObject<fvMesh>(nbrRegionName());
@@ -120,9 +119,7 @@ Foam::fv::variableHeatTransfer::calculateHtc()
 
     const scalarField htcNbrMapped(interpolate(htcNbr));
 
-    htc_.internalField() = htcNbrMapped*AoV_*meshInterp().V()/mesh_.V();
-
-    return htc_;
+    htc_.internalField() = htcNbrMapped*AoV_;
 }
 
 
diff --git a/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H
similarity index 98%
rename from src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H
rename to src/fvOptions/sources/interRegion/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H
index 8180f63a34..57e42be811 100644
--- a/src/fvOptions/sources/derived/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransferModel/variableHeatTransfer/variableHeatTransfer.H
@@ -108,7 +108,7 @@ public:
     // Public Functions
 
         //- Calculate the heat transfer coefficient
-        virtual const tmp<volScalarField> calculateHtc();
+        virtual void calculateHtc();
 
 
         // I-O
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
index 6d717bb174..6b59696660 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
@@ -383,15 +383,25 @@ template<class CloudType>
 inline Foam::scalar Foam::ThermoCloud<CloudType>::Tmax() const
 {
     scalar T = -GREAT;
+    scalar n = 0;
     forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
         T = max(T, p.T());
+        n++;
     }
 
     reduce(T, maxOp<scalar>());
+    reduce(n, sumOp<label>());
 
-    return max(0.0, T);
+    if (n > 0)
+    {
+        return T;
+    }
+    else
+    {
+        return 0.0;
+    }
 }
 
 
@@ -399,15 +409,25 @@ template<class CloudType>
 inline Foam::scalar Foam::ThermoCloud<CloudType>::Tmin() const
 {
     scalar T = GREAT;
+    scalar n = 0;
     forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
         T = min(T, p.T());
+        n++;
     }
 
     reduce(T, minOp<scalar>());
+    reduce(n, sumOp<label>());
 
-    return max(0.0, T);
+    if (n > 0)
+    {
+        return T;
+    }
+    else
+    {
+        return 0.0;
+    }
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index e9bc815304..ed60ded959 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -510,6 +510,9 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     const CompositionModel<reactingCloudType>& composition =
         td.cloud().composition();
 
+    const scalar TMax = td.cloud().phaseChange().TMax(pc_, this->Tc_);
+    const scalar Tdash = min(T, TMax);
+    const scalar Tsdash = min(Ts, TMax);
 
     // Calculate mass transfer due to phase change
     td.cloud().phaseChange().calculate
@@ -520,8 +523,8 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
         Pr,
         d,
         nus,
-        T,
-        Ts,
+        Tdash,
+        Tsdash,
         pc_,
         this->Tc_,
         YComponents,
@@ -541,7 +544,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
         const label idc = composition.localToGlobalCarrierId(idPhase, i);
         const label idl = composition.globalIds(idPhase)[i];
 
-        const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, T);
+        const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, Tdash);
         Sh -= dMassPC[i]*dh/dt;
     }
 
@@ -558,12 +561,12 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
             const label idc = composition.localToGlobalCarrierId(idPhase, i);
             const label idl = composition.globalIds(idPhase)[i];
 
-            const scalar Cp = composition.carrier().Cp(idc, pc_, Ts);
+            const scalar Cp = composition.carrier().Cp(idc, pc_, Tsdash);
             const scalar W = composition.carrier().W(idc);
             const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
 
             const scalar Dab =
-                composition.liquids().properties()[idl].D(pc_, Ts, Wc);
+                composition.liquids().properties()[idl].D(pc_, Tsdash, Wc);
 
             // Molar flux of species coming from the particle (kmol/m^2/s)
             N += Ni;
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
index f1dc7c4ba0..10bb48d887 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -404,20 +404,23 @@ void Foam::ParticleCollector<CloudType>::write()
     }
 
 
-    Field<scalar> faceMassTotal(mass_.size());
-    Field<scalar> faceMassFlowRate(massFlowRate_.size());
+    Field<scalar> faceMassTotal(mass_.size(), 0.0);
+    this->getModelProperty("massTotal", faceMassTotal);
+
+    Field<scalar> faceMassFlowRate(massFlowRate_.size(), 0.0);
+    this->getModelProperty("massFlowRate", faceMassFlowRate);
 
     forAll(faces_, faceI)
     {
         scalarList allProcMass(Pstream::nProcs());
         allProcMass[procI] = massTotal_[faceI];
         Pstream::gatherList(allProcMass);
-        faceMassTotal[faceI] = sum(allProcMass);
+        faceMassTotal[faceI] += sum(allProcMass);
 
         scalarList allProcMassFlowRate(Pstream::nProcs());
         allProcMassFlowRate[procI] = massFlowRate_[faceI];
         Pstream::gatherList(allProcMassFlowRate);
-        faceMassFlowRate[faceI] = sum(allProcMassFlowRate);
+        faceMassFlowRate[faceI] += sum(allProcMassFlowRate);
 
         Info<< "    face " << faceI
             << ": total mass = " << faceMassTotal[faceI]
@@ -470,20 +473,25 @@ void Foam::ParticleCollector<CloudType>::write()
 
     if (resetOnWrite_)
     {
-        forAll(faces_, faceI)
-        {
-            massFlowRate_[faceI] = 0.0;
-        }
+        Field<scalar> dummy(faceMassTotal.size(), 0.0);
+        this->setModelProperty("massTotal", dummy);
+        this->setModelProperty("massFlowRate", dummy);
+
         timeOld_ = timeNew;
         totalTime_ = 0.0;
     }
+    else
+    {
+        this->setModelProperty("massTotal", faceMassTotal);
+        this->setModelProperty("massFlowRate", faceMassFlowRate);
+    }
 
     forAll(faces_, faceI)
     {
         mass_[faceI] = 0.0;
+        massTotal_[faceI] = 0.0;
+        massFlowRate_[faceI] = 0.0;
     }
-
-    // writeProperties();
 }
 
 
@@ -552,8 +560,8 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
         (
             "Foam::ParticleCollector<CloudType>::ParticleCollector"
             "("
-                "const dictionary& dict,"
-                "CloudType& owner"
+                "const dictionary&,"
+                "CloudType&"
             ")"
         )
             << "Unknown mode " << mode << ".  Available options are "
@@ -565,8 +573,6 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
     massFlowRate_.setSize(faces_.size(), 0.0);
 
     makeLogFile(faces_, points_, area_);
-
-    // readProperties(); AND initialise mass... fields
 }
 
 
@@ -579,6 +585,7 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
     CloudFunctionObject<CloudType>(pc),
     mode_(pc.mode_),
     parcelType_(pc.parcelType_),
+    removeCollected_(pc.removeCollected_),
     points_(pc.points_),
     faces_(pc.faces_),
     faceTris_(pc.faceTris_),
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
index 97b84a6049..b46ea2b5bf 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -200,8 +200,8 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::dh
 (
     const label idc,
     const label idl,
-    const label p,
-    const label T
+    const scalar p,
+    const scalar T
 ) const
 {
     scalar dh = 0;
@@ -230,8 +230,8 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::dh
                 "("
                     "const label, "
                     "const label, "
-                    "const label, "
-                    "const label"
+                    "const scalar, "
+                    "const scalar"
                 ") const"
             )   << "Unknown enthalpyTransfer type" << abort(FatalError);
         }
@@ -241,4 +241,21 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::dh
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::LiquidEvaporation<CloudType>::TMax
+(
+    const scalar pIn,
+    const scalar TIn
+) const
+{
+    scalar T = -GREAT;
+    forAll(liquids_, i)
+    {
+        T = max(T, liquids_.properties()[i].pv(pIn, TIn));
+    }
+
+    return T;
+}
+
+
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
index 60db6c63a5..5a97dc7394 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,9 +127,12 @@ public:
         (
             const label idc,
             const label idl,
-            const label p,
-            const label T
+            const scalar p,
+            const scalar T
         ) const;
+
+        //- Return maximum/limiting temperature
+        virtual scalar TMax(const scalar pIn, const scalar TIn) const;
 };
 
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C
index 803b849fa4..8a88620992 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -293,8 +293,8 @@ Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::dh
 (
     const label idc,
     const label idl,
-    const label p,
-    const label T
+    const scalar p,
+    const scalar T
 ) const
 {
     scalar dh = 0;
@@ -329,8 +329,8 @@ Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::dh
                 "("
                     "const label, "
                     "const label, "
-                    "const label, "
-                    "const label"
+                    "const scalar, "
+                    "const scalar"
                 ") const"
             )   << "Unknown enthalpyTransfer type" << abort(FatalError);
         }
@@ -340,4 +340,21 @@ Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::dh
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::TMax
+(
+    const scalar pIn,
+    const scalar TIn
+) const
+{
+    scalar T = -GREAT;
+    forAll(liquids_, i)
+    {
+        T = max(T, liquids_.properties()[i].pv(pIn, TIn));
+    }
+
+    return T;
+}
+
+
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H
index 48d92292cd..bf9ca867ce 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -137,9 +137,12 @@ public:
         (
             const label idc,
             const label idl,
-            const label p,
-            const label T
+            const scalar p,
+            const scalar T
         ) const;
+
+        //- Return maximum/limiting temperature
+        virtual scalar TMax(const scalar pIn, const scalar TIn) const;
 };
 
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C
index c84c7b83f8..7141ff9814 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -170,14 +170,25 @@ Foam::scalar Foam::PhaseChangeModel<CloudType>::dh
 (
     const label idc,
     const label idl,
-    const label p,
-    const label T
+    const scalar p,
+    const scalar T
 ) const
 {
     return 0.0;
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::PhaseChangeModel<CloudType>::TMax
+(
+    const scalar,
+    const scalar
+) const
+{
+    return GREAT;
+}
+
+
 template<class CloudType>
 void Foam::PhaseChangeModel<CloudType>::addToPhaseChangeMass(const scalar dMass)
 {
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
index 3ee8d4dbd3..ac6c5c06b3 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -180,10 +180,12 @@ public:
         (
             const label idc,
             const label idl,
-            const label p,
-            const label T
+            const scalar p,
+            const scalar T
         ) const;
 
+        //- Return maximum/limiting temperature
+        virtual scalar TMax(const scalar pIn, const scalar TIn) const;
 
         //- Add to phase change mass
         void addToPhaseChangeMass(const scalar dMass);
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
index 3e1a29e988..2d47253de9 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -572,7 +572,7 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
     regionModels::surfaceFilmModels::surfaceFilmModel& filmModel =
         const_cast<regionModels::surfaceFilmModels::surfaceFilmModel&>
         (
-            this->owner().db().objectRegistry::template
+            this->owner().db().time().objectRegistry::template
                 lookupObject<regionModels::surfaceFilmModels::surfaceFilmModel>
                 (
                     "surfaceFilmProperties"
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.C
index 75ea3867d1..e3d36bfcb3 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,9 +43,11 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const pointData& wDist)
     }
 }
 
+
 Foam::Istream& Foam::operator>>(Istream& is, pointData& wDist)
 {
     return is >> static_cast<pointEdgePoint&>(wDist) >> wDist.s_ >> wDist.v_;
 }
 
+
 // ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H
index b215732e3e..c2a83752a2 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,6 +148,13 @@ public:
                 TrackingData& td
             );
 
+    // Member Operators
+
+        // Needed for List IO
+        inline bool operator==(const pointData&) const;
+        inline bool operator!=(const pointData&) const;
+
+
     // IOstream Operators
 
         friend Ostream& operator<<(Ostream&, const pointData&);
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H
index bab3c51e3c..b00ca55097 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -215,4 +215,23 @@ inline bool Foam::pointData::updateEdge
 }
 
 
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::pointData::operator==(const Foam::pointData& rhs)
+const
+{
+    return
+        pointEdgePoint::operator==(rhs)
+     && (s() == rhs.s())
+     && (v() == rhs.v());
+}
+
+
+inline bool Foam::pointData::operator!=(const Foam::pointData& rhs)
+const
+{
+    return !(*this == rhs);
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 69cdf234cd..04c5786ea6 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -258,6 +258,106 @@ void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 }
 
 
+void Foam::meshRefinement::testSyncPointList
+(
+    const string& msg,
+    const polyMesh& mesh,
+    const List<scalar>& fld
+)
+{
+    if (fld.size() != mesh.nPoints())
+    {
+        FatalErrorIn
+        (
+            "meshRefinement::testSyncPointList(const polyMesh&"
+            ", const List<scalar>&)"
+        )   << msg << endl
+            << "fld size:" << fld.size() << " mesh points:" << mesh.nPoints()
+            << abort(FatalError);
+    }
+
+    Pout<< "Checking field " << msg << endl;
+    scalarField minFld(fld);
+    syncTools::syncPointList
+    (
+        mesh,
+        minFld,
+        minEqOp<scalar>(),
+        GREAT
+    );
+    scalarField maxFld(fld);
+    syncTools::syncPointList
+    (
+        mesh,
+        maxFld,
+        maxEqOp<scalar>(),
+        -GREAT
+    );
+    forAll(minFld, pointI)
+    {
+        const scalar& minVal = minFld[pointI];
+        const scalar& maxVal = maxFld[pointI];
+        if (mag(minVal-maxVal) > SMALL)
+        {
+            Pout<< msg << " at:" << mesh.points()[pointI] << nl
+                << "    minFld:" << minVal << nl
+                << "    maxFld:" << maxVal << nl
+                << endl;
+        }
+    }
+}
+
+
+void Foam::meshRefinement::testSyncPointList
+(
+    const string& msg,
+    const polyMesh& mesh,
+    const List<point>& fld
+)
+{
+    if (fld.size() != mesh.nPoints())
+    {
+        FatalErrorIn
+        (
+            "meshRefinement::testSyncPointList(const polyMesh&"
+            ", const List<point>&)"
+        )   << msg << endl
+            << "fld size:" << fld.size() << " mesh points:" << mesh.nPoints()
+            << abort(FatalError);
+    }
+
+    Pout<< "Checking field " << msg << endl;
+    pointField minFld(fld);
+    syncTools::syncPointList
+    (
+        mesh,
+        minFld,
+        minMagSqrEqOp<point>(),
+        point(GREAT, GREAT, GREAT)
+    );
+    pointField maxFld(fld);
+    syncTools::syncPointList
+    (
+        mesh,
+        maxFld,
+        maxMagSqrEqOp<point>(),
+        vector::zero
+    );
+    forAll(minFld, pointI)
+    {
+        const point& minVal = minFld[pointI];
+        const point& maxVal = maxFld[pointI];
+        if (mag(minVal-maxVal) > SMALL)
+        {
+            Pout<< msg << " at:" << mesh.points()[pointI] << nl
+                << "    minFld:" << minVal << nl
+                << "    maxFld:" << maxVal << nl
+                << endl;
+        }
+    }
+}
+
+
 void Foam::meshRefinement::checkData()
 {
     Pout<< "meshRefinement::checkData() : Checking refinement structure."
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index c265f05257..899311d909 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -879,6 +879,20 @@ public:
             //- Debugging: check that all faces still obey start()>end()
             void checkData();
 
+            static void testSyncPointList
+            (
+                const string& msg,
+                const polyMesh& mesh,
+                const List<scalar>& fld
+            );
+
+            static void testSyncPointList
+            (
+                const string& msg,
+                const polyMesh& mesh,
+                const List<point>& fld
+            );
+
             //- Compare two lists over all boundary faces
             template<class T>
             void testSyncBoundaryFaceList
diff --git a/src/meshTools/algorithms/PointEdgeWave/pointEdgePoint.H b/src/meshTools/algorithms/PointEdgeWave/pointEdgePoint.H
index d4f4aa062f..5c323b0a3b 100644
--- a/src/meshTools/algorithms/PointEdgeWave/pointEdgePoint.H
+++ b/src/meshTools/algorithms/PointEdgeWave/pointEdgePoint.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,8 +25,8 @@ Class
     Foam::pointEdgePoint
 
 Description
-    Holds information regarding nearest wall point. Used in pointEdgeWave.
-    (so not standard meshWave)
+    Holds information regarding nearest wall point. Used in PointEdgeWave.
+    (so not standard FaceCellWave)
     To be used in wall distance calculation.
 
 SourceFiles
@@ -116,7 +116,7 @@ public:
             inline scalar distSqr() const;
 
 
-        // Needed by meshWave
+        // Needed by PointEdgeWave
 
             //- Check whether origin has been changed at all or
             //  still contains original (invalid) value.
diff --git a/src/meshTools/algorithms/PointEdgeWave/pointEdgePointI.H b/src/meshTools/algorithms/PointEdgeWave/pointEdgePointI.H
index f9fe6ca9e3..b466b26828 100644
--- a/src/meshTools/algorithms/PointEdgeWave/pointEdgePointI.H
+++ b/src/meshTools/algorithms/PointEdgeWave/pointEdgePointI.H
@@ -310,14 +310,14 @@ inline bool Foam::pointEdgePoint::equal
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::pointEdgePoint::operator==(const Foam::pointEdgePoint& rhs)
- const
+const
 {
-    return origin() == rhs.origin();
+    return (origin() == rhs.origin()) && (distSqr() == rhs.distSqr());
 }
 
 
 inline bool Foam::pointEdgePoint::operator!=(const Foam::pointEdgePoint& rhs)
- const
+const
 {
     return !(*this == rhs);
 }
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
index c0f192f7de..692096b8e9 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
@@ -54,9 +54,9 @@ Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol
     const tetPoints& tetB
 ) const
 {
-    tetPointRef::tetIntersectionList insideTets;
+    static tetPointRef::tetIntersectionList insideTets;
     label nInside = 0;
-    tetPointRef::tetIntersectionList cutInsideTets;
+    static tetPointRef::tetIntersectionList cutInsideTets;
     label nCutInside = 0;
 
     tetPointRef::storeOp inside(insideTets, nInside);
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
index a77459f82f..f8b689757d 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,6 +77,8 @@ Description
                 mean            on;
                 prime2Mean      on;
                 base            time;
+                window          10.0;
+                windowName      w1;
             }
             p
             {
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverageTemplates.C b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverageTemplates.C
index 0fdf344531..f4b03dd4b4 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverageTemplates.C
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverageTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,7 +42,15 @@ void Foam::fieldAverage::addMeanField
 
         const word& fieldName = faItems_[fieldI].fieldName();
 
-        const word meanFieldName = fieldName + EXT_MEAN;
+        word meanFieldName = fieldName + EXT_MEAN;
+        if
+        (
+            (faItems_[fieldI].window() > 0)
+         && (faItems_[fieldI].windowName() != "")
+        )
+        {
+            meanFieldName = meanFieldName + "_" + faItems_[fieldI].windowName();
+        }
 
         Info<< "Reading/calculating field " << meanFieldName << nl << endl;
 
@@ -100,7 +108,16 @@ void Foam::fieldAverage::addPrime2MeanField
 
         const word& fieldName = faItems_[fieldI].fieldName();
 
-        const word meanFieldName = fieldName + EXT_PRIME2MEAN;
+        word meanFieldName = fieldName + EXT_PRIME2MEAN;
+        if
+        (
+            (faItems_[fieldI].window() > 0)
+         && (faItems_[fieldI].windowName() != "")
+        )
+        {
+            meanFieldName = meanFieldName + "_" + faItems_[fieldI].windowName();
+        }
+
         Info<< "Reading/calculating field " << meanFieldName << nl << endl;
 
         if (obr_.foundObject<fieldType2>(meanFieldName))
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.C b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.C
index 39bbfb4fa9..3b01b21594 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.C
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -54,7 +54,8 @@ Foam::fieldAverageItem::fieldAverageItem()
     mean_(0),
     prime2Mean_(0),
     base_(ITER),
-    window_(-1.0)
+    window_(-1.0),
+    windowName_("")
 {}
 
 
@@ -64,7 +65,8 @@ Foam::fieldAverageItem::fieldAverageItem(const fieldAverageItem& faItem)
     mean_(faItem.mean_),
     prime2Mean_(faItem.prime2Mean_),
     base_(faItem.base_),
-    window_(faItem.window_)
+    window_(faItem.window_),
+    windowName_(faItem.windowName_)
 {}
 
 
@@ -94,6 +96,7 @@ void Foam::fieldAverageItem::operator=(const fieldAverageItem& rhs)
     prime2Mean_ = rhs.prime2Mean_;
     base_ = rhs.base_;
     window_ = rhs.window_;
+    windowName_ = rhs.windowName_;
 }
 
 
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H
index 992b5e393f..7fe21a9aef 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ Description
         prime2Mean      on;
         base            time; // iteration
         window          200;  // optional averaging window
+        windowName      w1;   // optional window name (default = "")
     }
     \endverbatim
 
@@ -107,6 +108,9 @@ private:
         //- Averaging window - defaults to -1 for 'all iters/time'
         scalar window_;
 
+        //- Averaging window name - defaults to 'window'
+        word windowName_;
+
 
 public:
 
@@ -171,6 +175,11 @@ public:
                 return window_;
             }
 
+            const word& windowName() const
+            {
+                return windowName_;
+            }
+
 
     // Member Operators
 
@@ -190,7 +199,8 @@ public:
              && a.mean_ == b.mean_
              && a.prime2Mean_ == b.prime2Mean_
              && a.base_ == b.base_
-             && a.window_ == b.window_;
+             && a.window_ == b.window_
+             && a.windowName_ == b.windowName_;
         }
 
         friend bool operator!=
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemIO.C b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemIO.C
index 4ebad485a7..af293ff989 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemIO.C
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemIO.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,6 +46,7 @@ Foam::fieldAverageItem::fieldAverageItem(Istream& is)
     entry.lookup("prime2Mean") >> prime2Mean_;
     base_ = baseTypeNames_[entry.lookup("base")];
     window_ = entry.lookupOrDefault<scalar>("window", -1.0);
+    windowName_ = entry.lookupOrDefault<word>("windowName", "");
 }
 
 
@@ -66,6 +67,7 @@ Foam::Istream& Foam::operator>>(Istream& is, fieldAverageItem& faItem)
     entry.lookup("prime2Mean") >> faItem.prime2Mean_;
     faItem.base_ = faItem.baseTypeNames_[entry.lookup("base")];
     faItem.window_ = entry.lookupOrDefault<scalar>("window", -1.0);
+    faItem.windowName_ = entry.lookupOrDefault<word>("windowName", "");
 
     return is;
 }
@@ -90,6 +92,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const fieldAverageItem& faItem)
     {
         os.writeKeyword("window") << faItem.window_
             << token::END_STATEMENT << nl;
+
+        if (faItem.windowName_ != "")
+        {
+            os.writeKeyword("windowName") << faItem.windowName_
+                << token::END_STATEMENT << nl;
+        }
     }
 
     os  << token::END_BLOCK << nl;
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
index 9344f75ba0..b99c0629cc 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
@@ -27,6 +27,9 @@ Class
 Description
     mesh to mesh interpolation class.
 
+Note
+    This class is due to be deprecated in favour of meshToMeshNew
+
 SourceFiles
     meshToMesh.C
     calculateMeshToMeshAddressing.C
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C
index b71e768f42..9ecda6cb4f 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.C
@@ -39,6 +39,20 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(meshToMeshNew, 0);
+
+    template<>
+    const char* Foam::NamedEnum
+    <
+        Foam::meshToMeshNew::interpolationMethod,
+        2
+    >::names[] =
+    {
+        "map",
+        "cellVolumeWeight"
+    };
+
+    const NamedEnum<meshToMeshNew::interpolationMethod, 2>
+        meshToMeshNew::interpolationMethodNames_;
 }
 
 Foam::scalar Foam::meshToMeshNew::tolerance_ = 1e-6;
@@ -707,22 +721,40 @@ void Foam::meshToMeshNew::calcAddressing
     }
 
 
-    if (directMapping_)
+    switch (method_)
     {
-        calcDirect(src, tgt, srcSeedI, tgtSeedI);
-    }
-    else
-    {
-        calcIndirect
-        (
-            src,
-            tgt,
-            srcSeedI,
-            tgtSeedI,
-            srcCellIDs,
-            mapFlag,
-            startSeedI
-        );
+        case imMap:
+        {
+            calcDirect(src, tgt, srcSeedI, tgtSeedI);
+            break;
+        }
+        case imCellVolumeWeight:
+        {
+            calcIndirect
+            (
+                src,
+                tgt,
+                srcSeedI,
+                tgtSeedI,
+                srcCellIDs,
+                mapFlag,
+                startSeedI
+            );
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "void Foam::meshToMeshNew::calcAddressing"
+                "("
+                    "const polyMesh&, "
+                    "const polyMesh&"
+                ")"
+            )
+                << "Unknown interpolation method"
+                << abort(FatalError);
+        }
     }
 
 
@@ -739,7 +771,7 @@ Foam::meshToMeshNew::meshToMeshNew
 (
     const polyMesh& src,
     const polyMesh& tgt,
-    const bool directMapping
+    const interpolationMethod& method
 )
 :
     srcRegionName_(src.name()),
@@ -748,11 +780,11 @@ Foam::meshToMeshNew::meshToMeshNew
     tgtToSrcCellAddr_(),
     srcToTgtCellWght_(),
     tgtToSrcCellWght_(),
+    method_(method),
     V_(0.0),
     singleMeshProc_(-1),
     srcMapPtr_(NULL),
-    tgtMapPtr_(NULL),
-    directMapping_(directMapping)
+    tgtMapPtr_(NULL)
 {
     Info<< "Creating mesh-to-mesh addressing for " << src.name()
         << " and " << tgt.name() << " regions" << endl;
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H
index 6b5211ef43..85f0c1332e 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNew.H
@@ -40,6 +40,7 @@ SourceFiles
 #include "boundBox.H"
 #include "mapDistribute.H"
 #include "volFieldsFwd.H"
+#include "NamedEnum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,6 +53,22 @@ namespace Foam
 
 class meshToMeshNew
 {
+public:
+
+    // Public data types
+
+        //- Enumeration specifying required accuracy
+        enum interpolationMethod
+        {
+            imMap,
+            imCellVolumeWeight
+        };
+
+        static const NamedEnum<interpolationMethod, 2>
+            interpolationMethodNames_;
+
+private:
+
     // Private data
 
         //- Name of source mesh region
@@ -72,6 +89,9 @@ class meshToMeshNew
         //- Target to source cell interpolation weights
         scalarListList tgtToSrcCellWght_;
 
+        //- Interpolation method
+        interpolationMethod method_;
+
         //- Cell total volume in overlap region [m3]
         scalar V_;
 
@@ -85,9 +105,6 @@ class meshToMeshNew
         //- Target map pointer - parallel running only
         autoPtr<mapDistribute> tgtMapPtr_;
 
-        //- Flag to indicate that direct (one-to-one) mapping should be applied
-        bool directMapping_;
-
         //- Tolerance used in volume overlap calculations
         static scalar tolerance_;
 
@@ -289,7 +306,7 @@ public:
     (
         const polyMesh& src,
         const polyMesh& tgt,
-        const bool directMapping
+        const interpolationMethod& method
     );
 
 
diff --git a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C
index 4587b8ea60..69bbc69667 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMeshNew/meshToMeshNewTemplates.C
@@ -395,28 +395,33 @@ void Foam::meshToMeshNew::interpolate
 
     if (interpPatches)
     {
-        if (directMapping_)
+        switch (method_)
         {
-            result.boundaryField() == field.boundaryField();
-        }
-        else
-        {
-            notImplemented
-            (
-                "void Foam::meshToMeshNew::interpolate"
-                "("
-                    "const GeometricField<Type, fvPatchField, volMesh>&, "
-                    "const CombineOp&, "
-                    "GeometricField<Type, fvPatchField, volMesh>&, "
-                    "const bool"
-                ") const - non-conformal patches"
-            )
+            case imMap:
+            {
+                result.boundaryField() == field.boundaryField();
+                break;
+            }
+            default:
+            {
+                notImplemented
+                (
+                    "void Foam::meshToMeshNew::interpolate"
+                    "("
+                        "const GeometricField<Type, fvPatchField, volMesh>&, "
+                        "const CombineOp&, "
+                        "GeometricField<Type, fvPatchField, volMesh>&, "
+                        "const bool"
+                    ") const - non-conformal patches"
+                )
 
-            // do something...
+                // do something...
+            }
         }
     }
 }
 
+
 template<class Type, class CombineOp>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
 Foam::meshToMeshNew::interpolate
diff --git a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
index 3edfd4ce38..08dd1e3764 100644
--- a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -312,7 +312,30 @@ Foam::radiation::greyMeanAbsorptionEmission::ECont(const label bandI) const
     {
         const volScalarField& dQ =
             mesh_.lookupObject<volScalarField>("dQ");
-        E().internalField() = EhrrCoeff_*dQ;
+
+        if (dQ.dimensions() == dimEnergy/dimTime)
+        {
+            E().internalField() = EhrrCoeff_*dQ/mesh_.V();
+        }
+        else if (dQ.dimensions() == dimEnergy/dimTime/dimVolume)
+        {
+            E().internalField() = EhrrCoeff_*dQ;
+        }
+        else
+        {
+            if (debug)
+            {
+                WarningIn
+                (
+                    "tmp<volScalarField>"
+                    "radiation::greyMeanAbsorptionEmission::ECont"
+                    "("
+                        "const label"
+                    ") const"
+                )
+                    << "Incompatible dimensions for dQ field" << endl;
+            }
+        }
     }
 
     return E;
diff --git a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.C b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.C
index ead9b263fc..99bb98acb4 100644
--- a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.C
+++ b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.C
@@ -137,8 +137,82 @@ gasSpecies() const
 template<class ReactionThermo>
 void Foam::solidReaction<ReactionThermo>::write(Ostream& os) const
 {
-    Reaction<ReactionThermo>::write(os);
+    OStringStream reaction;
+    os.writeKeyword("reaction") << solidReactionStr(reaction)
+        << token::END_STATEMENT << nl;
 }
 
 
+template<class ReactionThermo>
+Foam::string Foam::solidReaction<ReactionThermo>::solidReactionStr
+(
+    OStringStream& reaction
+) const
+{
+    this->reactionStrLeft(reaction);
+    reaction << " + ";
+    solidReactionStrLeft(reaction);
+    reaction << " = ";
+    this->reactionStrRight(reaction);
+    reaction << " + ";
+    solidReactionStrRight(reaction);
+    return reaction.str();
+
+}
+
+
+template<class ReactionThermo>
+void Foam::solidReaction<ReactionThermo>::solidReactionStrLeft
+(
+    OStringStream& reaction
+) const
+{
+    for (label i = 0; i < glhs().size(); ++i)
+    {
+        reaction << " + ";
+
+        if (i > 0)
+        {
+            reaction << " + ";
+        }
+        if (mag(glhs()[i].stoichCoeff - 1) > SMALL)
+        {
+            reaction << glhs()[i].stoichCoeff;
+        }
+        reaction << gasSpecies()[glhs()[i].index];
+        if (mag(glhs()[i].exponent - glhs()[i].stoichCoeff) > SMALL)
+        {
+            reaction << "^" << glhs()[i].exponent;
+        }
+    }
+}
+
+
+template<class ReactionThermo>
+void Foam::solidReaction<ReactionThermo>::solidReactionStrRight
+(
+    OStringStream& reaction
+) const
+{
+
+    for (label i = 0; i < grhs().size(); ++i)
+    {
+        reaction << " + ";
+
+        if (i > 0)
+        {
+            reaction << " + ";
+        }
+        if (mag(grhs()[i].stoichCoeff - 1) > SMALL)
+        {
+            reaction << grhs()[i].stoichCoeff;
+        }
+        reaction << gasSpecies()[grhs()[i].index];
+        if (mag(grhs()[i].exponent - grhs()[i].stoichCoeff) > SMALL)
+        {
+            reaction << "^" << grhs()[i].exponent;
+        }
+    }
+}
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H
index 1bc4a82243..cc9eb6a7b7 100644
--- a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H
+++ b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H
@@ -81,6 +81,15 @@ private:
     // Private Member Functions
 
 
+        //- Return string representation of reaction
+        string solidReactionStr(OStringStream&) const;
+
+        //- Return string representation of the left of the reaction
+        void solidReactionStrLeft(OStringStream&) const;
+
+        //- Return string representation of the right of the reaction
+        void solidReactionStrRight(OStringStream&) const;
+
         //- Disallow default bitwise assignment
         void operator=(const solidReaction&);
 
diff --git a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReactionI.H b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReactionI.H
index 534bb66fa7..9bad3ba778 100644
--- a/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReactionI.H
+++ b/src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReactionI.H
@@ -39,7 +39,8 @@ inline Ostream& operator<<
     const solidReaction<ReactionThermo>& r
 )
 {
-    r.write(os);
+    OStringStream reaction;
+    os << r.solidReactionStr(reaction)<< token::END_STATEMENT <<nl;
     return os;
 }
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
index 8f2c39c37e..2a8a5ca1e1 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
@@ -31,21 +31,14 @@ License
 template<class ReactionThermo>
 Foam::label Foam::Reaction<ReactionThermo>::nUnNamedReactions = 0;
 
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
 template<class ReactionThermo>
-Foam::label Foam::Reaction<ReactionThermo>::getNewReactionID()
+void Foam::Reaction<ReactionThermo>::reactionStrLeft
+(
+    OStringStream& reaction
+) const
 {
-    return nUnNamedReactions++;
-}
-
-
-template<class ReactionThermo>
-Foam::string Foam::Reaction<ReactionThermo>::reactionStr() const
-{
-    OStringStream reaction;
-
     for (label i = 0; i < lhs_.size(); ++i)
     {
         if (i > 0)
@@ -62,28 +55,15 @@ Foam::string Foam::Reaction<ReactionThermo>::reactionStr() const
             reaction << "^" << lhs_[i].exponent;
         }
     }
+}
 
-    for (label i = 0; i < glhs().size(); ++i)
-    {
-        reaction << " + ";
-
-        if (i > 0)
-        {
-            reaction << " + ";
-        }
-        if (mag(glhs()[i].stoichCoeff - 1) > SMALL)
-        {
-            reaction << glhs()[i].stoichCoeff;
-        }
-        reaction << gasSpecies()[glhs()[i].index];
-        if (mag(glhs()[i].exponent - glhs()[i].stoichCoeff) > SMALL)
-        {
-            reaction << "^" << glhs()[i].exponent;
-        }
-    }
-
-    reaction << " = ";
 
+template<class ReactionThermo>
+void Foam::Reaction<ReactionThermo>::reactionStrRight
+(
+    OStringStream& reaction
+) const
+{
     for (label i = 0; i < rhs_.size(); ++i)
     {
         if (i > 0)
@@ -100,26 +80,27 @@ Foam::string Foam::Reaction<ReactionThermo>::reactionStr() const
             reaction << "^" << rhs_[i].exponent;
         }
     }
+}
 
-    for (label i = 0; i < grhs().size(); ++i)
-    {
-        reaction << " + ";
 
-        if (i > 0)
-        {
-            reaction << " + ";
-        }
-        if (mag(grhs()[i].stoichCoeff - 1) > SMALL)
-        {
-            reaction << grhs()[i].stoichCoeff;
-        }
-        reaction << gasSpecies()[grhs()[i].index];
-        if (mag(grhs()[i].exponent - grhs()[i].stoichCoeff) > SMALL)
-        {
-            reaction << "^" << grhs()[i].exponent;
-        }
-    }
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<class ReactionThermo>
+Foam::label Foam::Reaction<ReactionThermo>::getNewReactionID()
+{
+    return nUnNamedReactions++;
+}
+
+
+template<class ReactionThermo>
+Foam::string Foam::Reaction<ReactionThermo>::reactionStr
+(
+    OStringStream& reaction
+) const
+{
+    reactionStrLeft(reaction);
+    reaction << " = ";
+    reactionStrRight(reaction);
     return reaction.str();
 }
 
@@ -464,7 +445,9 @@ Foam::Reaction<ReactionThermo>::New
 template<class ReactionThermo>
 void Foam::Reaction<ReactionThermo>::write(Ostream& os) const
 {
-    os.writeKeyword("reaction") << reactionStr() << token::END_STATEMENT << nl;
+    OStringStream reaction;
+    os.writeKeyword("reaction") << reactionStr(reaction)
+        << token::END_STATEMENT << nl;
 }
 
 
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
index b0c66fb86f..2679fc1bbb 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
@@ -66,6 +66,16 @@ class Reaction
 :
     public ReactionThermo
 {
+protected:
+
+    // Protected member functions
+
+        //- Return string representation of the left of the reaction
+        void reactionStrLeft(OStringStream& reaction) const;
+
+        //- Return string representation of the right of the reaction
+        void reactionStrRight(OStringStream& reaction) const;
+
 
 public:
 
@@ -134,7 +144,7 @@ private:
     // Private Member Functions
 
         //- Return string representation of reaction
-        string reactionStr() const;
+        string reactionStr(OStringStream& reaction) const;
 
         //- Construct reaction thermo
         void setThermo(const HashPtrTable<ReactionThermo>& thermoDatabase);
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H
index 7e380b4c80..6dece80ba7 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H
@@ -67,8 +67,8 @@ Reaction<ReactionThermo>::rhs() const
 template<class ReactionThermo>
 inline Ostream& operator<<(Ostream& os, const Reaction<ReactionThermo>& r)
 {
-   os << r.reactionStr()<< token::END_STATEMENT <<nl;
-
+    OStringStream reaction;
+    os << r.reactionStr(reaction)<< token::END_STATEMENT <<nl;
    return os;
 }
 
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/U b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/U
index 7d6875ae04..23454ad6de 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/U
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/0.org/porous/U
@@ -17,14 +17,14 @@ FoamFile
 
 dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform (0.1 0 0);
+internalField   uniform (0.01 0 0);
 
 boundaryField
 {
     inlet
     {
         type            fixedValue;
-        value           uniform (0.1 0 0);
+        value           uniform (0.01 0 0);
     }
     outlet
     {
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun-serial b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun
similarity index 100%
rename from tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun-serial
rename to tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/Allrun
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/polyMesh/blockMeshDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/polyMesh/blockMeshDict
index 49c36ebacf..d0439108c7 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/polyMesh/blockMeshDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/air/polyMesh/blockMeshDict
@@ -52,15 +52,15 @@ vertices
 
 blocks
 (
-    hex (0 1 5 4 12 13 17 16) (20 20 20) simpleGrading (1 1 1)
-    hex (1 2 6 5 13 14 18 17) (20 20 20) simpleGrading (1 1 1)
-    hex (2 3 7 6 14 15 19 18) (20 20 20) simpleGrading (1 1 1)
-    hex (3 0 4 7 15 12 16 19) (20 20 20) simpleGrading (1 1 1)
-    hex (4 5 9 8 16 17 21 20) cylinder (20 20 20) simpleGrading (1 1 1)
-    hex (5 6 10 9 17 18 22 21) cylinder (20 20 20) simpleGrading (1 1 1)
-    hex (6 7 11 10 18 19 23 22) cylinder (20 20 20) simpleGrading (1 1 1)
-    hex (7 4 8 11 19 16 20 23) cylinder (20 20 20) simpleGrading (1 1 1)
-    hex (8 9 10 11 20 21 22 23) innerCylinder (20 20 20) simpleGrading (1 1 1)
+    hex (0 1  5  4 12 13 17 16) (30 30 30) simpleGrading (1 1 1)
+    hex (1 2  6  5 13 14 18 17) (30 30 30) simpleGrading (1 1 1)
+    hex (2 3  7  6 14 15 19 18) (30 30 30) simpleGrading (1 1 1)
+    hex (3 0  4  7 15 12 16 19) (30 30 30) simpleGrading (1 1 1)
+    hex (4 5  9  8 16 17 21 20) cylinder (30 30 30) simpleGrading (1 1 1)
+    hex (5 6 10  9 17 18 22 21) cylinder (30 30 30) simpleGrading (1 1 1)
+    hex (6 7 11 10 18 19 23 22) cylinder (30 30 30) simpleGrading (1 1 1)
+    hex (7 4  8 11 19 16 20 23) cylinder (30 30 30) simpleGrading (1 1 1)
+    hex (8 9 10 11 20 21 22 23) innerCylinder (30 30 30) simpleGrading (1 1 1)
 );
 
 edges
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/polyMesh/blockMeshDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/polyMesh/blockMeshDict
index 86e616bd27..f32b1c7621 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/polyMesh/blockMeshDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/constant/porous/polyMesh/blockMeshDict
@@ -18,20 +18,20 @@ convertToMeters 0.05;
 
 vertices
 (
-    (-2 4 4)     // 0
-    (12 4 4)     // 1
-    (12 6 4)     // 2
-    (-2 6 4)     // 3
+    (-2 2 6)     // 0
+    (12 2 6)     // 1
+    (12 8 6)     // 2
+    (-2 8 6)     // 3
 
-    (-2 4 6)     // 4
-    (12 4 6)     // 5
-    (12 6 6)     // 6
-    (-2 6 6)     // 7
+    (-2 2 7)     // 4
+    (12 2 7)     // 5
+    (12 8 7)     // 6
+    (-2 8 7)     // 7
 );
 
 blocks
 (
-    hex (0 1 2 3 4 5 6 7) (20 4 4) simpleGrading (1 1 1)
+    hex (0 1 2 3 4 5 6 7) (40 30 30) simpleGrading (1 1 1)
 );
 
 edges
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvOptions b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvOptions
index 8f81bdae42..00978429d9 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvOptions
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvOptions
@@ -20,18 +20,44 @@ airToporous
     type            constantHeatTransfer;
     active          on;
     selectionMode   mapRegion;
-    consistentMeshes no;
-    nbrModelName    porousToair;
+    interpolationMethod cellVolumeWeight;
     nbrRegionName   porous;
     master          false;
 
     constantHeatTransferCoeffs
     {
+        nbrModelName    porousToair;
         fieldNames      (h);
         semiImplicit    no;
     }
 }
 
+porosityBlockage
+{
+    type            interRegionExplicitPorositySource;
+    active          on;
+    selectionMode   mapRegion;
+    interpolationMethod cellVolumeWeight;
+    nbrRegionName   porous;
+
+    interRegionExplicitPorositySourceCoeffs
+    {
+        type            DarcyForchheimer;
+
+        DarcyForchheimerCoeffs
+        {
+            d   d [0 -2 0 0 0] (10 -1000 -1000);
+            f   f [0 -1 0 0 0] (0 0 0);
+
+            coordinateSystem
+            {
+                e1  (0 1 0);
+                e2  (0 0 1);
+            }
+        }
+    }
+}
+
 MRF1
 {
     type            MRFSource;
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvSolution b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvSolution
index 4f83fc7823..36d404271f 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/air/fvSolution
@@ -55,14 +55,14 @@ relaxationFactors
     fields
     {
         rho             1;
-        p_rgh           0.5;
+        p_rgh           0.7;
     }
     equations
     {
-        U               0.2;
-        "(h|e)"         0.2;
-        k               0.2;
-        epsilon         0.2;
+        U               0.3;
+        "(h|e)"         0.3;
+        k               0.3;
+        epsilon         0.3;
     }
 }
 
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvOptions b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvOptions
index 3c5bd05176..f2f64e4882 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvOptions
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvOptions
@@ -20,13 +20,13 @@ porousToair
     type            constantHeatTransfer;
     active          on;
     selectionMode   mapRegion;
-    consistentMeshes no;
-    nbrModelName    airToporous;
+    interpolationMethod cellVolumeWeight;
     nbrRegionName   air;
     master          true;
 
     constantHeatTransferCoeffs
     {
+        nbrModelName    airToporous;
         fieldNames      (h);
         semiImplicit    no;
     }
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvSolution b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvSolution
index 1b3b041a2a..fd6da4016c 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/heatExchanger/system/porous/fvSolution
@@ -53,14 +53,14 @@ relaxationFactors
     fields
     {
         rho             1;
-        p_rgh           0.5;
+        p_rgh           0.7;
     }
     equations
     {
-        U               0.2;
-        "(h|e)"         0.2;
-        k               0.2;
-        epsilon         0.2;
+        U               0.3;
+        "(h|e)"         0.3;
+        k               0.3;
+        epsilon         0.3;
     }
 }
 
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
index 42ddef8235..bcb3e3e80d 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
@@ -133,7 +133,24 @@ subModels
 
 
 cloudFunctions
-{}
+{
+    particleCollector
+    {
+        mode            concentricCircle;
+
+        origin          (0.05 0.025 0.005);
+        radius          (0.01 0.025 0.05);
+        nSector         10;
+        refDir          (1 0 0);
+        normal          (0 0 1);
+
+        negateParcelsOppositeNormal no;
+        removeCollected no;
+        surfaceFormat   vtk;
+        resetOnWrite    no;
+        log             yes;
+    }
+}
 
 
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/H2O b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/H2O
new file mode 100644
index 0000000000..baa8dd2152
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/H2O
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      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/simpleReactingParcelFoam/verticalChannel/0.org/T b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/T
new file mode 100644
index 0000000000..df744edb03
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/T
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      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/simpleReactingParcelFoam/verticalChannel/0.org/U b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/U
new file mode 100644
index 0000000000..401e7ced02
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/U
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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;
+        rhoInlet        1.2;
+        massFlowRate    constant 0.00379;
+        value           uniform (0 14.68 0);
+    }
+    inletSides
+    {
+        type            flowRateInletVelocity;
+        rhoInlet        1.2;
+        massFlowRate    constant 0.00832;
+        value           uniform (0 17.79 0);
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+    }
+    walls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/air b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/air
new file mode 100644
index 0000000000..81541afef9
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/air
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      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/simpleReactingParcelFoam/verticalChannel/0.org/alphat b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/alphat
new file mode 100644
index 0000000000..33d260bf7c
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/alphat
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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/simpleReactingParcelFoam/verticalChannel/0.org/k b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/k
new file mode 100644
index 0000000000..2a76936893
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/k
@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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-9;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.15;
+        value           uniform 3.75e-9;
+    }
+    inletSides
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.16;
+        value           uniform 3.75e-9;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 3.75e-9;
+    }
+    walls
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/mut b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/mut
new file mode 100644
index 0000000000..7cfeaae133
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/mut
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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/simpleReactingParcelFoam/verticalChannel/0.org/omega b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/omega
new file mode 100644
index 0000000000..efd2924775
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/omega
@@ -0,0 +1,62 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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 4.5e-3;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            compressible::turbulentMixingLengthFrequencyInlet;
+        mixingLength    0.007;
+        k               k;
+        value           uniform 4.5e-3;
+    }
+    inletSides
+    {
+        type            compressible::turbulentMixingLengthFrequencyInlet;
+        mixingLength    0.007;
+        k               k;
+        value           uniform 4.5e-3;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 4.5e-3;
+    }
+    walls
+    {
+        type            compressible::omegaWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/p b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/p
new file mode 100644
index 0000000000..921f06d7a1
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0.org/p
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      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/simpleReactingParcelFoam/verticalChannel/0/H2O b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/H2O
new file mode 100644
index 0000000000..baa8dd2152
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/H2O
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      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/simpleReactingParcelFoam/verticalChannel/0/T b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/T
new file mode 100644
index 0000000000..df744edb03
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/T
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      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/simpleReactingParcelFoam/verticalChannel/0/U b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/U
new file mode 100644
index 0000000000..401e7ced02
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/U
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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;
+        rhoInlet        1.2;
+        massFlowRate    constant 0.00379;
+        value           uniform (0 14.68 0);
+    }
+    inletSides
+    {
+        type            flowRateInletVelocity;
+        rhoInlet        1.2;
+        massFlowRate    constant 0.00832;
+        value           uniform (0 17.79 0);
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+    }
+    walls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/air b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/air
new file mode 100644
index 0000000000..81541afef9
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/air
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      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/simpleReactingParcelFoam/verticalChannel/0/alphat b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/alphat
new file mode 100644
index 0000000000..33d260bf7c
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/alphat
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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/simpleReactingParcelFoam/verticalChannel/0/k b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/k
new file mode 100644
index 0000000000..2a76936893
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/k
@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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-9;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.15;
+        value           uniform 3.75e-9;
+    }
+    inletSides
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.16;
+        value           uniform 3.75e-9;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 3.75e-9;
+    }
+    walls
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/mut b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/mut
new file mode 100644
index 0000000000..7cfeaae133
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/mut
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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/simpleReactingParcelFoam/verticalChannel/0/omega b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/omega
new file mode 100644
index 0000000000..efd2924775
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/omega
@@ -0,0 +1,62 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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 4.5e-3;
+
+boundaryField
+{
+    back
+    {
+        type            symmetryPlane;
+    }
+    front
+    {
+        type            symmetryPlane;
+    }
+    inletCentral
+    {
+        type            compressible::turbulentMixingLengthFrequencyInlet;
+        mixingLength    0.007;
+        k               k;
+        value           uniform 4.5e-3;
+    }
+    inletSides
+    {
+        type            compressible::turbulentMixingLengthFrequencyInlet;
+        mixingLength    0.007;
+        k               k;
+        value           uniform 4.5e-3;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 4.5e-3;
+    }
+    walls
+    {
+        type            compressible::omegaWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/p b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/p
new file mode 100644
index 0000000000..921f06d7a1
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/0/p
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      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/simpleReactingParcelFoam/verticalChannel/Allclean b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allclean
new file mode 100755
index 0000000000..d0839b002e
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allclean
@@ -0,0 +1,15 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+# remove old time and post-processing folders
+rm -rf 0 *[1-9]* processor* postProcessing
+
+# copy 0.org to 0
+cp -r 0.org 0
+
+cleanCase
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allrun b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allrun
new file mode 100755
index 0000000000..36e8545b84
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/Allrun
@@ -0,0 +1,19 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+rm -rf 0
+cp -rf 0.org 0
+
+runApplication blockMesh
+
+runApplication potentialFoam
+
+# remove incompatible (volumetric) flux field
+rm -f 0/phi
+
+runApplication $(getApplication)
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/RASProperties b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/RASProperties
new file mode 100644
index 0000000000..568821b00d
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/RASProperties
@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      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/simpleReactingParcelFoam/verticalChannel/constant/chemistryProperties b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/chemistryProperties
new file mode 100644
index 0000000000..0fbd9de5dc
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/chemistryProperties
@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      chemistryProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+chemistryType
+{
+    chemistrySolver   noChemistrySolver;
+    chemistryThermo   rho;
+}
+
+chemistry       off;
+
+initialChemicalTimeStep 1e-07;  // NOT USED
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/combustionProperties b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/combustionProperties
new file mode 100644
index 0000000000..5ede7572c3
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/combustionProperties
@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      combustionProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+combustionModel  PaSR<rhoChemistryCombustion>;
+
+active  false;
+
+PaSRCoeffs
+{
+    Cmix                0.1;
+    turbulentReaction   off;
+    useReactionRate     true;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/g b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/g
new file mode 100644
index 0000000000..e0ac2653b5
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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 -9.81 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/particleTrackDict b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/particleTrackDict
new file mode 100644
index 0000000000..bf4c329856
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/particleTrackDict
@@ -0,0 +1,28 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      particleTrackDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+cloudName       reactingCloud1Tracks;
+
+fields
+(
+    d
+    U
+    T
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/polyMesh/blockMeshDict b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000..47e0eb9f8c
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/polyMesh/blockMeshDict
@@ -0,0 +1,223 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      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;
+}
+
+boundary
+(
+    back
+    {
+        type symmetryPlane;
+        faces
+        (
+            (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)
+        );
+    }
+
+    front
+    {
+        type symmetryPlane;
+        faces
+        (
+            (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)
+        );
+    }
+
+    inletCentral
+    {
+        type patch;
+        faces
+        (
+            (13 14 32 31)
+        );
+    }
+
+    inletSides
+    {
+        type patch;
+        faces
+        (
+            (17 0 18 35)
+            (9 10 28 27)
+        );
+    }
+
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (3 4 22 21)
+            (4 5 23 22)
+            (5 6 24 23)
+        );
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/polyMesh/boundary b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/polyMesh/boundary
new file mode 100644
index 0000000000..1319b62326
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/polyMesh/boundary
@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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;
+        inGroups        1(symmetryPlane);
+        nFaces          9340;
+        startFace       265900;
+    }
+    front
+    {
+        type            symmetryPlane;
+        inGroups        1(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/simpleReactingParcelFoam/verticalChannel/constant/radiationProperties b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/radiationProperties
new file mode 100644
index 0000000000..0ebf90015a
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/radiationProperties
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation       off;
+
+radiationModel  none;
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
new file mode 100644
index 0000000000..aa8b02241a
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
@@ -0,0 +1,206 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      reactingCloud1Properties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solution
+{
+    active          yes;
+
+    transient       no; // yes;
+    calcFrequency   10;
+    maxTrackTime    5.0;
+    maxCo           0.3;
+
+    coupled         true;
+    cellValueSourceCorrection off;
+
+    sourceTerms
+    {
+        resetOnStartup  false;
+        schemes
+        {
+            rho             semiImplicit 1;
+            U               semiImplicit 1;
+            Yi              semiImplicit 1;
+            h               semiImplicit 1;
+            radiation       semiImplicit 1;
+        }
+    }
+
+    interpolationSchemes
+    {
+        rho             cell;
+        U               cellPoint;
+        thermo:mu       cell;
+        T               cell;
+        Cp              cell;
+        p               cell;
+    }
+
+    integrationSchemes
+    {
+        U               Euler;
+        T               analytical;
+    }
+}
+
+
+constantProperties
+{
+    rho0            1000;
+    T0              350;
+    Cp0             4100;
+
+    youngsModulus   1e9;
+    poissonsRatio   0.35;
+
+    epsilon0        1;
+    f0              0.5;
+
+    Tvap            273;
+    Tbp             373;
+    Pr              0.7;
+    LDevol          0;
+    hRetentionCoeff 1;
+
+    constantVolume  false;
+}
+
+
+subModels
+{
+    particleForces
+    {
+        sphereDrag;
+        gravity;
+    }
+
+    injectionModels
+    {
+        model1
+        {
+            type            patchInjection;
+            massFlowRate    0.8e-03;
+            parcelBasisType mass;
+            patchName       inletCentral;
+            parcelsPerSecond 100;
+            duration        1; // NOTE: set to 1 for steady state
+            U0              (0 40 0);
+            flowRateProfile constant 1;
+            sizeDistribution
+            {
+                type        general;
+                generalDistribution
+                {
+                    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)
+                    );
+                }
+            }
+        }
+    }
+
+    dispersionModel stochasticDispersionRAS;
+
+    patchInteractionModel standardWallInteraction;
+
+    heatTransferModel RanzMarshall;
+
+    compositionModel singleMixtureFraction;
+
+    phaseChangeModel liquidEvaporation;
+
+    devolatilisationModel none;
+
+    surfaceReactionModel none;
+
+    surfaceFilmModel none;
+
+    radiation       off;
+
+    standardWallInteractionCoeffs
+    {
+        type            rebound;
+    }
+
+    RanzMarshallCoeffs
+    {
+        BirdCorrection  off;
+    }
+
+    singleMixtureFractionCoeffs
+    {
+        phases
+        (
+            gas
+            {
+            }
+            liquid
+            {
+                H2O 1;
+            }
+            solid
+            {
+            }
+        );
+        YGasTot0        0;
+        YLiquidTot0     1;
+        YSolidTot0      0;
+    }
+
+    liquidEvaporationCoeffs
+    {
+        enthalpyTransfer enthalpyDifference;
+        activeLiquids   ( H2O );
+    }
+}
+
+
+cloudFunctions
+{
+    patchPostProcessing
+    {
+        maxStoredParcels 100;
+        patches         ( outlet );
+    }
+
+    particleTracks
+    {
+        trackInterval   5;
+        maxSamples      1000000;
+        resetOnWrite    yes;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/reactions b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/reactions
new file mode 100644
index 0000000000..228f5f836b
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/reactions
@@ -0,0 +1,8 @@
+species
+(
+    air
+    H2O
+);
+
+reactions
+{}
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/thermo.incompressiblePoly b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/thermo.incompressiblePoly
new file mode 100644
index 0000000000..3b8bf27dfe
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/thermo.incompressiblePoly
@@ -0,0 +1,91 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermo.incompressiblePoly;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+N2
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       28.0134;
+    }
+    equationOfState
+    {
+        rhoCoeffs<8>    ( 3.8936 -0.016463 3.2101e-05 -2.9174e-08 9.9889e-12 0 0 0 );
+    }
+    thermodynamics
+    {
+        Hf              0;
+        Sf              0;
+        CpCoeffs<8>     ( 979.08 0.41787 -0.0011761 1.6742e-06 -7.2559e-10 0 0 0 );
+    }
+    transport
+    {
+        muCoeffs<8>     ( 1.5068e-06 6.1598e-08 -1.8188e-11 0 0 0 0 0 );
+        kappaCoeffs<8>  ( 0.0031494 8.4997e-05 -1.2621e-08 0 0 0 0 0 );
+    }
+}
+
+H2O
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       18.0153;
+    }
+    equationOfState
+    {
+        rhoCoeffs<8>    ( 2.5039 -0.010587 2.0643e-05 -1.8761e-08 6.4237e-12 0 0 0 );
+    }
+    thermodynamics
+    {
+        Hf              -13423000;
+        Sf              10482;
+        CpCoeffs<8>     ( 1563.1 1.604 -0.0029334 3.2168e-06 -1.1571e-09 0 0 0 );
+    }
+    transport
+    {
+        muCoeffs<8>     ( 1.5068e-06 6.1598e-08 -1.8188e-11 0 0 0 0 0 );
+        kappaCoeffs<8>  ( 0.0037972 0.00015336 -1.1859e-08 0 0 0 0 0 );
+    }
+}
+
+air
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       28.85;
+    }
+    equationOfState
+    {
+        rhoCoeffs<8>    ( 4.0097 -0.016954 3.3057e-05 -3.0042e-08 1.0286e-11 0 0 0 );
+    }
+    thermodynamics
+    {
+        Hf              0;
+        Sf              0;
+        CpCoeffs<8>     ( 948.76 0.39171 -0.00095999 1.393e-06 -6.2029e-10 0 0 0 );
+    }
+    transport
+    {
+        muCoeffs<8>     ( 1.5061e-06 6.16e-08 -1.819e-11 0 0 0 0 0 );
+        kappaCoeffs<8>  ( 0.0025219 8.506e-05 -1.312e-08 0 0 0 0 0 );
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/thermophysicalProperties b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/thermophysicalProperties
new file mode 100644
index 0000000000..9b63bc18b3
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/thermophysicalProperties
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         reactingMixture;
+    transport       polynomial;
+    thermo          hPolynomial;
+    energy          sensibleEnthalpy;
+    equationOfState icoPolynomial;
+    specie          specie;
+}
+
+chemistryReader foamChemistryReader;
+
+foamChemistryFile "$FOAM_CASE/constant/reactions";
+
+foamChemistryThermoFile "$FOAM_CASE/constant/thermo.incompressiblePoly";
+
+inertSpecie     air;
+
+liquids
+{
+    H2O
+    {
+        defaultCoeffs   yes;
+    }
+}
+
+solids
+{
+    // none
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/turbulenceProperties b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/turbulenceProperties
new file mode 100644
index 0000000000..aaccd5feb0
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/constant/turbulenceProperties
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RASModel;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/controlDict b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/controlDict
new file mode 100644
index 0000000000..273cb3b23d
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/controlDict
@@ -0,0 +1,72 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     simpleReactingParcelFoam;
+
+startFoam       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         500;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   20;
+
+purgeWrite      10;
+
+writeFormat     ascii;
+
+writePrecision  10;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+
+functions
+{
+    faceSource1
+    {
+        type            faceSource;
+        functionObjectLibs ("libfieldFunctionObjects.so");
+        enabled         yes;
+        outputControl   outputTime;
+        log             yes;
+        valueOutput     no;
+        source          patch;
+        sourceName      outlet;
+        operation       weightedAverage;
+        weightField     phi;
+        fields
+        (
+            H2O
+            T
+        );
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/fvOptions b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/fvOptions
new file mode 100644
index 0000000000..82e5c6b4b0
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/fvOptions
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvOptions;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// none
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/fvSchemes b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/fvSchemes
new file mode 100644
index 0000000000..45390e50dc
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/fvSchemes
@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      bounded Gauss upwind;
+    div(phid,p)     bounded Gauss upwind;
+    div(phi,K)      bounded Gauss linear;
+    div(phi,h)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,omega) bounded Gauss upwind;
+    div(phi,Yi_h)   Gauss upwind;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/fvSolution b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/fvSolution
new file mode 100644
index 0000000000..a1c90c36c6
--- /dev/null
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/fvSolution
@@ -0,0 +1,79 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    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.05;
+
+        smoother         DICGaussSeidel;
+        nPreSweeps       0;
+        nPostSweeps      2;
+
+        cacheAgglomeration true;
+
+        nCellsInCoarsestLevel 10;
+        agglomerator     faceAreaPair;
+        mergeLevels      1;
+
+        maxIter          50;
+    };
+
+    "(U|Yi|h|k|omega)"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       0;
+        relTol          0.1;
+        maxIter         20;
+    }
+}
+
+potentialFlow
+{
+    nNonOrthogonalCorrectors 5;
+}
+
+SIMPLE
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    rhoMin          rhoMin [1 -3 0 0 0] 0.1;
+    rhoMax          rhoMax [1 -3 0 0 0] 1.5;
+}
+
+relaxationFactors
+{
+    fields
+    {
+        p               0.3;
+        rho             1;
+    }
+    equations
+    {
+        U               0.7;
+        h               0.7;
+        ".*"            0.7;
+    }
+}
+
+
+// ************************************************************************* //