diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options b/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options
index 11053f31a9..502938c53c 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options
@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I../rhoPorousMRFPimpleFoam \
     -I.. \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
index 65cfcb4516..4067c125eb 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
@@ -104,13 +104,9 @@ int main(int argc, char *argv[])
             }
 
             rho = thermo.rho();
+        }
 
-            runTime.write();
-        }
-        else
-        {
-            runTime.write();
-        }
+        runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
diff --git a/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C b/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C
index 44a02e234a..2a958b3a3c 100644
--- a/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C
+++ b/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C
@@ -269,9 +269,10 @@ void selectCurvatureCells
         querySurf.surface(),
         querySurf,
         pointField(1, mesh.cellCentres()[0]),
-        false,
-        false,
-        false,
+        false,              // includeCut
+        false,              // includeInside
+        false,              // includeOutside
+        false,              // geometricOnly
         nearDist,
         curvature
     );
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
index e42c52451f..bc4a363517 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
@@ -64,6 +64,7 @@ initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C
 initialPointsMethod/faceCentredCubic/faceCentredCubic.C
 initialPointsMethod/pointFile/pointFile.C
 initialPointsMethod/autoDensity/autoDensity.C
+initialPointsMethod/rayShooting/rayShooting.C
 
 relaxationModel/relaxationModel/relaxationModel.C
 relaxationModel/adaptiveLinear/adaptiveLinear.C
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 6f137249c1..e580febdd8 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -1249,7 +1249,7 @@ void Foam::conformalVoronoiMesh::move()
                     if
                     (
                         (
-                            (vA->internalPoint() && vB->internalPoint())
+                            (vA->internalPoint() || vB->internalPoint())
                          && (!vA->referred() || !vB->referred())
 //                         ||
 //                            (
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 20a8b278eb..0360cf0d83 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -670,6 +670,32 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
                 continue;
             }
 
+            const searchableSurface& surface(allGeometry_[surfaces_[s]]);
+
+            if (!surface.hasVolumeType())
+            {
+                pointField sample(1, samplePts[i]);
+                scalarField nearestDistSqr(1, GREAT);
+                List<pointIndexHit> info;
+
+                surface.findNearest(sample, nearestDistSqr, info);
+
+                vector hitDir = info[0].rawPoint() - samplePts[i];
+                hitDir /= mag(hitDir) + SMALL;
+
+                if
+                (
+                    findSurfaceAnyIntersection
+                    (
+                        samplePts[i],
+                        info[0].rawPoint() - 1e-3*mag(hitDir)*(hitDir)
+                    )
+                )
+                {
+                    continue;
+                }
+            }
+
             if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
             {
                 if
@@ -694,44 +720,6 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
                     break;
                 }
             }
