diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
index 16ef1e517b..46c0171e5a 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
@@ -1210,12 +1210,13 @@ int main(int argc, char *argv[])
         << endl;
 
     label nSide = 0;
+
     forAll(zoneSidePatch, zoneI)
     {
         if (oneD)
         {
-            // Always add empty patches, one per zone.
-            word patchName = faceZones[zoneI].name() + "_" + "side";
+            // Reuse single empty patch.
+            word patchName = "oneDEmptPatch";
 
             zoneSidePatch[zoneI] = addPatch<emptyPolyPatch>
             (
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index 90e6a379e0..996e6fcb37 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -6,6 +6,7 @@
 #include "faceSet.H"
 #include "pointSet.H"
 #include "IOmanip.H"
+#include "emptyPolyPatch.H"
 
 Foam::label Foam::checkTopology
 (
@@ -21,6 +22,29 @@ Foam::label Foam::checkTopology
     // Check if the boundary definition is unique
     mesh.boundaryMesh().checkDefinition(true);
 
+    // Check that empty patches cover all sides of the mesh
+    {
+        label nEmpty = 0;
+        forAll(mesh.boundaryMesh(), patchI)
+        {
+            if (isA<emptyPolyPatch>(mesh.boundaryMesh()[patchI]))
+            {
+                nEmpty += mesh.boundaryMesh()[patchI].size();
+            }
+        }
+        reduce(nEmpty, sumOp<label>());
+        label nTotCells = returnReduce(mesh.cells().size(), sumOp<label>());
+
+        // These are actually warnings, not errors.
+        if (nEmpty % nTotCells)
+        {
+            Info<< " ***Total number of faces on empty patches"
+                << " is not divisible by the number of cells in the mesh."
+                << " Hence this mesh is not 1D or 2D."
+                << endl;
+        }
+    }
+
     // Check if the boundary processor patches are correct
     mesh.boundaryMesh().checkParallelSync(true);
 
@@ -41,6 +65,8 @@ Foam::label Foam::checkTopology
         noFailedChecks++;
     }
 
+
+
     {
         pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
         if (mesh.checkPoints(true, &points))
@@ -74,6 +100,22 @@ Foam::label Foam::checkTopology
         }
     }
 
+    {
+        faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
+        if (mesh.checkFaceVertices(true, &faces))
+        {
+            noFailedChecks++;
+
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nFaces
+                << " faces with out-of-range or duplicate vertices to set "
+                << faces.name() << endl;
+            faces.instance() = mesh.pointsInstance();
+            faces.write();
+        }
+    }
+
     if (allTopology)
     {
         cellSet cells(mesh, "zipUpCells", mesh.nCells()/100);
@@ -91,22 +133,6 @@ Foam::label Foam::checkTopology
         }
     }
 
-    {
-        faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
-        if (mesh.checkFaceVertices(true, &faces))
-        {
-            noFailedChecks++;
-
-            label nFaces = returnReduce(faces.size(), sumOp<label>());
-
-            Info<< "  <<Writing " << nFaces
-                << " faces with out-of-range or duplicate vertices to set "
-                << faces.name() << endl;
-            faces.instance() = mesh.pointsInstance();
-            faces.write();
-        }
-    }
-
     if (allTopology)
     {
         faceSet faces(mesh, "edgeFaces", mesh.nFaces()/100);
diff --git a/src/OpenFOAM/db/Time/Time.C b/src/OpenFOAM/db/Time/Time.C
index d276de95ce..8e509a738b 100644
--- a/src/OpenFOAM/db/Time/Time.C
+++ b/src/OpenFOAM/db/Time/Time.C
@@ -228,6 +228,8 @@ Foam::Time::Time
 
     objectRegistry(*this),
 
+    libs_(),
+
     controlDict_
     (
         IOobject
@@ -257,9 +259,10 @@ Foam::Time::Time
     graphFormat_("raw"),
     runTimeModifiable_(true),
 
-    libs_(controlDict_, "libs"),
     functionObjects_(*this)
 {
+    libs_.open(controlDict_, "libs");
+
     // Explicitly set read flags on objectRegistry so anything constructed
     // from it reads as well (e.g. fvSolution).
     readOpt() = IOobject::MUST_READ_IF_MODIFIED;
@@ -313,6 +316,8 @@ Foam::Time::Time
 
     objectRegistry(*this),
 
+    libs_(),
+
     controlDict_
     (
         IOobject
@@ -343,9 +348,11 @@ Foam::Time::Time
     graphFormat_("raw"),
     runTimeModifiable_(true),
 
-    libs_(controlDict_, "libs"),
     functionObjects_(*this)
 {
+    libs_.open(controlDict_, "libs");
+
+
     // Explicitly set read flags on objectRegistry so anything constructed
     // from it reads as well (e.g. fvSolution).
     readOpt() = IOobject::MUST_READ_IF_MODIFIED;
@@ -401,6 +408,8 @@ Foam::Time::Time
 
     objectRegistry(*this),
 
+    libs_(),
+
     controlDict_
     (
         IOobject
@@ -430,9 +439,10 @@ Foam::Time::Time
     graphFormat_("raw"),
     runTimeModifiable_(true),
 
-    libs_(controlDict_, "libs"),
     functionObjects_(*this)
-{}
+{
+    libs_.open(controlDict_, "libs");
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/db/Time/Time.H b/src/OpenFOAM/db/Time/Time.H
index 28ff775a68..b9da549895 100644
--- a/src/OpenFOAM/db/Time/Time.H
+++ b/src/OpenFOAM/db/Time/Time.H
@@ -75,6 +75,10 @@ class Time
         //- file-change monitor for all registered files
         mutable autoPtr<fileMonitor> monitorPtr_;
 
+        //- Any loaded dynamic libraries. Make sure to construct before
+        //  reading controlDict.
+        dlLibraryTable libs_;
+
         //- The controlDict
         IOdictionary controlDict_;
 
@@ -166,9 +170,6 @@ private:
         //- Is runtime modification of dictionaries allowed?
         Switch runTimeModifiable_;
 
-        //- Any loaded dynamic libraries
-        dlLibraryTable libs_;
-
         //- Function objects executed at start and on ++, +=
         mutable functionObjectList functionObjects_;
 
diff --git a/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C b/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C
index e0a1bbbdd4..8eefcb86d9 100644
--- a/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C
+++ b/src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C
@@ -131,11 +131,7 @@ bool Foam::functionEntries::codeStream::execute
 
     // see if library is loaded
     void* lib = NULL;
-    if
-    (
-        isA<IOdictionary>(topDict(parentDict))
-     && parentDict.dictName() != Time::controlDictName
-    )
+    if (isA<IOdictionary>(topDict(parentDict)))
     {
         lib = libs(parentDict).findLibrary(libPath);
     }
@@ -150,11 +146,7 @@ bool Foam::functionEntries::codeStream::execute
     // avoid compilation if possible by loading an existing library
     if (!lib)
     {
-        if
-        (
-            isA<IOdictionary>(topDict(parentDict))
-         && parentDict.dictName() != Time::controlDictName
-        )
+        if (isA<IOdictionary>(topDict(parentDict)))
         {
             // Cached access to dl libs. Guarantees clean up upon destruction
             // of Time.
@@ -223,11 +215,7 @@ bool Foam::functionEntries::codeStream::execute
         // all processes must wait for compile to finish
         reduce(create, orOp<bool>());
 
-        if
-        (
-            isA<IOdictionary>(topDict(parentDict))
-         && parentDict.dictName() != Time::controlDictName
-        )
+        if (isA<IOdictionary>(topDict(parentDict)))
         {
             // Cached access to dl libs. Guarantees clean up upon destruction
             // of Time.
diff --git a/src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C b/src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C
index deca23890d..93d934a1b5 100644
--- a/src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C
+++ b/src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C
@@ -25,6 +25,12 @@ License
 
 #include "dlLibraryTable.H"
 #include "OSspecific.H"
+#include "long.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::dlLibraryTable, 0);
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -55,7 +61,11 @@ Foam::dlLibraryTable::~dlLibraryTable()
         // bug in dlclose - does not call static destructors of
         // loaded library when actually unloading the library.
         // See https://bugzilla.novell.com/show_bug.cgi?id=680125 and 657627.
-        // Seems related to using a non-system compiler!
+        if (debug)
+        {
+            Info<< "dlLibraryTable::~dlLibraryTable() : closing " << iter()
+                << " with handle " << long(iter.key()) << endl;
+        }
         dlClose(iter.key());
     }
 }
@@ -73,6 +83,12 @@ bool Foam::dlLibraryTable::open
     {
         void* functionLibPtr = dlOpen(functionLibName);
 
+        if (debug)
+        {
+            Info<< "dlLibraryTable::open : opened " << functionLibName
+                << " resulting in handle " << long(functionLibPtr) << endl;
+        }
+
         if (!functionLibPtr)
         {
             if (verbose)
@@ -107,6 +123,12 @@ bool Foam::dlLibraryTable::close
     void* libPtr = findLibrary(functionLibName);
     if (libPtr)
     {
+        if (debug)
+        {
+            Info<< "dlLibraryTable::close : closing " << functionLibName
+                << " with handle " << long(libPtr) << endl;
+        }
+
         erase(libPtr);
 
         if (!dlClose(libPtr))
diff --git a/src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.H b/src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.H
index c3388907d1..f2925fe511 100644
--- a/src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.H
+++ b/src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.H
@@ -63,6 +63,9 @@ class dlLibraryTable
 
 public:
 
+    // Declare name of the class and its debug switch
+    ClassName("dlLibraryTable");
+
     // Constructors
 
         //- Construct null
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
index f8f7ab41b0..92b64e85a7 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
@@ -81,30 +81,6 @@ void Foam::addPatchCellLayer::addVertex
             f[fp++] = pointI;
         }
     }
-
-    // Check for duplicates.
-    if (debug)
-    {
-        label n = 0;
-        for (label i = 0; i < fp; i++)
-        {
-            if (f[i] == pointI)
-            {
-                n++;
-
-                if (n == 2)
-                {
-                    f.setSize(fp);
-                    FatalErrorIn
-                    (
-                        "addPatchCellLayer::addVertex(const label, face&"
-                        ", label&)"
-                    )   << "Point " << pointI << " already present in face "
-                        << f << abort(FatalError);
-                }
-            }
-        }
-    }
 }
 
 
@@ -1641,11 +1617,11 @@ void Foam::addPatchCellLayer::setRefinement
                         {
                             label offset =
                                 addedPoints_[vStart].size() - numEdgeSideFaces;
-                            for (label ioff = offset; ioff > 0; ioff--)
+                            for (label ioff = offset-1; ioff >= 0; ioff--)
                             {
                                 addVertex
                                 (
-                                    addedPoints_[vStart][ioff-1],
+                                    addedPoints_[vStart][ioff],
                                     newFace,
                                     newFp
                                 );
@@ -1660,6 +1636,51 @@ void Foam::addPatchCellLayer::setRefinement
 
                         newFace.setSize(newFp);
 
+                        if (debug)
+                        {
+                            labelHashSet verts(2*newFace.size());
+                            forAll(newFace, fp)
+                            {
+                                if (!verts.insert(newFace[fp]))
+                                {
+                                    FatalErrorIn
+                                    (
+                                        "addPatchCellLayer::setRefinement(..)"
+                                    )   << "Duplicate vertex in face"
+                                        << " to be added." << nl
+                                        << "newFace:" << newFace << nl
+                                        << "points:"
+                                        <<  UIndirectList<point>
+                                            (
+                                                meshMod.points(),
+                                                newFace
+                                            ) << nl
+                                        << "Layer:" << i
+                                        << " out of:" << numEdgeSideFaces << nl
+                                        << "ExtrudeEdge:" << meshEdgeI
+                                        << " at:"
+                                        <<  mesh_.edges()[meshEdgeI].line
+                                            (
+                                                mesh_.points()
+                                            ) << nl
+                                        << "string:" << stringedVerts
+                                        << "stringpoints:"
+                                        << UIndirectList<point>
+                                            (
+                                                pp.localPoints(),
+                                                stringedVerts
+                                            ) << nl
+                                        << "stringNLayers:"
+                                        <<  UIndirectList<label>
+                                            (
+                                                nPointLayers,
+                                                stringedVerts
+                                            ) << nl
+                                        << abort(FatalError);
+                                }
+                            }
+                        }
+
                         label nbrFaceI = nbrFace
                         (
                             pp.edgeFaces(),
diff --git a/src/finiteVolume/cfdTools/general/pimpleLoop/pimpleLoop.H b/src/finiteVolume/cfdTools/general/pimpleLoop/pimpleLoop.H
deleted file mode 100644
index 94e179f16f..0000000000
--- a/src/finiteVolume/cfdTools/general/pimpleLoop/pimpleLoop.H
+++ /dev/null
@@ -1,128 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 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::pimpleLoop
-
-Description
-    PIMPLE loop class to formalise the iteration and automate the handling
-    of the "finalIteration" mesh data entry.
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef pimpleLoop_H
-#define pimpleLoop_H
-
-#include "fvMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                         Class pimpleLoop Declaration
-\*---------------------------------------------------------------------------*/
-
-class pimpleLoop
-{
-    // Private data
-
-        //- Reference to the mesh
-        fvMesh& mesh_;
-
-        //- Number of PIMPLE correctors
-        const int nCorr_;
-
-        //- Current PIMPLE corrector
-        int corr_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        pimpleLoop(const pimpleLoop&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const pimpleLoop&);
-
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        pimpleLoop(fvMesh& mesh, const int nCorr)
-        :
-            mesh_(mesh),
-            nCorr_(nCorr),
-            corr_(0)
-        {}
-
-
-    //- Destructor
-    ~pimpleLoop()
-    {}
-
-
-    // Member Functions
-
-        bool loop()
-        {
-            if (finalIter())
-            {
-                mesh_.data::add("finalIteration", true);
-            }
-
-            return corr_ < nCorr_;
-        }
-
-        bool finalIter() const
-        {
-            return corr_ == nCorr_-1;
-        }
-
-
-    // Member Operators
-
-        void operator++(int)
-        {
-            if (finalIter())
-            {
-                mesh_.data::remove("finalIteration");
-            }
-
-            corr_++;
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
index cd5d46c824..3402992f46 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
@@ -137,7 +137,7 @@ Foam::pimpleControl::pimpleControl(fvMesh& mesh)
             {
                 Info<< "    field " << residualControl_[i].name << token::TAB
                     << ": relTol " << residualControl_[i].relTol
-                    << ", absTol " << residualControl_[i].absTol
+                    << ", tolerance " << residualControl_[i].absTol
                     << nl;
             }
             Info<< endl;
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C
index b0cb59e927..e6104de03d 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C
@@ -68,7 +68,7 @@ bool Foam::simpleControl::criteriaSatisfied()
             {
                 Info<< algorithmName_ << " solution statistics:" << endl;
 
-                Info<< "    " << variableName << ": abs tol = " << residual
+                Info<< "    " << variableName << ": tolerance = " << residual
                     << " (" << residualControl_[fieldI].absTol << ")"
                     << endl;
             }
@@ -96,7 +96,7 @@ Foam::simpleControl::simpleControl(fvMesh& mesh)
         forAll(residualControl_, i)
         {
             Info<< "    field " << residualControl_[i].name << token::TAB
-                << " absTol " << residualControl_[i].absTol
+                << " tolerance " << residualControl_[i].absTol
                 << nl;
         }
         Info<< endl;
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
index 55fa22435b..e68c6a2d4d 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
@@ -68,7 +68,7 @@ void Foam::solutionControl::read(const bool absTolOnly)
                 if (iter().isDict())
                 {
                     const dictionary& fieldDict(iter().dict());
-                    fd.absTol = readScalar(fieldDict.lookup("absTol"));
+                    fd.absTol = readScalar(fieldDict.lookup("tolerance"));
                     fd.relTol = readScalar(fieldDict.lookup("relTol"));
                     fd.initialResidual = 0.0;
                 }
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/empty/emptyFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/empty/emptyFvPatchField.C
index 9848e2e9f5..5488a03885 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/empty/emptyFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/empty/emptyFvPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -138,18 +138,21 @@ emptyFvPatchField<Type>::emptyFvPatchField
 template<class Type>
 void emptyFvPatchField<Type>::updateCoeffs()
 {
-    if
-    (
-        this->patch().patch().size()
-      % this->dimensionedInternalField().mesh().nCells()
-    )
-    {
-        FatalErrorIn("emptyFvPatchField<Type>::updateCoeffs()")
-            << "This mesh contains patches of type empty but is not 1D or 2D\n"
-               "    by virtue of the fact that the number of faces of this\n"
-               "    empty patch is not divisible by the number of cells."
-            << exit(FatalError);
-    }
+    //- Check moved to checkMesh. Test here breaks down if multiple empty
+    //  patches.
+    //if
+    //(
+    //    this->patch().patch().size()
+    //  % this->dimensionedInternalField().mesh().nCells()
+    //)
+    //{
+    //    FatalErrorIn("emptyFvPatchField<Type>::updateCoeffs()")
+    //        << "This mesh contains patches of type empty but is not"
+    //        << "1D or 2D\n"
+    //           "    by virtue of the fact that the number of faces of this\n"
+    //           "    empty patch is not divisible by the number of cells."
+    //        << exit(FatalError);
+    //}
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index 754677803c..06e2b1e462 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -518,7 +518,7 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::Re
     const scalar muc
 ) const
 {
-    return rhoc*mag(U - Uc_)*d/muc;
+    return rhoc*mag(U - Uc_)*d/(muc + ROOTVSMALL);
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 0d605da15d..664fac9781 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -63,7 +63,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::CpEff
 
 template<class ParcelType>
 template<class TrackData>
-Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
+Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HsEff
 (
     TrackData& td,
     const scalar p,
@@ -74,9 +74,9 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
 ) const
 {
     return
-        this->Y_[GAS]*td.cloud().composition().H(idG, YGas_, p, T)
-      + this->Y_[LIQ]*td.cloud().composition().H(idL, YLiquid_, p, T)
-      + this->Y_[SLD]*td.cloud().composition().H(idS, YSolid_, p, T);
+        this->Y_[GAS]*td.cloud().composition().Hs(idG, YGas_, p, T)
+      + this->Y_[LIQ]*td.cloud().composition().Hs(idL, YLiquid_, p, T)
+      + this->Y_[SLD]*td.cloud().composition().Hs(idS, YSolid_, p, T);
 }
 
 
@@ -326,7 +326,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
 
 
-
     // Heat transfer
     // ~~~~~~~~~~~~~
 
@@ -368,7 +367,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             d0,
             U0,
             rho0,
-            0.5*(mass0 + mass1),
+            mass0,
             Su,
             dUTrans,
             Spu
@@ -383,25 +382,41 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         // Transfer mass lost from particle to carrier mass source
         forAll(YGas_, i)
         {
+            scalar dm = np0*dMassGas[i];
             label gid = composition.localToGlobalCarrierId(GAS, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
+            scalar hs = composition.carrier().Hs(gid, T0);
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().UTrans()[cellI] += dm*U0;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
         forAll(YLiquid_, i)
         {
+            scalar dm = np0*dMassLiquid[i];
             label gid = composition.localToGlobalCarrierId(LIQ, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
+            scalar hs = composition.carrier().Hs(gid, T0);
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().UTrans()[cellI] += dm*U0;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 /*
         // No mapping between solid components and carrier phase
         forAll(YSolid_, i)
         {
+            scalar dm = np0*dMassSolid[i];
             label gid = composition.localToGlobalCarrierId(SLD, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
+            scalar hs = composition.carrier().Hs(gid, T0);
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().UTrans()[cellI] += dm*U0;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 */
         forAll(dMassSRCarrier, i)
         {
-            td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
+            scalar dm = np0*dMassSRCarrier[i];
+            scalar hs = composition.carrier().Hs(i, T0);
+            td.cloud().rhoTrans(i)[cellI] += dm;
+            td.cloud().UTrans()[cellI] += dm*U0;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 
         // Update momentum transfer
@@ -421,36 +436,38 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // Remove the particle when mass falls below minimum threshold
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (mass1 < td.cloud().constProps().minParticleMass())
+    if (np0*mass1 < td.cloud().constProps().minParticleMass())
     {
         td.keepParticle = false;
 
         if (td.cloud().solution().coupled())
         {
+            scalar dm = np0*mass1;
+
             // Absorb parcel into carrier phase
             forAll(YGas_, i)
             {
                 label gid = composition.localToGlobalCarrierId(GAS, i);
-                td.cloud().rhoTrans(gid)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[GAS]*YGas_[i];
             }
             forAll(YLiquid_, i)
             {
                 label gid = composition.localToGlobalCarrierId(LIQ, i);
-                td.cloud().rhoTrans(gid)[cellI] +=
-                    np0*mass1*YMix[LIQ]*YLiquid_[i];
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[LIQ]*YLiquid_[i];
             }
 /*
             // No mapping between solid components and carrier phase
             forAll(YSolid_, i)
             {
                 label gid = composition.localToGlobalCarrierId(SLD, i);
-                td.cloud().rhoTrans(gid)[cellI] +=
-                    np0*mass1*YMix[SLD]*YSolid_[i];
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[SLD]*YSolid_[i];
             }
 */
-            td.cloud().UTrans()[cellI] += np0*mass1*U1;
+            td.cloud().UTrans()[cellI] += dm*U1;
 
-            // enthalpy transfer accounted for via change in mass fractions
+            td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T1, idG, idL, idS);
+
+            td.cloud().addToMassPhaseChange(dm);
         }
     }
 
@@ -531,28 +548,35 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 
     Sh -= dMassTot*td.cloud().constProps().LDevol()/dt;
 
-    // Molar average molecular weight of carrier mix
-    const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
-
     // Update molar emissions
-    forAll(dMassDV, i)
+    if (td.cloud().heatTransfer().BirdCorrection())
     {
+        // Molar average molecular weight of carrier mix
+        const scalar Wc =
+            max(SMALL, this->rhoc_*specie::RR*this->Tc_/this->pc_);
+
         // Note: hardcoded gaseous diffusivities for now
         // TODO: add to carrier thermo
         const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
-        const label id = composition.localToGlobalCarrierId(GAS, i);
-        const scalar Cp = composition.carrier().Cp(id, Ts);
-        const scalar W = composition.carrier().W(id);
-        const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 
-        // Dab calc'd using API vapour mass diffusivity function
-        const scalar Dab =
-            3.6059e-3*(pow(1.8*Ts, 1.75))*sqrt(1.0/W + 1.0/Wc)/(this->pc_*beta);
+        forAll(dMassDV, i)
+        {
+            const label id = composition.localToGlobalCarrierId(GAS, i);
+            const scalar Cp = composition.carrier().Cp(id, Ts);
+            const scalar W = composition.carrier().W(id);
+            const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 
-        N += Ni;
-        NCpW += Ni*Cp*W;
-        Cs[id] += Ni*d/(2.0*Dab);
-     }
+            // Dab calc'd using API vapour mass diffusivity function
+            const scalar Dab =
+                3.6059e-3*(pow(1.8*Ts, 1.75))
+               *sqrt(1.0/W + 1.0/Wc)
+               /(this->pc_*beta);
+
+            N += Ni;
+            NCpW += Ni*Cp*W;
+            Cs[id] += Ni*d/(2.0*Dab);
+        }
+    }
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index cbb6bb8ead..7e9e16bbc7 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -134,9 +134,9 @@ private:
             const label idS
         ) const;
 
-        //- Return the mixture effective enthalpy
+        //- Return the mixture effective sensible enthalpy
         template<class TrackData>
-        scalar HEff
+        scalar HsEff
         (
             TrackData& td,
             const scalar p,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index b5a2706732..7133014464 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -213,9 +213,17 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
         sumYiCbrtW += Ys[i]*cbrtW;
     }
 
+    Cps = max(Cps, ROOTVSMALL);
+
     rhos *= pc_/(specie::RR*T);
+    rhos = max(rhos, ROOTVSMALL);
+
     mus /= sumYiSqrtW;
+    mus = max(mus, ROOTVSMALL);
+
     kappas /= sumYiCbrtW;
+    kappas = max(kappas, ROOTVSMALL);
+
     Prs = Cps*mus/kappas;
 }
 
@@ -335,7 +343,9 @@ void Foam::ReactingParcel<ParcelType>::calc
     Res = this->Re(U0, d0, rhos, mus);
 
     // Update particle component mass and mass fractions
-    scalar mass1 = updateMassFraction(mass0, dMassPC, Y_);
+    scalarField dMass(dMassPC);
+
+    scalar mass1 = updateMassFraction(mass0, dMass, Y_);
 
 
     // Heat transfer
@@ -379,7 +389,7 @@ void Foam::ReactingParcel<ParcelType>::calc
             d0,
             U0,
             rho0,
-            0.5*(mass0 + mass1),
+            mass0,
             Su,
             dUTrans,
             Spu
@@ -390,11 +400,16 @@ void Foam::ReactingParcel<ParcelType>::calc
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if (td.cloud().solution().coupled())
     {
-        // Transfer mass lost from particle to carrier mass source
-        forAll(dMassPC, i)
+        // Transfer mass lost to carrier mass and enthalpy sources
+        forAll(dMass, i)
         {
+            scalar dm = np0*dMass[i];
             label gid = composition.localToGlobalCarrierId(0, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
+            scalar hs = composition.carrier().Hs(gid, T0);
+
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().UTrans()[cellI] += dm*U0;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 
         // Update momentum transfer
@@ -413,21 +428,27 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     // Remove the particle when mass falls below minimum threshold
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (mass1 < td.cloud().constProps().minParticleMass())
+    if (np0*mass1 < td.cloud().constProps().minParticleMass())
     {
         td.keepParticle = false;
 
         if (td.cloud().solution().coupled())
         {
+            scalar dm = np0*mass1;
+
             // Absorb parcel into carrier phase
             forAll(Y_, i)
             {
+                scalar dmi = dm*Y_[i];
                 label gid = composition.localToGlobalCarrierId(0, i);
-                td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
-            }
-            td.cloud().UTrans()[cellI] += np0*mass1*U1;
+                scalar hs = composition.carrier().Hs(gid, T1);
 
-            // enthalpy transfer accounted for via change in mass fractions
+                td.cloud().rhoTrans(gid)[cellI] += dmi;
+                td.cloud().hsTrans()[cellI] += dmi*hs;
+            }
+            td.cloud().UTrans()[cellI] += dm*U1;
+
+            td.cloud().addToMassPhaseChange(dm);
         }
     }
 
@@ -514,33 +535,44 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     // Add to cumulative phase change mass
     td.cloud().addToMassPhaseChange(this->nParticle_*dMassTot);
 
-    // Average molecular weight of carrier mix - assumes perfect gas
-    const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
-
-    forAll(YComponents, i)
+    forAll(dMassPC, i)
     {
         const label idc = composition.localToGlobalCarrierId(idPhase, i);
         const label idl = composition.globalIds(idPhase)[i];
 
         const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, T);
         Sh -= dMassPC[i]*dh/dt;
+    }
 
-        // Update particle surface thermo properties
-        const scalar Dab =
-            composition.liquids().properties()[idl].D(pc_, Ts, Wc);
 
-        const scalar Cp = composition.carrier().Cp(idc, Ts);
-        const scalar W = composition.carrier().W(idc);
-        const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
+    // Update molar emissions
+    if (td.cloud().heatTransfer().BirdCorrection())
+    {
+        // Average molecular weight of carrier mix - assumes perfect gas
+        const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
 
-        // Molar flux of species coming from the particle (kmol/m^2/s)
-        N += Ni;
 
-        // Sum of Ni*Cpi*Wi of emission species
-        NCpW += Ni*Cp*W;
+        forAll(dMassPC, i)
+        {
+            const label idc = composition.localToGlobalCarrierId(idPhase, i);
+            const label idl = composition.globalIds(idPhase)[i];
 
-        // Concentrations of emission species
-        Cs[idc] += Ni*d/(2.0*Dab);
+            const scalar Cp = composition.carrier().Cp(idc, Ts);
+            const scalar W = composition.carrier().W(idc);
+            const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
+
+            const scalar Dab =
+                composition.liquids().properties()[idl].D(pc_, Ts, Wc);
+
+            // Molar flux of species coming from the particle (kmol/m^2/s)
+            N += Ni;
+
+            // Sum of Ni*Cpi*Wi of emission species
+            NCpW += Ni*Cp*W;
+
+            // Concentrations of emission species
+            Cs[idc] += Ni*d/(2.0*Dab);
+        }
     }
 }
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 6397e29437..dbc63e7c2e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -153,7 +153,10 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
     mus = td.muInterp().interpolate(this->position(), tetIs)/TRatio;
 
     Pr = td.cloud().constProps().Pr();
+    Pr = max(ROOTVSMALL, Pr);
+
     kappas = Cpc_*mus/Pr;
+    kappas = max(ROOTVSMALL, kappas);
 }
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
index 840a2da729..dbe23ae859 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
@@ -129,7 +129,8 @@ Foam::refinementFeatures::refinementFeatures
         }
 
         Info<< "Detected " << featurePoints.size()
-            << " featurePoints out of " << points.size() << endl;
+            << " featurePoints out of " << points.size()
+            << " on feature " << eMesh.name() << endl;
 
         pointTrees_.set
         (
@@ -164,26 +165,30 @@ void Foam::refinementFeatures::findNearestEdge
     forAll(edgeTrees_, featI)
     {
         const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
-        forAll(samples, sampleI)
+
+        if (tree.shapes().size() > 0)
         {
-            const point& sample = samples[sampleI];
-
-            scalar distSqr;
-            if (nearInfo[sampleI].hit())
+            forAll(samples, sampleI)
             {
-                distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
-            }
-            else
-            {
-                distSqr = nearestDistSqr[sampleI];
-            }
+                const point& sample = samples[sampleI];
 
-            pointIndexHit info = tree.findNearest(sample, distSqr);
+                scalar distSqr;
+                if (nearInfo[sampleI].hit())
+                {
+                    distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
+                }
+                else
+                {
+                    distSqr = nearestDistSqr[sampleI];
+                }
 
-            if (info.hit())
-            {
-                nearInfo[sampleI] = info;
-                nearFeature[sampleI] = featI;
+                pointIndexHit info = tree.findNearest(sample, distSqr);
+
+                if (info.hit())
+                {
+                    nearInfo[sampleI] = info;
+                    nearFeature[sampleI] = featI;
+                }
             }
         }
     }
@@ -206,33 +211,37 @@ void Foam::refinementFeatures::findNearestPoint
     forAll(pointTrees_, featI)
     {
         const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
-        forAll(samples, sampleI)
+
+        if (tree.shapes().pointLabels().size() > 0)
         {
-            const point& sample = samples[sampleI];
-
-            scalar distSqr;
-            if (nearFeature[sampleI] != -1)
+            forAll(samples, sampleI)
             {
-                label nearFeatI = nearFeature[sampleI];
-                const indexedOctree<treeDataPoint>& nearTree =
-                    pointTrees_[nearFeatI];
-                label featPointI =
-                    nearTree.shapes().pointLabels()[nearIndex[sampleI]];
-                const point& featPt =
-                    operator[](nearFeatI).points()[featPointI];
-                distSqr = magSqr(featPt-sample);
-            }
-            else
-            {
-                distSqr = nearestDistSqr[sampleI];
-            }
+                const point& sample = samples[sampleI];
 
-            pointIndexHit info = tree.findNearest(sample, distSqr);
+                scalar distSqr;
+                if (nearFeature[sampleI] != -1)
+                {
+                    label nearFeatI = nearFeature[sampleI];
+                    const indexedOctree<treeDataPoint>& nearTree =
+                        pointTrees_[nearFeatI];
+                    label featPointI =
+                        nearTree.shapes().pointLabels()[nearIndex[sampleI]];
+                    const point& featPt =
+                        operator[](nearFeatI).points()[featPointI];
+                    distSqr = magSqr(featPt-sample);
+                }
+                else
+                {
+                    distSqr = nearestDistSqr[sampleI];
+                }
 
-            if (info.hit())
-            {
-                nearFeature[sampleI] = featI;
-                nearIndex[sampleI] = info.index();
+                pointIndexHit info = tree.findNearest(sample, distSqr);
+
+                if (info.hit())
+                {
+                    nearFeature[sampleI] = featI;
+                    nearIndex[sampleI] = info.index();
+                }
             }
         }
     }
diff --git a/src/regionModels/regionModel/regionModel1D/regionModel1D.C b/src/regionModels/regionModel/regionModel1D/regionModel1D.C
index 19bf9f2ece..033db12366 100644
--- a/src/regionModels/regionModel/regionModel1D/regionModel1D.C
+++ b/src/regionModels/regionModel/regionModel1D/regionModel1D.C
@@ -114,6 +114,7 @@ void Foam::regionModels::regionModel1D::initialise()
             boundaryFaceCells_[localPyrolysisFaceI].transfer(cellIDs);
 
             localPyrolysisFaceI++;
+            nLayers_ = nCells;
         }
     }
 
@@ -266,6 +267,7 @@ Foam::regionModels::regionModel1D::regionModel1D(const fvMesh& mesh)
     boundaryFaceFaces_(),
     boundaryFaceCells_(),
     boundaryFaceOppositeFace_(),
+    nLayers_(0),
     nMagSfPtr_(NULL),
     moveMesh_(false)
 {}
@@ -283,6 +285,7 @@ Foam::regionModels::regionModel1D::regionModel1D
     boundaryFaceFaces_(regionMesh().nCells()),
     boundaryFaceCells_(regionMesh().nCells()),
     boundaryFaceOppositeFace_(regionMesh().nCells()),
+    nLayers_(0),
     nMagSfPtr_(NULL),
     moveMesh_(true)
 {
diff --git a/src/regionModels/regionModel/regionModel1D/regionModel1D.H b/src/regionModels/regionModel/regionModel1D/regionModel1D.H
index 650dea510f..55812aab2d 100644
--- a/src/regionModels/regionModel/regionModel1D/regionModel1D.H
+++ b/src/regionModels/regionModel/regionModel1D/regionModel1D.H
@@ -88,6 +88,9 @@ protected:
             //- Global boundary face IDs oppossite coupled patch
             labelList boundaryFaceOppositeFace_;
 
+            //- Number of layers in the region
+            label nLayers_;
+
 
         // Geometry
 
@@ -157,6 +160,9 @@ public:
 
                 //- Return the face area magnitudes / [m2]
                 inline const surfaceScalarField& nMagSf() const;
+
+                //- Return the number of layers in the region
+                inline label nLayers() const;
 };
 
 
diff --git a/src/regionModels/regionModel/regionModel1D/regionModel1DI.H b/src/regionModels/regionModel/regionModel1D/regionModel1DI.H
index 55e61d3eef..a115fdf456 100644
--- a/src/regionModels/regionModel/regionModel1D/regionModel1DI.H
+++ b/src/regionModels/regionModel/regionModel1D/regionModel1DI.H
@@ -65,4 +65,10 @@ Foam::regionModels::regionModel1D::nMagSf() const
 }
 
 
+inline Foam::label Foam::regionModels::regionModel1D::nLayers() const
+{
+    return nLayers_;
+}
+
+
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSolution
index 5fdb36dac6..e39ae484a2 100644
--- a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSolution
+++ b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSolution
@@ -63,7 +63,7 @@ PIMPLE
         "(U|k|epsilon)"
         {
             relTol          0;
-            absTol          0.0001;
+            tolerance       0.0001;
         }
     }
 }