-//            else
-//            {
-//                // Surface volume type is unknown
-//                Info<< "UNKNOWN" << endl;
-//                // Get nearest face normal
-//
-//                pointField sample(1, samplePts[i]);
-//                scalarField nearestDistSqr(1, GREAT);
-//                List<pointIndexHit> info;
-//                vectorField norms(1);
-//
-//                surface.findNearest(sample, nearestDistSqr, info);
-//                surface.getNormal(info, norms);
-//
-//                vector fN = norms[0];
-//                fN /= mag(fN);
-//
-//                vector hitDir = info[0].rawPoint() - samplePts[i];
-//                hitDir /= mag(hitDir);
-//
-//                if ((fN & hitDir) < 0)
-//                {
-//                    // Point is OUTSIDE
-//
-//                    if
-//                    (
-//                        normalVolumeTypes_[regionI]
-//                     == extendedFeatureEdgeMesh::OUTSIDE
-//                    )
-//                    {
-//                    }
-//                    else
-//                    {
-//                        insidePoint[i] = false;
-//                        break;
-//                    }
-//                }
-//            }
         }
     }
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
new file mode 100644
index 0000000000..ebecad88e1
--- /dev/null
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rayShooting.H"
+#include "addToRunTimeSelectionTable.H"
+#include "triSurfaceMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(rayShooting, 0);
+addToRunTimeSelectionTable(initialPointsMethod, rayShooting, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+rayShooting::rayShooting
+(
+    const dictionary& initialPointsDict,
+    const Time& runTime,
+    Random& rndGen,
+    const conformationSurfaces& geometryToConformTo,
+    const cellShapeControl& cellShapeControls,
+    const autoPtr<backgroundMeshDecomposition>& decomposition
+)
+:
+    initialPointsMethod
+    (
+        typeName,
+        initialPointsDict,
+        runTime,
+        rndGen,
+        geometryToConformTo,
+        cellShapeControls,
+        decomposition
+    ),
+    randomiseInitialGrid_(detailsDict().lookup("randomiseInitialGrid")),
+    randomPerturbationCoeff_
+    (
+        readScalar(detailsDict().lookup("randomPerturbationCoeff"))
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+List<Vb::Point> rayShooting::initialPoints() const
+{
+    // Loop over surface faces
+    const searchableSurfaces& surfaces = geometryToConformTo().geometry();
+    const labelList& surfacesToConformTo = geometryToConformTo().surfaces();
+
+    const scalar maxRayLength = surfaces.bounds().mag();
+
+    // Initialise points list
+    label initialPointsSize = 0;
+    forAll(surfaces, surfI)
+    {
+        initialPointsSize += surfaces[surfI].size();
+    }
+
+    DynamicList<Vb::Point> initialPoints(initialPointsSize);
+
+    forAll(surfacesToConformTo, surfI)
+    {
+        const searchableSurface& s = surfaces[surfacesToConformTo[surfI]];
+
+        tmp<pointField> faceCentresTmp(s.coordinates());
+        const pointField& faceCentres = faceCentresTmp();
+
+        Info<< "    Shoot rays from " << s.name() << nl
+            << "    nRays = " << faceCentres.size() << endl;
+
+
+        forAll(faceCentres, fcI)
+        {
+            const Foam::point& fC = faceCentres[fcI];
+
+            if
+            (
+                Pstream::parRun()
+             && !decomposition().positionOnThisProcessor(fC)
+            )
+            {
+                continue;
+            }
+
+            const scalar pert =
+                randomPerturbationCoeff_
+               *cellShapeControls().cellSize(fC);
+
+            pointIndexHit surfHitStart;
+            label hitSurfaceStart;
+
+            // Face centres should be on the surface so search distance can be
+            // small
+            geometryToConformTo().findSurfaceNearest
+            (
+                 fC,
+                 sqr(pert),
+                 surfHitStart,
+                 hitSurfaceStart
+            );
+
+            vectorField normStart(1, vector::min);
+            geometryToConformTo().getNormal
+            (
+                hitSurfaceStart,
+                List<pointIndexHit>(1, surfHitStart),
+                normStart
+            );
+
+            pointIndexHit surfHitEnd;
+            label hitSurfaceEnd;
+
+            geometryToConformTo().findSurfaceNearestIntersection
+            (
+                fC - normStart[0]*pert,
+                fC - normStart[0]*maxRayLength,
+                surfHitEnd,
+                hitSurfaceEnd
+            );
+
+            if (surfHitEnd.hit())
+            {
+                vectorField normEnd(1, vector::min);
+                geometryToConformTo().getNormal
+                (
+                    hitSurfaceEnd,
+                    List<pointIndexHit>(1, surfHitEnd),
+                    normEnd
+                );
+
+                if ((normStart[0] & normEnd[0]) < 0)
+                {
+                    line<point, point> l(fC, surfHitEnd.hitPoint());
+
+                    if (Pstream::parRun())
+                    {
+                        // Clip the line in parallel
+                        pointIndexHit procIntersection =
+                            decomposition().findLine
+                            (
+                                l.start(),
+                                l.end()
+                            );
+
+                        if (procIntersection.hit())
+                        {
+                            l =
+                                line<point, point>
+                                (
+                                    l.start(),
+                                    procIntersection.hitPoint()
+                                );
+                        }
+                    }
+
+                    Foam::point midPoint(l.centre());
+
+                    const scalar minDistFromSurfaceSqr =
+                        minimumSurfaceDistanceCoeffSqr_
+                       *sqr(cellShapeControls().cellSize(midPoint));
+
+                    if (randomiseInitialGrid_)
+                    {
+                        midPoint.x() += pert*(rndGen().scalar01() - 0.5);
+                        midPoint.y() += pert*(rndGen().scalar01() - 0.5);
+                        midPoint.z() += pert*(rndGen().scalar01() - 0.5);
+                    }
+
+                    if
+                    (
+                        magSqr(midPoint - l.start()) > minDistFromSurfaceSqr
+                     && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr
+                    )
+                    {
+                        initialPoints.append(toPoint(midPoint));
+                    }
+                }
+            }
+        }
+    }
+
+    return initialPoints.shrink();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H
new file mode 100644
index 0000000000..e66df62b7c
--- /dev/null
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H
@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::rayShooting
+
+Description
+
+SourceFiles
+    rayShooting.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rayShooting_H
+#define rayShooting_H
+
+#include "initialPointsMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class rayShooting Declaration
+\*---------------------------------------------------------------------------*/
+
+class rayShooting
+:
+    public initialPointsMethod
+{
+
+private:
+
+    // Private data
+
+        //- Should the initial positions be randomised
+        Switch randomiseInitialGrid_;
+
+        //- Randomise the initial positions by fraction of the initialCellSize_
+        scalar randomPerturbationCoeff_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("rayShooting");
+
+    // Constructors
+
+        //- Construct from components
+        rayShooting
+        (
+            const dictionary& initialPointsDict,
+            const Time& runTime,
+            Random& rndGen,
+            const conformationSurfaces& geometryToConformTo,
+            const cellShapeControl& cellShapeControls,
+            const autoPtr<backgroundMeshDecomposition>& decomposition
+        );
+
+
+    //- Destructor
+    virtual ~rayShooting()
+    {}
+
+
+    // Member Functions
+
+        //- Return the initial points for the conformalVoronoiMesh
+        virtual List<Vb::Point> initialPoints() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSetDict b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
index 11dccf942b..9dee4d1df9 100644
--- a/applications/utilities/mesh/manipulation/topoSet/topoSetDict
+++ b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
@@ -143,6 +143,9 @@ FoamFile
 //    sourceInfo
 //    {
 //        file            "www.avl.com-geometry.stl";
+//        useSurfaceOrientation false;  // use closed surface inside/outside
+//                                      // test (ignores includeCut,
+//                                      // outsidePoints)
 //        outsidePoints   ((-99 -99 -59));    // definition of outside
 //        includeCut      false;              // cells cut by surface
 //        includeInside   false;              // cells not on outside of surf
diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C
index 5a6fad75a6..4d2ab943ac 100644
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.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
@@ -184,6 +184,7 @@ int main(int argc, char *argv[])
             }
         }
 
+
         if (nTracks == 0)
         {
             Info<< "\n    No track data" << endl;
@@ -202,11 +203,13 @@ int main(int argc, char *argv[])
 
             // particle "age" property used to sort the tracks
             List<SortableList<scalar> > agePerTrack(nTracks);
+            List<List<label> > particleMap(nTracks);
 
             forAll(trackLengths, i)
             {
                 const label length = trackLengths[i];
                 agePerTrack[i].setSize(length);
+                particleMap[i].setSize(length);
             }
 
             // store the particle age per track
@@ -228,6 +231,7 @@ int main(int argc, char *argv[])
                     const label trackI = particleToTrack[i];
                     const label sampleI = trackSamples[trackI];
                     agePerTrack[trackI][sampleI] = age[i];
+                    particleMap[trackI][sampleI] = i;
                     trackSamples[trackI]++;
                 }
                 tage.clear();
@@ -251,21 +255,30 @@ int main(int argc, char *argv[])
                 Info<< "\n    Writing points" << endl;
 
                 {
-                    label offset = 0;
                     forAll(agePerTrack, i)
                     {
                         agePerTrack[i].sort();
+
                         const labelList& ids = agePerTrack[i].indices();
+                        labelList& particleIds = particleMap[i];
+
+                        {
+                            // update addressing
+                            List<label> sortedIds(ids);
+                            forAll(sortedIds, j)
+                            {
+                                sortedIds[j] = particleIds[ids[j]];
+                            }
+                            particleIds = sortedIds;
+                        }
 
                         forAll(ids, j)
                         {
-                            const label localId = offset + ids[j];
+                            const label localId = particleIds[j];
                             const vector& pos = particles[localId].position();
                             os  << pos.x() << ' ' << pos.y() << ' ' << pos.z()
                                 << nl;
                         }
-
-                        offset += trackLengths[i];
                     }
                 }
 
@@ -278,13 +291,14 @@ int main(int argc, char *argv[])
                 // Write ids of track points to file
                 {
                     label globalPtI = 0;
-                    forAll(agePerTrack, i)
+                    forAll(particleMap, i)
                     {
-                        os  << agePerTrack[i].size() << nl;
+                        os  << particleMap[i].size() << nl;
 
-                        forAll(agePerTrack[i], j)
+                        forAll(particleMap[i], j)
                         {
                             os  << ' ' << globalPtI++;
+
                             if (((j + 1) % 10 == 0) && (j != 0))
                             {
                                 os  << nl;
@@ -303,14 +317,14 @@ int main(int argc, char *argv[])
 
                 Info<< "\n    Processing fields" << nl << endl;
 
-                processFields<label>(os, agePerTrack, userFields, cloudObjs);
-                processFields<scalar>(os, agePerTrack, userFields, cloudObjs);
-                processFields<vector>(os, agePerTrack, userFields, cloudObjs);
+                processFields<label>(os, particleMap, userFields, cloudObjs);
+                processFields<scalar>(os, particleMap, userFields, cloudObjs);
+                processFields<vector>(os, particleMap, userFields, cloudObjs);
                 processFields<sphericalTensor>
-                    (os, agePerTrack, userFields, cloudObjs);
+                    (os, particleMap, userFields, cloudObjs);
                 processFields<symmTensor>
-                    (os, agePerTrack, userFields, cloudObjs);
-                processFields<tensor>(os, agePerTrack, userFields, cloudObjs);
+                    (os, particleMap, userFields, cloudObjs);
+                processFields<tensor>(os, particleMap, userFields, cloudObjs);
 
             }
         }
diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C
index 89e499b17d..3dcc365c1d 100644
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.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
@@ -57,39 +57,57 @@ tmp<Field<Type> > readParticleField
         return tmp<Field<Type> >(new Field<Type>(newField.xfer()));
     }
 
-    Info<< "error: cloud field name " << name << " not found" << endl;
+    FatalErrorIn
+    (
+        "template<class Type>"
+        "void readParticleField"
+        "("
+            "const word&, "
+            "const IOobjectList"
+        ")"
+    )
+        << "error: cloud field name " << name << " not found"
+        << abort(FatalError);
 
     return Field<Type>::null();
 }
 
 
 template<class Type>
-PtrList<List<Type> > readFields
+void readFields
 (
     PtrList<List<Type> >& values,
-    const List<word>& fields,
+    const List<word>& fieldNames,
     const IOobjectList& cloudObjs
 )
 {
     IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
 
-    label fieldI = 0;
-    forAllConstIter(IOobjectList, objects, iter)
+    forAll(fieldNames, j)
     {
-        const IOobject& obj = *iter();
-        forAll(fields, j)
+        const IOobject* obj = objects.lookup(fieldNames[j]);
+        if (obj != NULL)
         {
-            if (obj.name() == fields[j])
-            {
-                Info<< "        reading field " << obj.name() << endl;
-                IOField<Type> newField(obj);
-                values.set(fieldI++, new List<Type>(newField.xfer()));
-                break;
-            }
+            Info<< "        reading field " << fieldNames[j] << endl;
+            IOField<Type> newField(*obj);
+            values.set(j, new List<Type>(newField.xfer()));
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "template<class Type>"
+                "void readFields"
+                "("
+                    "PtrList<List<Type> >&, "
+                    "const List<word>&, "
+                    "const IOobjectList&"
+                ")"
+            )
+                << "Unable to read field " << fieldNames[j]
+                << abort(FatalError);
         }
     }
-
-    return values;
 }
 
 
@@ -109,7 +127,7 @@ void writeVTKFields
 (
     OFstream& os,
     const PtrList<List<Type> >& values,
-    const List<SortableList<scalar> >& agePerTrack,
+    const List<List<label> >& addr,
     const List<word>& fieldNames
 )
 {
@@ -121,9 +139,9 @@ void writeVTKFields
         os  << nl << fieldNames[fieldI] << ' ' << pTraits<Type>::nComponents
             << ' ' << values[fieldI].size() << " float" << nl;
         label offset = 0;
-        forAll(agePerTrack, trackI)
+        forAll(addr, trackI)
         {
-            const List<label> ids = agePerTrack[trackI].indices() + offset;
+            const List<label> ids(addr[trackI]);
 
             List<Type> data(UIndirectList<Type>(values[fieldI], ids));
             label nData = data.size() - 1;
@@ -149,7 +167,7 @@ template<class Type>
 void processFields
 (
     OFstream& os,
-    const List<SortableList<scalar> >& agePerTrack,
+    const List<List<label> >& addr,
     const List<word>& userFieldNames,
     const IOobjectList& cloudObjs
 )
@@ -176,12 +194,14 @@ void processFields
         (
             os,
             values,
-            agePerTrack,
-            fieldNames.xfer()
+            addr,
+            fieldNames
         );
     }
 }
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 } // End namespace Foam
 
 // ************************************************************************* //
diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.H b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.H
index 973e6da255..f3b5d48d2d 100644
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.H
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.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
@@ -49,7 +49,7 @@ namespace Foam
     );
 
     template<class Type>
-    PtrList<List<Type> > readFields
+    void readFields
     (
         PtrList<List<Type> >& values,
         const List<word>& fields,
diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict
index 86ba1b4849..60e1d7d0ee 100644
--- a/applications/utilities/postProcessing/sampling/sample/sampleDict
+++ b/applications/utilities/postProcessing/sampling/sample/sampleDict
@@ -165,6 +165,10 @@ surfaces
 
         //- Optional: restrict to a particular zone
         // zone        zone1;
+
+        //- Optional: do not triangulate (only for surfaceFormats that support
+        //            polygons)
+        //triangulate     false;
     }
 
     interpolatedPlane
diff --git a/applications/utilities/postProcessing/turbulence/R/Make/options b/applications/utilities/postProcessing/turbulence/R/Make/options
index a136b16661..27b70cae0a 100644
--- a/applications/utilities/postProcessing/turbulence/R/Make/options
+++ b/applications/utilities/postProcessing/turbulence/R/Make/options
@@ -1,13 +1,19 @@
 EXE_INC = \
-    -I$(LIB_SRC)/postProcessing/postCalc \
-    -I$(LIB_SRC)/turbulenceModels \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
-    $(FOAM_LIBBIN)/postCalc.o \
-    -lincompressibleRASModels \
     -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lcompressibleRASModels \
     -lfiniteVolume \
-    -lgenericPatchFields
+    -lgenericPatchFields \
+    -lmeshTools \
+    -lsampling
diff --git a/applications/utilities/postProcessing/turbulence/R/R.C b/applications/utilities/postProcessing/turbulence/R/R.C
index b96bb03947..708cda7214 100644
--- a/applications/utilities/postProcessing/turbulence/R/R.C
+++ b/applications/utilities/postProcessing/turbulence/R/R.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
@@ -29,35 +29,140 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "calc.H"
 #include "fvCFD.H"
 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
-#include "RASModel.H"
+#include "incompressible/turbulenceModel/turbulenceModel.H"
+
+#include "fluidThermo.H"
+#include "compressible/turbulenceModel/turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
+void calcIncompressibleR
+(
+    const fvMesh& mesh,
+    const Time& runTime,
+    const volVectorField& U
+)
 {
-    #include "createFields.H"
+    #include "createPhi.H"
 
-    Info<< "\nCalculating the Reynolds Stress R\n" << endl;
+    singlePhaseTransportModel laminarTransport(U, phi);
 
-    volSymmTensorField R
+    autoPtr<incompressible::turbulenceModel> model
     (
-        IOobject
-        (
-            "R",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        RASModel->R()
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
     );
 
-    R.write();
+    Info<< "Writing R field" << nl << endl;
 
-    Info<< "End" << endl;
+    model->R()().write();
+}
+
+
+void calcCompressibleR
+(
+    const fvMesh& mesh,
+    const Time& runTime,
+    const volVectorField& U
+)
+{
+    IOobject rhoHeader
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (!rhoHeader.headerOk())
+    {
+        Info<< "    no " << rhoHeader.name() <<" field" << endl;
+        return;
+    }
+
+    Info<< "Reading field rho\n" << endl;
+    volScalarField rho(rhoHeader, mesh);
+
+    #include "compressibleCreatePhi.H"
+
+    autoPtr<fluidThermo> pThermo(fluidThermo::New(mesh));
+    fluidThermo& thermo = pThermo();
+
+    autoPtr<compressible::turbulenceModel> model
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    Info<< "Writing R field" << nl << endl;
+
+    model->R()().write();
+}
+
+
+int main(int argc, char *argv[])
+{
+    timeSelector::addOptions();
+
+    #include "addRegionOption.H"
+
+    argList::addBoolOption
+    (
+        "compressible",
+        "calculate compressible R"
+    );
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    instantList timeDirs = timeSelector::select0(runTime, args);
+    #include "createNamedMesh.H"
+
+    const bool compressible = args.optionFound("compressible");
+
+    forAll(timeDirs, timeI)
+    {
+        runTime.setTime(timeDirs[timeI], timeI);
+        Info<< "Time = " << runTime.timeName() << endl;
+
+        IOobject UHeader
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+
+        if (UHeader.headerOk())
+        {
+            Info<< "Reading field " << UHeader.name() << nl << endl;
+            volVectorField U(UHeader, mesh);
+
+            if (compressible)
+            {
+                calcCompressibleR(mesh, runTime, U);
+            }
+            else
+            {
+                calcIncompressibleR(mesh, runTime, U);
+            }
+        }
+        else
+        {
+            Info<< "    no " << UHeader.name() << " field" << endl;
+        }
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
 }
 
 
diff --git a/applications/utilities/postProcessing/turbulence/R/createFields.H b/applications/utilities/postProcessing/turbulence/R/createFields.H
deleted file mode 100644
index 7c07f44f8b..0000000000
--- a/applications/utilities/postProcessing/turbulence/R/createFields.H
+++ /dev/null
@@ -1,22 +0,0 @@
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    #include "createPhi.H"
-
-    singlePhaseTransportModel laminarTransport(U, phi);
-
-    autoPtr<incompressible::RASModel> RASModel
-    (
-        incompressible::RASModel::New(U, phi, laminarTransport)
-    );
diff --git a/applications/utilities/preProcessing/changeDictionary/changeDictionary.C b/applications/utilities/preProcessing/changeDictionary/changeDictionary.C
index a53a6142cf..a6b18b6114 100644
--- a/applications/utilities/preProcessing/changeDictionary/changeDictionary.C
+++ b/applications/utilities/preProcessing/changeDictionary/changeDictionary.C
@@ -602,8 +602,13 @@ int main(int argc, char *argv[])
                 dictList.set(sz++, iter().clone());
             }
 
-            Info<< "Writing modified fieldDict " << fieldName << endl;
-            dictList.write();
+            Info<< "Writing modified " << fieldName << endl;
+            dictList.writeObject
+            (
+                runTime.writeFormat(),
+                runTime.writeFormat(),
+                IOstream::UNCOMPRESSED
+            );
         }
         else
         {
diff --git a/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C b/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C
index 66065c4528..1c774c352b 100644
--- a/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C
+++ b/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C
@@ -39,6 +39,14 @@ namespace functionEntries
 {
     defineTypeNameAndDebug(codeStream, 0);
 
+    addToMemberFunctionSelectionTable
+    (
+        functionEntry,
+        codeStream,
+        execute,
+        dictionaryIstream
+    );
+
     addToMemberFunctionSelectionTable
     (
         functionEntry,
@@ -364,6 +372,38 @@ bool Foam::functionEntries::codeStream::execute
     IStringStream resultStream(os.str());
     entry.read(parentDict, resultStream);
 
+    return true;
+}
+
+
+bool Foam::functionEntries::codeStream::execute
+(
+    dictionary& parentDict,
+    Istream& is
+)
+{
+    Info<< "Using #codeStream at line " << is.lineNumber()
+        << " in file " <<  parentDict.name() << endl;
+
+    dynamicCode::checkSecurity
+    (
+        "functionEntries::codeStream::execute(..)",
+        parentDict
+    );
+
+    // get code dictionary
+    // must reference parent for stringOps::expand to work nicely
+    dictionary codeDict("#codeStream", parentDict, is);
+
+    streamingFunctionType function = getFunction(parentDict, codeDict);
+
+    // use function to write stream
+    OStringStream os(is.format());
+    (*function)(os, parentDict);
+
+    // get the entry from this stream
+    IStringStream resultStream(os.str());
+    parentDict.read(resultStream);
 
     return true;
 }
diff --git a/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.H b/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.H
index b8baa0183c..1dc2b57c87 100644
--- a/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.H
+++ b/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.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
@@ -160,6 +160,10 @@ public:
 
     // Member Functions
 
+        //- Execute the functionEntry in a sub-dict context
+        static bool execute(dictionary& parentDict, Istream&);
+
+        //- Execute the functionEntry in a primitiveEntry context
         static bool execute
         (
             const dictionary& parentDict,
diff --git a/src/OpenFOAM/primitives/functions/DataEntry/CSV/CSVIO.C b/src/OpenFOAM/primitives/functions/DataEntry/CSV/CSVIO.C
index 1dd395cee0..778de9a92f 100644
--- a/src/OpenFOAM/primitives/functions/DataEntry/CSV/CSVIO.C
+++ b/src/OpenFOAM/primitives/functions/DataEntry/CSV/CSVIO.C
@@ -50,10 +50,7 @@ Foam::Ostream& Foam::operator<<
     }
 
     // Check state of Ostream
-    os.check
-    (
-        "Ostream& operator<<(Ostream&, const CSV<Type>&)"
-    );
+    os.check("Ostream& operator<<(Ostream&, const CSV<Type>&)");
 
     return os;
 }
@@ -83,11 +80,15 @@ void Foam::CSV<Type>::writeData(Ostream& os) const
         os  << componentColumns_;
         os.format(IOstream::BINARY);
     }
+    else
+    {
+        os  << componentColumns_;
+    }
     os  << token::END_STATEMENT << nl;
 
     os.writeKeyword("separator") << string(separator_)
         << token::END_STATEMENT << nl;
-    os.writeKeyword("mergeSeparators") << string(mergeSeparators_)
+    os.writeKeyword("mergeSeparators") << mergeSeparators_
         << token::END_STATEMENT << nl;
     os.writeKeyword("fileName") << fName_ << token::END_STATEMENT << nl;
     os  << decrIndent << indent << token::END_BLOCK << endl;
diff --git a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C
index b712a51c7f..b052b4de86 100644
--- a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C
+++ b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.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
@@ -96,6 +96,7 @@ void Foam::layerAdditionRemoval::checkDefinition()
     }
 }
 
+
 Foam::scalar Foam::layerAdditionRemoval::readOldThickness
 (
     const dictionary& dict
@@ -115,7 +116,6 @@ void Foam::layerAdditionRemoval::clearAddressing() const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::layerAdditionRemoval::layerAdditionRemoval
 (
     const word& name,
@@ -140,7 +140,6 @@ Foam::layerAdditionRemoval::layerAdditionRemoval
 }
 
 
-// Construct from dictionary
 Foam::layerAdditionRemoval::layerAdditionRemoval
 (
     const word& name,
@@ -377,7 +376,7 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
         // Clear addressing.  This also resets the addition/removal data
         if (debug)
         {
-            Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
+            Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
                 << " for object " << name() << " : "
                 << "Clearing addressing after layer removal. " << endl;
         }
@@ -393,7 +392,7 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
         // Clear addressing.  This also resets the addition/removal data
         if (debug)
         {
-            Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
+            Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
                 << " for object " << name() << " : "
                 << "Clearing addressing after layer addition. " << endl;
         }
@@ -431,16 +430,14 @@ void Foam::layerAdditionRemoval::updateMesh(const mapPolyMesh&)
 
 void Foam::layerAdditionRemoval::setMinLayerThickness(const scalar t) const
 {
-    if
-    (
-        t < VSMALL
-     || maxLayerThickness_ < t
-    )
+    if (t < VSMALL || maxLayerThickness_ < t)
     {
         FatalErrorIn
         (
-            "void layerAdditionRemoval::setMinLayerThickness("
-            "const scalar t) const"
+            "void layerAdditionRemoval::setMinLayerThickness"
+            "("
+                "const scalar"
+            ") const"
         )   << "Incorrect layer thickness definition."
             << abort(FatalError);
     }
@@ -451,15 +448,14 @@ void Foam::layerAdditionRemoval::setMinLayerThickness(const scalar t) const
 
 void Foam::layerAdditionRemoval::setMaxLayerThickness(const scalar t) const
 {
-    if
-    (
-        t < minLayerThickness_
-    )
+    if (t < minLayerThickness_)
     {
         FatalErrorIn
         (
-            "void layerAdditionRemoval::setMaxLayerThickness("
-            "const scalar t) const"
+            "void layerAdditionRemoval::setMaxLayerThickness"
+            "("
+                "const scalar"
+            ") const"
         )   << "Incorrect layer thickness definition."
             << abort(FatalError);
     }
diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index 8ffb988d45..c71559d042 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.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
@@ -34,6 +34,8 @@ const Foam::scalar Foam::particle::trackingCorrectionTol = 1e-5;
 
 const Foam::scalar Foam::particle::lambdaDistanceToleranceCoeff = 1e3*SMALL;
 
+const Foam::scalar Foam::particle::minStepFractionTol = 1e5*SMALL;
+
 namespace Foam
 {
     defineTypeNameAndDebug(particle, 0);
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index 15742cacc9..0275cc5fd1 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -307,6 +307,9 @@ public:
         //  for the denominator and numerator of lambda
         static const scalar lambdaDistanceToleranceCoeff;
 
+        //- Minimum stepFraction tolerance
+        static const scalar minStepFractionTol;
+
 
     // Constructors
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index acf9e8c235..7cab777680 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -652,7 +652,11 @@ void Foam::KinematicCloud<CloudType>::scaleSources()
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::preEvolve()
 {
-    Info<< "\nSolving cloud " << this->name() << endl;
+    // force calculaion of mesh dimensions - needed for parallel runs
+    // with topology change due to lazy evaluation of valid mesh dimensions
+    label nGeometricD = mesh_.nGeometricD();
+
+    Info<< "\nSolving " << nGeometricD << "-D cloud " << this->name() << endl;
 
     this->dispersion().cacheFields(true);
     forces_.cacheFields(true);
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
index fd7bf31706..c2bea38236 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.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
@@ -110,11 +110,11 @@ void Foam::cloudSolution::read()
     dict_.lookup("transient") >> transient_;
     dict_.lookup("coupled") >> coupled_;
     dict_.lookup("cellValueSourceCorrection") >> cellValueSourceCorrection_;
+    dict_.lookup("maxCo") >> maxCo_;
 
     if (steadyState())
     {
         dict_.lookup("calcFrequency") >> calcFrequency_;
-        dict_.lookup("maxCo") >> maxCo_;
         dict_.lookup("maxTrackTime") >> maxTrackTime_;
 
         if (coupled_)
diff --git a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.C b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.C
index f04ac77a15..f422f7ea37 100644
--- a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.C
@@ -69,13 +69,6 @@ bool Foam::CollidingParcel<ParcelType>::move
     typename TrackData::cloudType::parcelType& p =
         static_cast<typename TrackData::cloudType::parcelType&>(*this);
 
-    td.switchProcessor = false;
-    td.keepParticle = true;
-
-    const polyMesh& mesh = td.cloud().pMesh();
-    const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
-    const scalarField& V = mesh.cellVolumes();
-
     switch (td.part())
     {
         case TrackData::tpVelocityHalfStep:
@@ -92,70 +85,7 @@ bool Foam::CollidingParcel<ParcelType>::move
 
         case TrackData::tpLinearTrack:
         {
-            scalar tEnd = (1.0 - p.stepFraction())*trackTime;
-            const scalar dtMax = tEnd;
-
-            while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
-            {
-                // Apply correction to position for reduced-D cases
-                meshTools::constrainToMeshCentre(mesh, p.position());
-
-                const point start(p.position());
-
-                // Set the Lagrangian time-step
-                scalar dt = min(dtMax, tEnd);
-
-                // Remember which cell the parcel is in since this
-                // will change if a face is hit
-                const label cellI = p.cell();
-
-                const scalar magU = mag(p.U());
-                if (p.active() && magU > ROOTVSMALL)
-                {
-                    const scalar d = dt*magU;
-                    const scalar maxCo = td.cloud().solution().maxCo();
-                    const scalar dCorr = min(d, maxCo*cbrt(V[cellI]));
-                    dt *=
-                        dCorr/d
-                       *p.trackToFace(p.position() + dCorr*p.U()/magU, td);
-                }
-
-                tEnd -= dt;
-                p.stepFraction() = 1.0 - tEnd/trackTime;
-
-                // Avoid problems with extremely small timesteps
-                if (dt > ROOTVSMALL)
-                {
-                    // Update cell based properties
-                    p.setCellValues(td, dt, cellI);
-
-                    if (td.cloud().solution().cellValueSourceCorrection())
-                    {
-                        p.cellValueSourceCorrection(td, dt, cellI);
-                    }
-
-                    p.calc(td, dt, cellI);
-                }
-
-                if (p.onBoundary() && td.keepParticle)
-                {
-                    if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
-                    {
-                        td.switchProcessor = true;
-                    }
-                }
-
-                p.age() += dt;
-
-                td.cloud().functions().postMove
-                (
-                    p,
-                    cellI,
-                    dt,
-                    start,
-                    td.keepParticle
-                );
-            }
+            ParcelType::move(td, trackTime);
 
             break;
         }
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 8cf9317c86..b833962825 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -271,6 +271,8 @@ bool Foam::KinematicParcel<ParcelType>::move
     scalar tEnd = (1.0 - p.stepFraction())*trackTime;
     const scalar dtMax = tEnd;
 
+    bool moving = true;
+
     while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
     {
         // Apply correction to position for reduced-D cases
@@ -281,12 +283,11 @@ bool Foam::KinematicParcel<ParcelType>::move
         // Set the Lagrangian time-step
         scalar dt = min(dtMax, tEnd);
 
-        // Remember which cell the parcel is in since this will change if
-        // a face is hit
+        // Cache the parcel current cell as this will change if a face is hit
         const label cellI = p.cell();
 
         const scalar magU = mag(U_);
-        if (p.active() && magU > ROOTVSMALL)
+        if (p.active() && moving && (magU > ROOTVSMALL))
         {
             const scalar d = dt*magU;
             const scalar dCorr = min(d, maxCo*cbrt(V[cellI]));
@@ -296,7 +297,19 @@ bool Foam::KinematicParcel<ParcelType>::move
         }
 
         tEnd -= dt;
-        p.stepFraction() = 1.0 - tEnd/trackTime;
+
+        scalar newStepFraction = 1.0 - tEnd/trackTime;
+
+        if
+        (
+            mag(p.stepFraction() - newStepFraction)
+          < particle::minStepFractionTol
+        )
+        {
+            moving = false;
+        }
+
+        p.stepFraction() = newStepFraction;
 
         // Avoid problems with extremely small timesteps
         if (dt > ROOTVSMALL)
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
index 68d19890dd..a8003d43be 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
@@ -64,31 +64,40 @@ void Foam::ParticleCollector<CloudType>::makeLogFile
 
             outputFilePtr_()
                 << "# Source     : " << type() << nl
+                << "# Bins       : " << faces.size() << nl
                 << "# Total area : " << sum(area) << nl;
 
-            outputFilePtr_() << "# Geometry   :" << nl;
+            outputFilePtr_()
+                << "# Geometry   :" << nl
+                << '#'
+                << tab << "Bin"
+                << tab << "(Centre_x Centre_y Centre_z)"
+                << tab << "Area"
+                << nl;
+
             forAll(faces, i)
             {
-                word id = Foam::name(i);
-
                 outputFilePtr_()
-                    << '#' << tab << "point[" << id << "] = "
-                    << faces[i].centre(points) << nl;
+                    << '#'
+                    << tab << i
+                    << tab << faces[i].centre(points)
+                    << tab << area[i]
+                    << nl;
             }
 
-            outputFilePtr_()<< '#' << nl << "# Time";
+            outputFilePtr_()
+                << '#' << nl
+                << "# Output format:" << nl;
 
             forAll(faces, i)
             {
                 word id = Foam::name(i);
-
-                if (i != 0)
-                {
-                    outputFilePtr_() << "#";
-                }
+                word binId = "bin_" + id;
 
                 outputFilePtr_()
-                    << tab << "area[" << id << "]"
+                    << '#'
+                    << tab << "Time"
+                    << tab << binId
                     << tab << "mass[" << id << "]"
                     << tab << "massFlowRate[" << id << "]"
                     << endl;
@@ -406,12 +415,6 @@ void Foam::ParticleCollector<CloudType>::write()
 
     Info<< type() << " output:" << nl;
 
-    if (outputFilePtr_.valid())
-    {
-        outputFilePtr_() << time.timeName();
-    }
-
-
     Field<scalar> faceMassTotal(mass_.size(), 0.0);
     this->getModelProperty("massTotal", faceMassTotal);
 
@@ -439,7 +442,8 @@ void Foam::ParticleCollector<CloudType>::write()
         if (outputFilePtr_.valid())
         {
             outputFilePtr_()
-                << tab << area_[faceI]
+                << time.timeName()
+                << tab << faceI
                 << tab << faceMassTotal[faceI]
                 << tab << faceMassFlowRate[faceI]
                 << endl;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
index d9b38daccf..fa86f67ca9 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.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
@@ -57,7 +57,7 @@ Foam::StandardWallInteraction<CloudType>::StandardWallInteraction
                 "StandardWallInteraction<CloudType>::StandardWallInteraction"
                 "("
                     "const dictionary&, "
-                    "CloudType& cloud"
+                    "CloudType&"
                 ")"
             )   << "Unknown interaction result type "
                 << interactionTypeName
@@ -174,10 +174,11 @@ bool Foam::StandardWallInteraction<CloudType>::correct
                 (
                     "bool StandardWallInteraction<CloudType>::correct"
                     "("
+                        "typename CloudType::parcelType&, "
                         "const polyPatch&, "
-                        "const label, "
-                        "bool&, "
-                        "vector&"
+                        "bool& keepParticle, "
+                        "const scalar, "
+                        "const tetIndices&"
                     ") const"
                 )   << "Unknown interaction type "
                     << this->interactionTypeToWord(interactionType_)
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.C b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.C
index abfdc581ac..870dc31ffc 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.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
@@ -59,7 +59,6 @@ Foam::HeatTransferModel<CloudType>::HeatTransferModel
 {}
 
 
-
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class CloudType>
diff --git a/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C b/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
index ff37610fad..f4a0009028 100644
--- a/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
+++ b/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
@@ -157,7 +157,30 @@ void Foam::surfaceToCell::combine(topoSet& set, const bool add) const
 {
     cpuTime timer;
 
-    if (includeCut_ || includeInside_ || includeOutside_)
+
+    if (useSurfaceOrientation_ && (includeInside_ || includeOutside_))
+    {
+        const meshSearch queryMesh(mesh_);
+
+        //- Calculate for each searchPoint inside/outside status.
+        boolList isInside(querySurf().calcInside(mesh_.cellCentres()));
+
+        Info<< "    Marked inside/outside using surface orientation in = "
+            << timer.cpuTimeIncrement() << " s" << endl << endl;
+
+        forAll(isInside, cellI)
+        {
+            if (isInside[cellI] && includeInside_)
+            {
+                addOrDelete(set, cellI, add);
+            }
+            else if (!isInside[cellI] && includeOutside_)
+            {
+                addOrDelete(set, cellI, add);
+            }
+        }
+    }
+    else if (includeCut_ || includeInside_ || includeOutside_)
     {
         //
         // Cut cells with surface and classify cells
@@ -166,7 +189,7 @@ void Foam::surfaceToCell::combine(topoSet& set, const bool add) const
 
         // Construct search engine on mesh
 
-        meshSearch queryMesh(mesh_);
+        const meshSearch queryMesh(mesh_);
 
 
         // Check all 'outside' points
@@ -196,10 +219,10 @@ void Foam::surfaceToCell::combine(topoSet& set, const bool add) const
         );
 
 
-        Info<< "    Marked inside/outside in = "
+        Info<< "    Marked inside/outside using surface intersection in = "
             << timer.cpuTimeIncrement() << " s" << endl << endl;
 
-
+        //- Add/remove cells using set
         forAll(cellType, cellI)
         {
             label cType = cellType[cellI];
@@ -326,6 +349,18 @@ void Foam::surfaceToCell::checkSettings() const
             << " or set curvature to a value -1 .. 1."
             << exit(FatalError);
     }
+
+    if (useSurfaceOrientation_ && includeCut_)
+    {
+        FatalErrorIn
+        (
+            "surfaceToCell:checkSettings()"
+        )   << "Illegal include cell specification."
+            << " You cannot specify both 'useSurfaceOrientation'"
+            << " and 'includeCut'"
+            << " since 'includeCut' specifies a topological split"
+            << exit(FatalError);
+    }
 }
 
 
@@ -340,6 +375,7 @@ Foam::surfaceToCell::surfaceToCell
     const bool includeCut,
     const bool includeInside,
     const bool includeOutside,
+    const bool useSurfaceOrientation,
     const scalar nearDist,
     const scalar curvature
 )
@@ -350,6 +386,7 @@ Foam::surfaceToCell::surfaceToCell
     includeCut_(includeCut),
     includeInside_(includeInside),
     includeOutside_(includeOutside),
+    useSurfaceOrientation_(useSurfaceOrientation),
     nearDist_(nearDist),
     curvature_(curvature),
     surfPtr_(new triSurface(surfName_)),
@@ -371,6 +408,7 @@ Foam::surfaceToCell::surfaceToCell
     const bool includeCut,
     const bool includeInside,
     const bool includeOutside,
+    const bool useSurfaceOrientation,
     const scalar nearDist,
     const scalar curvature
 )
@@ -381,6 +419,7 @@ Foam::surfaceToCell::surfaceToCell
     includeCut_(includeCut),
     includeInside_(includeInside),
     includeOutside_(includeOutside),
+    useSurfaceOrientation_(useSurfaceOrientation),
     nearDist_(nearDist),
     curvature_(curvature),
     surfPtr_(&surf),
@@ -404,6 +443,10 @@ Foam::surfaceToCell::surfaceToCell
     includeCut_(readBool(dict.lookup("includeCut"))),
     includeInside_(readBool(dict.lookup("includeInside"))),
     includeOutside_(readBool(dict.lookup("includeOutside"))),
+    useSurfaceOrientation_
+    (
+        dict.lookupOrDefault<bool>("useSurfaceOrientation", false)
+    ),
     nearDist_(readScalar(dict.lookup("nearDistance"))),
     curvature_(readScalar(dict.lookup("curvature"))),
     surfPtr_(new triSurface(surfName_)),
@@ -427,6 +470,7 @@ Foam::surfaceToCell::surfaceToCell
     includeCut_(readBool(checkIs(is))),
     includeInside_(readBool(checkIs(is))),
     includeOutside_(readBool(checkIs(is))),
+    useSurfaceOrientation_(false),
     nearDist_(readScalar(checkIs(is))),
     curvature_(readScalar(checkIs(is))),
     surfPtr_(new triSurface(surfName_)),
diff --git a/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.H b/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.H
index cfaf80dd30..931f412a2a 100644
--- a/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.H
+++ b/src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.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
@@ -29,6 +29,8 @@ Description
 
     Selects:
     - all cells inside/outside/cut by surface
+    - all cells inside/outside surface ('useSurfaceOrientation', requires closed
+      surface)
     - cells with centre nearer than XXX to surface
     - cells with centre nearer than XXX to surface \b and with normal
       at nearest point to centre and cell-corners differing by
@@ -67,27 +69,31 @@ class surfaceToCell
         static addToUsageTable usage_;
 
         //- Name of surface file
-        fileName surfName_;
+        const fileName surfName_;
 
         //- Points which are outside
-        pointField outsidePoints_;
+        const pointField outsidePoints_;
 
         //- Include cut cells
-        bool includeCut_;
+        const bool includeCut_;
 
         //- Include inside cells
-        bool includeInside_;
+        const bool includeInside_;
 
         //- Include outside cells
-        bool includeOutside_;
+        const bool includeOutside_;
+
+        //- Determine inside/outside purely using geometric test
+        //  (does not allow includeCut)
+        const bool useSurfaceOrientation_;
 
         //- if > 0 : include cells with distance from cellCentre to surface
         //  less than nearDist.
-        scalar nearDist_;
+        const scalar nearDist_;
 
         //- if > -1 : include cells with normals at nearest surface points
         //  varying more than curvature_.
-        scalar curvature_;
+        const scalar curvature_;
 
         //- triSurface to search on. On pointer since can be external.
         const triSurface* surfPtr_;
@@ -97,7 +103,7 @@ class surfaceToCell
 
         //- whether I allocated above surface ptrs or whether they are
         //  external.
-        bool IOwnPtrs_;
+        const bool IOwnPtrs_;
 
 
     // Private Member Functions
@@ -155,6 +161,7 @@ public:
             const bool includeCut,
             const bool includeInside,
             const bool includeOutside,
+            const bool useSurfaceOrientation,
             const scalar nearDist,
             const scalar curvature
         );
@@ -170,6 +177,7 @@ public:
             const bool includeCut,
             const bool includeInside,
             const bool includeOutside,
+            const bool useSurfaceOrientation,
             const scalar nearDist,
             const scalar curvature
         );
diff --git a/src/meshTools/sets/topoSets/topoSet.C b/src/meshTools/sets/topoSets/topoSet.C
index 37fe9a0595..691339102a 100644
--- a/src/meshTools/sets/topoSets/topoSet.C
+++ b/src/meshTools/sets/topoSets/topoSet.C
@@ -472,7 +472,6 @@ void Foam::topoSet::invert(const label maxLen)
             insert(cellI);
         }
     }
-
 }
 
 
@@ -550,20 +549,6 @@ void Foam::topoSet::writeDebug(Ostream& os, const label maxLen) const
 }
 
 
-//void Foam::topoSet::writeDebug
-//(
-//    Ostream&,
-//    const primitiveMesh&,
-//    const label
-//) const
-//{
-//    notImplemented
-//    (
-//        "topoSet::writeDebug(Ostream&, const primitiveMesh&, const label)"
-//    );
-//}
-
-
 bool Foam::topoSet::writeData(Ostream& os) const
 {
     return (os << *this).good();
@@ -576,14 +561,6 @@ void Foam::topoSet::updateMesh(const mapPolyMesh&)
 }
 
 
-////- Return max index+1.
-//label topoSet::maxSize(const polyMesh&) const
-//{
-//    notImplemented("topoSet::maxSize(const polyMesh&)");
-//
-//    return -1;
-//}
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 void Foam::topoSet::operator=(const topoSet& rhs)
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
index 10c8a4b4ac..1e6a30627e 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
@@ -32,6 +32,7 @@ License
 #include "fvcVolumeIntegrate.H"
 #include "fvMatrices.H"
 #include "absorptionEmissionModel.H"
+#include "fvcLaplacian.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,6 +62,9 @@ void reactingOneDim::readReactingOneDimControls()
 
     coeffs().lookup("gasHSource") >> gasHSource_;
     coeffs().lookup("QrHSource") >> QrHSource_;
+    useChemistrySolvers_ =
+        coeffs().lookupOrDefault<bool>("useChemistrySolvers", true);
+
 }
 
 
@@ -321,6 +325,8 @@ void reactingOneDim::solveEnergy()
     (
         fvm::ddt(rho_, h_)
       - fvm::laplacian(alpha, h_)
+      + fvc::laplacian(alpha, h_)
+      - fvc::laplacian(kappa(), T())
      ==
         chemistrySh_
       - fvm::Sp(solidChemistry_->RRg(), h_)
@@ -462,7 +468,8 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
     totalGasMassFlux_(0.0),
     totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0)),
     gasHSource_(false),
-    QrHSource_(false)
+    QrHSource_(false),
+    useChemistrySolvers_(true)
 {
     if (active_)
     {
@@ -560,7 +567,8 @@ reactingOneDim::reactingOneDim
     totalGasMassFlux_(0.0),
     totalHeatRR_(dimensionedScalar("zero", dimEnergy/dimTime, 0.0)),
     gasHSource_(false),
-    QrHSource_(false)
+    QrHSource_(false),
+    useChemistrySolvers_(true)
 {
     if (active_)
     {
@@ -681,11 +689,18 @@ void reactingOneDim::evolveRegion()
 {
     Info<< "\nEvolving pyrolysis in region: " << regionMesh().name() << endl;
 
-    solidChemistry_->solve
-    (
-        time().value() - time().deltaTValue(),
-        time().deltaTValue()
-    );
+    if (useChemistrySolvers_)
+    {
+        solidChemistry_->solve
+        (
+            time().value() - time().deltaTValue(),
+            time().deltaTValue()
+        );
+    }
+    else
+    {
+        solidChemistry_->calculate();
+    }
 
     solveContinuity();
 
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
index 4a92fd6967..a5950a9804 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
@@ -154,6 +154,9 @@ protected:
             //- Add in depth radiation source term
             bool  QrHSource_;
 
+            //- Use chemistry solvers (ode or sequential)
+            bool useChemistrySolvers_;
+
 
     // Protected member functions
 
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
index 2f99a1872d..50d1c7c1ef 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
@@ -1049,6 +1049,7 @@ void kinematicSingleLayer::info() const
 {
     Info<< "\nSurface film: " << type() << endl;
 
+    const scalarField& deltaInternal = delta_.internalField();
     const vectorField& Uinternal = U_.internalField();
 
     Info<< indent << "added mass         = "
@@ -1057,8 +1058,8 @@ void kinematicSingleLayer::info() const
         << gSum((deltaRho_*magSf())()) << nl
         << indent << "min/max(mag(U))    = " << gMin(mag(Uinternal)) << ", "
         << gMax(mag(Uinternal)) << nl
-        << indent << "min/max(delta)     = " << min(delta_).value() << ", "
-        << max(delta_).value() << nl
+        << indent << "min/max(delta)     = " << gMin(deltaInternal) << ", "
+        << gMax(deltaInternal) << nl
         << indent << "coverage           = "
         << gSum(alpha_.internalField()*magSf())/gSum(magSf()) <<  nl;
 
diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
index 3873962cca..87167b572c 100644
--- a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
+++ b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
@@ -45,12 +45,14 @@ Foam::sampledPlane::sampledPlane
     const word& name,
     const polyMesh& mesh,
     const plane& planeDesc,
-    const keyType& zoneKey
+    const keyType& zoneKey,
+    const bool triangulate
 )
 :
     sampledSurface(name, mesh),
     cuttingPlane(planeDesc),
     zoneKey_(zoneKey),
+    triangulate_(triangulate),
     needsUpdate_(true)
 {
     if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) < 0)
@@ -71,6 +73,7 @@ Foam::sampledPlane::sampledPlane
     sampledSurface(name, mesh, dict),
     cuttingPlane(plane(dict.lookup("basePoint"), dict.lookup("normalVector"))),
     zoneKey_(keyType::null),
+    triangulate_(dict.lookupOrDefault("triangulate", true)),
     needsUpdate_(true)
 {
     // make plane relative to the coordinateSystem (Cartesian)
@@ -138,11 +141,11 @@ bool Foam::sampledPlane::update()
 
     if (selectedCells.empty())
     {
-        reCut(mesh(), true);    // always triangulate. Note:Make option?
+        reCut(mesh(), triangulate_);
     }
     else
     {
-        reCut(mesh(), true, selectedCells);
+        reCut(mesh(), triangulate_, selectedCells);
     }
 
     if (debug)
@@ -250,6 +253,7 @@ void Foam::sampledPlane::print(Ostream& os) const
     os  << "sampledPlane: " << name() << " :"
         << "  base:" << refPoint()
         << "  normal:" << normal()
+        << "  triangulate:" << triangulate_
         << "  faces:" << faces().size()
         << "  points:" << points().size();
 }
diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlane.H b/src/sampling/sampledSurface/sampledPlane/sampledPlane.H
index 1c6a1b8058..3b24a09e52 100644
--- a/src/sampling/sampledSurface/sampledPlane/sampledPlane.H
+++ b/src/sampling/sampledSurface/sampledPlane/sampledPlane.H
@@ -25,7 +25,7 @@ Class
     Foam::sampledPlane
 
 Description
-    A sampledSurface defined by a cuttingPlane. Always triangulated.
+    A sampledSurface defined by a cuttingPlane. Triangulated by default.
 
 Note
     Does not actually cut until update() called.
@@ -60,6 +60,9 @@ class sampledPlane
         //- If restricted to zones, name of this zone or a regular expression
         keyType zoneKey_;
 
+        //- Triangulated faces or keep faces as is
+        const bool triangulate_;
+
         //- Track if the surface needs an update
         mutable bool needsUpdate_;
 
@@ -92,7 +95,8 @@ public:
             const word& name,
             const polyMesh& mesh,
             const plane& planeDesc,
-            const keyType& zoneKey = word::null
+            const keyType& zoneKey = word::null,
+            const bool triangulate = true
         );
 
         //- Construct from dictionary
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
index 67704bb9df..75ea14eb67 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.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
@@ -140,6 +140,12 @@ public:
                     const label i
                 ) const = 0;
 
+                //- Return access to chemical source terms [kg/m3/s]
+                virtual DimensionedField<scalar, volMesh>& RR
+                (
+                    const label i
+                ) = 0;
+
 
             // Chemistry solution
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
index 7b3d1c663a..475fc7c8fd 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
@@ -211,6 +211,12 @@ public:
                 const label i
             ) const;
 
+            //- Return non const access to chemical source terms [kg/m3/s]
+            virtual DimensionedField<scalar, volMesh>& RR
+            (
+                const label i
+            );
+
             //- Solve the reaction system for the given start time and time
             //  step and return the characteristic time
             virtual scalar solve(const scalar t0, const scalar deltaT);
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H
index d128443fe9..c009a11e46 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.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
@@ -78,5 +78,15 @@ Foam::chemistryModel<CompType, ThermoType>::RR
     return RR_[i];
 }
 
+template<class CompType, class ThermoType>
+Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::chemistryModel<CompType, ThermoType>::RR
+(
+    const label i
+)
+{
+    return RR_[i];
+}
+
 
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C b/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C
index dbdc804d8b..f10a1ce48b 100644
--- a/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C
+++ b/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.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
@@ -49,7 +49,7 @@ const ThermoType& Foam::multiComponentMixture<ThermoType>::constructSpeciesData
 template<class ThermoType>
 void Foam::multiComponentMixture<ThermoType>::correctMassFractions()
 {
-    // It changes Yt patches to "calculated"
+    // Multiplication by 1.0 changes Yt patches to "calculated"
     volScalarField Yt("Yt", 1.0*Y_[0]);
 
     for (label n=1; n<Y_.size(); n++)
@@ -57,6 +57,17 @@ void Foam::multiComponentMixture<ThermoType>::correctMassFractions()
         Yt += Y_[n];
     }
 
+    if (mag(max(Yt).value()) < ROOTVSMALL)
+    {
+        FatalErrorIn
+        (
+            "void Foam::multiComponentMixture<ThermoType>::"
+            "correctMassFractions()"
+        )
+            << "Sum of mass fractions is zero for species " << this->species()
+            << exit(FatalError);
+    }
+
     forAll(Y_, n)
     {
         Y_[n] /= Yt;
diff --git a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C
index cfff71b9ed..f94db5946b 100644
--- a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C
@@ -53,4 +53,35 @@ Foam::basicSolidChemistryModel::~basicSolidChemistryModel()
 {}
 
 
+const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::basicSolidChemistryModel::RR(const label i) const
+{
+    notImplemented
+    (
+        "const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&"
+        "basicSolidChemistryModel::RR(const label)"
+    );
+    return (DimensionedField<scalar, volMesh>::null());
+}
+
+
+Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::basicSolidChemistryModel::RR(const label i)
+{
+    notImplemented
+    (
+        "Foam::DimensionedField<Foam::scalar, Foam::volMesh>&"
+        "basicSolidChemistryModel::RR(const label)"
+    );
+
+    return dynamic_cast<DimensionedField<scalar, volMesh>&>
+    (
+        const_cast<DimensionedField<scalar, volMesh>& >
+        (
+            DimensionedField<scalar, volMesh>::null()
+        )
+    );
+}
+
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H
index 2134bffca1..50d6a2e96b 100644
--- a/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H
+++ b/src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H
@@ -149,6 +149,15 @@ public:
 
         //- Calculates the reaction rates
         virtual void calculate() = 0;
+
+        //- Return const access to the total source terms
+        virtual const DimensionedField<scalar, volMesh>& RR
+        (
+            const label i
+        ) const;
+
+        //- Return non-const access to the total source terms
+        virtual DimensionedField<scalar, volMesh>& RR(const label i);
 };
 
 
diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
index e301a1b85e..4962bf716b 100644
--- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
@@ -535,14 +535,21 @@ updateConcsInReactionI
         c[si] = max(0.0, c[si]);
     }
 
+    scalar sr = 0.0;
     forAll(R.rhs(), s)
     {
         label si = R.rhs()[s].index;
         const scalar rhoR = this->solidThermo_[si].rho(p, T);
-        const scalar sr = rhoR/rhoL;
+        sr = rhoR/rhoL;
         c[si] += dt*sr*omeg;
         c[si] = max(0.0, c[si]);
     }
+
+    forAll(R.grhs(), g)
+    {
+        label gi = R.grhs()[g].index;
+        c[gi + this->nSolids_] += (1.0 - sr)*omeg*dt;
+    }
 }
 
 
@@ -561,24 +568,11 @@ updateRRInReactionI
     simpleMatrix<scalar>& RR
 ) const
 {
-    const Reaction<SolidThermo>& R = this->reactions_[index];
-    scalar rhoL = 0.0;
-    forAll(R.lhs(), s)
-    {
-        label si = R.lhs()[s].index;
-        rhoL = this->solidThermo_[si].rho(p, T);
-        RR[si][rRef] -= pr*corr;
-        RR[si][lRef] += pf*corr;
-    }
-
-    forAll(R.rhs(), s)
-    {
-        label si = R.rhs()[s].index;
-        const scalar rhoR = this->solidThermo_[si].rho(p, T);
-        const scalar sr = rhoR/rhoL;
-        RR[si][lRef] -= sr*pf*corr;
-        RR[si][rRef] += sr*pr*corr;
-    }
+    notImplemented
+    (
+        "void Foam::pyrolysisChemistryModel<CompType, SolidThermo,GasThermo>::"
+        "updateRRInReactionI"
+    );
 }
 
 
@@ -666,7 +660,6 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
 
                 scalar newCp = 0.0;
                 scalar newhi = 0.0;
-                scalar invRho = 0.0;
                 scalarList dcdt = (c - c0)/dt;
 
                 for (label i=0; i<this->nSolids_; i++)
@@ -675,7 +668,6 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
                     scalar Yi = c[i]/cTot;
                     newCp += Yi*this->solidThermo_[i].Cp(pi, Ti);
                     newhi -= dYi*this->solidThermo_[i].Hc();
-                    invRho += Yi/this->solidThermo_[i].rho(pi, Ti);
                 }
 
                 scalar dTi = (newhi/newCp)*dt;
diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H
index 2d4a11cc90..7ae76c5233 100644
--- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H
+++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H
@@ -137,8 +137,9 @@ public:
             const bool updateC0 = false
         ) const;
 
-        //- Return the reaction rate for reaction r and the reference
-        //  species and charateristic times
+        //- Return the reaction rate for reaction r
+        // NOTE: Currently does not calculate reference specie and
+        // characteristic times (pf, cf,l Ref, etc.)
         virtual scalar omega
         (
             const Reaction<SolidThermo>& r,
@@ -153,8 +154,10 @@ public:
             label& rRef
         ) const;
 
-        //- Return the reaction rate for iReaction and the reference
-        //  species and charateristic times
+
+        //- Return the reaction rate for iReaction
+        // NOTE: Currently does not calculate reference specie and
+        // characteristic times (pf, cf,l Ref, etc.)
         virtual scalar omegaI
         (
             label iReaction,
@@ -169,6 +172,7 @@ public:
             label& rRef
         ) const;
 
+
         //- Calculates the reaction rates
         virtual void calculate();
 
@@ -186,6 +190,7 @@ public:
 
 
         //- Update matrix RR for reaction i. Used by EulerImplicit
+        // (Not implemented)
         virtual void updateRRInReactionI
         (
             const label i,
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
index 03076d79ce..bf67f1ca4d 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,6 +65,9 @@ class solidChemistryModel
 {
     // Private Member Functions
 
+        //- Disallow copy constructor
+        solidChemistryModel(const solidChemistryModel&);
+
         //- Disallow default bitwise assignment
         void operator=(const solidChemistryModel&);
 
@@ -151,6 +154,7 @@ public:
             label& rRef
         ) const = 0;
 
+
         //- Return the reaction rate for iReaction and the reference
         //  species and charateristic times
         virtual scalar omegaI
@@ -167,6 +171,10 @@ public:
             label& rRef
         ) const = 0;
 
+        //- Calculates the reaction rates
+        virtual void calculate() = 0;
+
+
         //- Update concentrations in reaction i given dt and reaction rate
         // omega used by sequential solver
         virtual void updateConcsInReactionI
@@ -194,11 +202,8 @@ public:
             simpleMatrix<scalar>& RR
         ) const = 0;
 
-        //- Calculates the reaction rates
-        virtual void calculate() = 0;
 
-
-        // Chemistry model functions
+        // Solid Chemistry model functions
 
             //- Return const access to the chemical source terms for solids
             inline const DimensionedField<scalar, volMesh>& RRs
@@ -209,13 +214,6 @@ public:
             //- Return total solid source term
             inline tmp<DimensionedField<scalar, volMesh> > RRs() const;
 
-            //- Return const access to the total source terms
-            inline const DimensionedField<scalar, volMesh>& RR
-            (
-                const label i
-            ) const;
-
-
             //- Solve the reaction system for the given start time and time
             //  step and return the characteristic time
             virtual scalar solve(const scalar t0, const scalar deltaT) = 0;
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
index 92b461e11e..f2d0dad79a 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,16 +95,4 @@ Foam::solidChemistryModel<CompType, SolidThermo>::RRs() const
 }
 
 
-template<class CompType, class SolidThermo>
-inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::solidChemistryModel<CompType, SolidThermo>::RR
-(
-    const label i
-) const
-{
-    notImplemented("solidChemistryModel::RR(const label)");
-    return (DimensionedField<scalar, volMesh>::null());
-}
-
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H
index 12f842a054..447020e2f7 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H
@@ -52,9 +52,8 @@ namespace Foam
     defineTemplateTypeNameAndDebugWithName                                    \
     (                                                                         \
         SS##Schem##Comp##SThermo##GThermo,                                    \
-        (#SS"<" + word(Schem::typeName_())                                    \
-      + "<"#Comp"," + SThermo::typeName()                                     \
-      + "," + GThermo::typeName() + ">>").c_str(),                            \
+        (#SS"<"#Schem"<"#Comp"," + SThermo::typeName() + ","                  \
+      + GThermo::typeName() + ">>").c_str(),                                  \
         0                                                                     \
     );                                                                        \
                                                                               \
@@ -77,14 +76,6 @@ namespace Foam
         GThermo                                                               \
     );                                                                        \
                                                                               \
-    makeSolidChemistrySolverType                                              \
-    (                                                                         \
-        EulerImplicit,                                                        \
-        SolidChem,                                                            \
-        Comp,                                                                 \
-        SThermo,                                                              \
-        GThermo                                                               \
-    );                                                                        \
                                                                               \
     makeSolidChemistrySolverType                                              \
     (                                                                         \
@@ -104,7 +95,6 @@ namespace Foam
         GThermo                                                               \
     );
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H b/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
index 5193a4b9de..07b6ffb050 100644
--- a/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
+++ b/src/thermophysicalModels/solidThermo/solidThermo/makeSolidThermo.H
@@ -30,7 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #ifndef makeSolidThermo_H
-#define makesolidThermo_H
+#define makeSolidThermo_H
 
 #include "addToRunTimeSelectionTable.H"
 
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties
index 5408eb7928..9228607e91 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         yes;
     transient       yes;
     cellValueSourceCorrection yes;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties
index 5408eb7928..9228607e91 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         yes;
     transient       yes;
     cellValueSourceCorrection yes;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties
index 5408eb7928..9228607e91 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         yes;
     transient       yes;
     cellValueSourceCorrection yes;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
index 7c71c8971e..e164395342 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
@@ -21,6 +21,7 @@ solution
     transient       yes;
     coupled         true;
     cellValueSourceCorrection on;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
index 4b5b5393e4..36b287137d 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         true;
     transient       yes;
     cellValueSourceCorrection on;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
index 06571ffdd1..150adb2953 100644
--- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
+++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
@@ -21,6 +21,7 @@ solution
     coupled         false;
     transient       yes;
     cellValueSourceCorrection off;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
index 874c386653..2fa48b6da0 100644
--- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
+++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
@@ -21,6 +21,7 @@ solution
     coupled         false;
     transient       yes;
     cellValueSourceCorrection off;
+    maxCo           0.3;
 
     interpolationSchemes
     {
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties
index c96b27a4b4..6e01c75332 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         no;
     transient       yes;
     cellValueSourceCorrection no;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties
index e6afd3b55d..4e8221412f 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         yes;
     transient       yes;
     cellValueSourceCorrection yes;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
index 5652c17be9..e38561ee8f 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         no;
     transient       yes;
     cellValueSourceCorrection no;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties
index 94ede2729a..2eba294927 100644
--- a/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         true;
     transient       yes;
     cellValueSourceCorrection on;
+    maxCo           0.3;
 
     sourceTerms
     {
@@ -41,7 +42,7 @@ solution
         thermo:mu       cell;
         T               cell;
         Cp              cell;
-        kapppa          cell;
+        kappa           cell;
         p               cell;
     }
 
diff --git a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
index 41fe1084d4..1ce6679e09 100644
--- a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         false;
     transient       yes;
     cellValueSourceCorrection off;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
index 70524de8a8..c6e386958c 100644
--- a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         true;
     transient       yes;
     cellValueSourceCorrection on;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties b/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties
index b0989cc8b9..4d9cfe2711 100644
--- a/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties
+++ b/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties
@@ -21,6 +21,7 @@ solution
     coupled         true;
     transient       yes;
     cellValueSourceCorrection on;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphawater b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.air
similarity index 97%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphawater
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.air
index 0f3dcd4639..210222182a 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphawater
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.air
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha1;
+    object      alpha.air;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphaoil b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.mercury
similarity index 97%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphaoil
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.mercury
index 0f3dcd4639..cbb4647b47 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphaoil
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.mercury
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha1;
+    object      alpha.mercury;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphamercury b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.oil
similarity index 97%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphamercury
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.oil
index 0f3dcd4639..e8b4d42497 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphamercury
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.oil
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha1;
+    object      alpha.oil;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphaair b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.water
similarity index 100%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphaair
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.water
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphas b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphas
deleted file mode 100644
index 0f3dcd4639..0000000000
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphas
+++ /dev/null
@@ -1,44 +0,0 @@
-/*--------------------------------*- 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;
-    object      alpha1;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    rotor
-    {
-        type            zeroGradient;
-    }
-
-    stator
-    {
-        type            zeroGradient;
-    }
-
-    front
-    {
-        type            empty;
-    }
-
-    back
-    {
-        type            empty;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaair b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.air
similarity index 99%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaair
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.air
index b25f228620..a53034d33d 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaair
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.air
@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      alphaair;
+    object      alpha.air;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphamercury b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.mercury
similarity index 99%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphamercury
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.mercury
index 0e6d124fd0..76002901ed 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphamercury
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.mercury
@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      alphamercury;
+    object      alpha.mercury;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaoil b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.oil
similarity index 99%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaoil
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.oil
index 2311855110..56d642d0d9 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaoil
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.oil
@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      alphaoil;
+    object      alpha.oil;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphawater b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.water
similarity index 99%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphawater
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.water
index d61eb80ddf..766d8c1462 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphawater
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.water
@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      alphawater;
+    object      alpha.water;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/polyMesh/boundary b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/polyMesh/boundary
index 292f25b806..188a0f0c58 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/polyMesh/boundary
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/polyMesh/boundary
@@ -32,12 +32,14 @@ FoamFile
     front
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          3072;
         startFace       6336;
     }
     back
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          3072;
         startFace       9408;
     }
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/transportProperties b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/transportProperties
index e338328156..3cfdbabfaf 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/transportProperties
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/transportProperties
@@ -46,8 +46,6 @@ phases
      }
 );
 
-refPhase        air;
-
 sigmas
 (
     (air water) 0.07
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/setFieldsDict b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/setFieldsDict
index d5dd8aca43..82d9ffc783 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/setFieldsDict
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/setFieldsDict
@@ -17,10 +17,10 @@ FoamFile
 
 defaultFieldValues
 (
-    volScalarFieldValue alphaair 1
-    volScalarFieldValue alphawater 0
-    volScalarFieldValue alphaoil 0
-    volScalarFieldValue alphamercury 0
+    volScalarFieldValue alpha.air 1
+    volScalarFieldValue alpha.water 0
+    volScalarFieldValue alpha.oil 0
+    volScalarFieldValue alpha.mercury 0
 );
 
 regions
@@ -30,10 +30,10 @@ regions
         box (0 0 -1) (1 1 1);
         fieldValues
         (
-            volScalarFieldValue alphawater 1
-            volScalarFieldValue alphaoil 0
-            volScalarFieldValue alphamercury 0
-            volScalarFieldValue alphaair 0
+            volScalarFieldValue alpha.water 1
+            volScalarFieldValue alpha.oil 0
+            volScalarFieldValue alpha.mercury 0
+            volScalarFieldValue alpha.air 0
         );
     }
     boxToCell
@@ -41,10 +41,10 @@ regions
         box (0 -1 -1) (1 0 1);
         fieldValues
         (
-            volScalarFieldValue alphawater 0
-            volScalarFieldValue alphaoil 1
-            volScalarFieldValue alphamercury 0
-            volScalarFieldValue alphaair 0
+            volScalarFieldValue alpha.water 0
+            volScalarFieldValue alpha.oil 1
+            volScalarFieldValue alpha.mercury 0
+            volScalarFieldValue alpha.air 0
         );
     }
     boxToCell
@@ -52,10 +52,10 @@ regions
         box (-1 -1 -1) (0 0 1);
         fieldValues
         (
-            volScalarFieldValue alphawater 0
-            volScalarFieldValue alphaoil 0
-            volScalarFieldValue alphamercury 1
-            volScalarFieldValue alphaair 0
+            volScalarFieldValue alpha.water 0
+            volScalarFieldValue alpha.oil 0
+            volScalarFieldValue alpha.mercury 1
+            volScalarFieldValue alpha.air 0
         );
     }
 );