diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index 9b09d73bbe..858d83749e 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    rhoSimpleFoam
+    rhoPimpleFoam
 
 Description
     Transient solver for turbulent flow of compressible fluids for
diff --git a/applications/solvers/incompressible/pimpleFoam/Make/files b/applications/solvers/incompressible/pimpleFoam/Make/files
index 92943d3370..f4ec4391df 100644
--- a/applications/solvers/incompressible/pimpleFoam/Make/files
+++ b/applications/solvers/incompressible/pimpleFoam/Make/files
@@ -1,3 +1,3 @@
 pimpleFoam.C
 
-EXE = $(FOAM_USER_APPBIN)/pimpleFoam
+EXE = $(FOAM_APPBIN)/pimpleFoam
diff --git a/applications/test/speed/Make/options b/applications/test/speed/Make/options
index e69de29bb2..e680561987 100644
--- a/applications/test/speed/Make/options
+++ b/applications/test/speed/Make/options
@@ -0,0 +1 @@
+EXE_INC = /* -ffast-math -mtune=core2 */
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 1df11e06da..4737f38e53 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -370,7 +370,7 @@ int main(int argc, char *argv[])
         // Refinement parameters
         refinementParameters refineParams(refineDict);
 
-        refineDriver.doRefine(refineDict, refineParams, wantSnap);
+        refineDriver.doRefine(refineDict, refineParams, wantSnap, motionDict);
 
         writeMesh
         (
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml
index 93bb9ffe9d..644b794a44 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml
@@ -27,6 +27,92 @@
       </Documentation>
     </DoubleVectorProperty>
 
+
+    <!-- Global settings -->
+
+    <!-- Extrapolate Walls check-box -->
+    <IntVectorProperty
+      name="ExtrapolateWalls"
+      command="SetExtrapolateWalls"
+      number_of_elements="1"
+      default_values="0"
+      animateable="0">
+      <BooleanDomain name="bool"/>
+      <Documentation>
+        Extrapolate internalField to wall and empty patches
+      </Documentation>
+    </IntVectorProperty>
+
+    <!-- Include Sets check-box -->
+    <IntVectorProperty
+      name="IncludeSets"
+      command="SetIncludeSets"
+      number_of_elements="1"
+      default_values="0"
+      animateable="0">
+      <Documentation>
+        Searches the polyMesh/sets/ directory
+      </Documentation>
+      <BooleanDomain name="bool"/>
+    </IntVectorProperty>
+
+    <!-- Include Zones check-box -->
+    <IntVectorProperty
+      name="IncludeZones"
+      command="SetIncludeZones"
+      number_of_elements="1"
+      default_values="0"
+      animateable="0">
+      <Documentation>
+        ZoneMesh information is used to find {cell,face,point}Zones.
+        The polyMesh/ directory is only checked on startup.
+      </Documentation>
+      <BooleanDomain name="bool"/>
+    </IntVectorProperty>
+
+    <!-- Show Patch Names check-box -->
+    <IntVectorProperty
+      name="ShowPatchNames"
+      command="SetShowPatchNames"
+      number_of_elements="1"
+      default_values="0"
+      animateable="0">
+      <BooleanDomain name="bool"/>
+      <Documentation>
+        Show patch names in render window.
+      </Documentation>
+    </IntVectorProperty>
+
+    <!-- Cache Mesh check-box -->
+    <IntVectorProperty
+      name="CacheMesh"
+      command="SetCacheMesh"
+      number_of_elements="1"
+      default_values="1"
+      animateable="0">
+      <BooleanDomain name="bool"/>
+      <Documentation>
+        Cache the fvMesh in memory.
+      </Documentation>
+    </IntVectorProperty>
+
+
+    <!-- Update GUI check box -->
+    <IntVectorProperty
+      name="UpdateGUI"
+      command="SetUpdateGUI"
+      number_of_elements="1"
+      default_values="0"
+      animateable="0">
+      <BooleanDomain name="bool"/>
+      <Documentation>
+        A simple way cause a reader GUI modification.
+      </Documentation>
+    </IntVectorProperty>
+
+
+    <!-- Selections -->
+
     <!-- Available Parts (volume, patches, lagrangian) array -->
     <StringVectorProperty
       name="PartArrayInfo"
@@ -118,87 +204,6 @@
         </RequiredProperties>
       </ArraySelectionDomain>
     </StringVectorProperty>
-
-    <!-- Extrapolate Walls check-box -->
-    <IntVectorProperty
-      name="ExtrapolateWalls"
-      command="SetExtrapolateWalls"
-      number_of_elements="1"
-      default_values="0"
-      animateable="0">
-      <BooleanDomain name="bool"/>
-      <Documentation>
-        Extrapolate internalField to wall and empty patches
-      </Documentation>
-    </IntVectorProperty>
-
-    <!-- Include Sets check-box -->
-    <IntVectorProperty
-      name="IncludeSets"
-      command="SetIncludeSets"
-      number_of_elements="1"
-      default_values="0"
-      animateable="0">
-      <Documentation>
-        Searches the polyMesh/sets/ directory
-      </Documentation>
-      <BooleanDomain name="bool"/>
-    </IntVectorProperty>
-
-    <!-- Include Zones check-box -->
-    <IntVectorProperty
-      name="IncludeZones"
-      command="SetIncludeZones"
-      number_of_elements="1"
-      default_values="0"
-      animateable="0">
-      <Documentation>
-        ZoneMesh information is used to find {cell,face,point}Zones.
-        The polyMesh/ directory is only checked on startup.
-      </Documentation>
-      <BooleanDomain name="bool"/>
-    </IntVectorProperty>
-
-    <!-- Show Patch Names check-box -->
-    <IntVectorProperty
-      name="ShowPatchNames"
-      command="SetShowPatchNames"
-      number_of_elements="1"
-      default_values="0"
-      animateable="0">
-      <BooleanDomain name="bool"/>
-      <Documentation>
-        Show patch names in render window.
-      </Documentation>
-    </IntVectorProperty>
-
-    <!-- Cache Mesh check-box -->
-    <IntVectorProperty
-      name="CacheMesh"
-      command="SetCacheMesh"
-      number_of_elements="1"
-      default_values="1"
-      animateable="0">
-      <BooleanDomain name="bool"/>
-      <Documentation>
-        Cache the fvMesh in memory.
-      </Documentation>
-    </IntVectorProperty>
-
-
-    <!-- Update GUI check box -->
-    <IntVectorProperty
-      name="UpdateGUI"
-      command="SetUpdateGUI"
-      number_of_elements="1"
-      default_values="0"
-      animateable="0">
-      <BooleanDomain name="bool"/>
-      <Documentation>
-        A simple way cause a reader GUI modification.
-      </Documentation>
-    </IntVectorProperty>
-
   </SourceProxy>
   </ProxyGroup>
 </ServerManagerConfiguration>
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
index 2a78ad0133..8fc518d745 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
@@ -623,6 +623,7 @@ void vtkPV3FoamReader::SetLagrangianFieldArrayStatus
     }
 }
 
+
 // ----------------------------------------------------------------------
 
 void vtkPV3FoamReader::SelectionModifiedCallback
@@ -658,4 +659,5 @@ int vtkPV3FoamReader::FillOutputPortInformation
     return 1;
 }
 
+
 // ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
index 48f52d7880..c07b950d0a 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
@@ -491,6 +491,8 @@ void Foam::autoHexMeshDriver::doMesh()
 
     if (wantRefine)
     {
+        const dictionary& motionDict = dict_.subDict("motionDict");
+
         autoRefineDriver refineDriver
         (
             meshRefinerPtr_(),
@@ -502,7 +504,7 @@ void Foam::autoHexMeshDriver::doMesh()
         // Get all the refinement specific params
         refinementParameters refineParams(dict_, 1);
 
-        refineDriver.doRefine(dict_, refineParams, wantSnap);
+        refineDriver.doRefine(dict_, refineParams, wantSnap, motionDict);
 
         // Write mesh
         writeMesh("Refined mesh");
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index 7bbe58ff8f..ca0c0cbddc 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -497,7 +497,8 @@ Foam::label Foam::autoRefineDriver::shellRefine
 void Foam::autoRefineDriver::baffleAndSplitMesh
 (
     const refinementParameters& refineParams,
-    const bool handleSnapProblems
+    const bool handleSnapProblems,
+    const dictionary& motionDict
 )
 {
     Info<< nl
@@ -514,6 +515,7 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
     (
         handleSnapProblems,
         !handleSnapProblems,            // merge free standing baffles?
+        motionDict,
         const_cast<Time&>(mesh.time()),
         globalToPatch_,
         refineParams.keepPoints()[0]
@@ -568,7 +570,8 @@ void Foam::autoRefineDriver::zonify
 void Foam::autoRefineDriver::splitAndMergeBaffles
 (
     const refinementParameters& refineParams,
-    const bool handleSnapProblems
+    const bool handleSnapProblems,
+    const dictionary& motionDict
 )
 {
     Info<< nl
@@ -588,6 +591,7 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
     (
         handleSnapProblems,
         false,                  // merge free standing baffles?
+        motionDict,
         const_cast<Time&>(mesh.time()),
         globalToPatch_,
         refineParams.keepPoints()[0]
@@ -685,7 +689,8 @@ void Foam::autoRefineDriver::doRefine
 (
     const dictionary& refineDict,
     const refinementParameters& refineParams,
-    const bool prepareForSnapping
+    const bool prepareForSnapping,
+    const dictionary& motionDict
 )
 {
     Info<< nl
@@ -734,13 +739,13 @@ void Foam::autoRefineDriver::doRefine
     );
 
     // Introduce baffles at surface intersections
-    baffleAndSplitMesh(refineParams, prepareForSnapping);
+    baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);
 
     // Mesh is at its finest. Do optional zoning.
     zonify(refineParams);
 
     // Pull baffles apart
-    splitAndMergeBaffles(refineParams, prepareForSnapping);
+    splitAndMergeBaffles(refineParams, prepareForSnapping, motionDict);
 
     // Do something about cells with refined faces on the boundary
     if (prepareForSnapping)
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
index 8aa33855de..f1747a29be 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
@@ -46,7 +46,6 @@ namespace Foam
 class featureEdgeMesh;
 class refinementParameters;
 
-
 /*---------------------------------------------------------------------------*\
                            Class autoRefineDriver Declaration
 \*---------------------------------------------------------------------------*/
@@ -112,7 +111,8 @@ class autoRefineDriver
         void baffleAndSplitMesh
         (
             const refinementParameters& refineParams,
-            const bool handleSnapProblems
+            const bool handleSnapProblems,
+            const dictionary& motionDict
         );
 
         //- Add zones
@@ -121,7 +121,8 @@ class autoRefineDriver
         void splitAndMergeBaffles
         (
             const refinementParameters& refineParams,
-            const bool handleSnapProblems
+            const bool handleSnapProblems,
+            const dictionary& motionDict
         );
 
         //- Merge refined boundary faces (from exposing coarser cell)
@@ -163,7 +164,8 @@ public:
         (
             const dictionary& refineDict,
             const refinementParameters& refineParams,
-            const bool prepareForSnapping
+            const bool prepareForSnapping,
+            const dictionary& motionDict
         );
 };
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index cc866d13ca..520e6f69ef 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -1294,6 +1294,169 @@ void Foam::autoSnapDriver::scaleMesh
 }
 
 
+// After snapping: correct patching according to nearest surface.
+// Code is very similar to calcNearestSurface.
+// - calculate face-wise snap distance as max of point-wise
+// - calculate face-wise nearest surface point
+// - repatch face according to patch for surface point.
+Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
+(
+    const snapParameters& snapParams,
+    const labelList& adaptPatchIDs
+)
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
+
+    Info<< "Repatching faces according to nearest surface ..." << endl;
+
+    // Get the labels of added patches.
+    autoPtr<indirectPrimitivePatch> ppPtr
+    (
+        meshRefinement::makePatch
+        (
+            mesh,
+            adaptPatchIDs
+        )
+    );
+    indirectPrimitivePatch& pp = ppPtr();
+
+    // Divide surfaces into zoned and unzoned
+    labelList zonedSurfaces;
+    labelList unzonedSurfaces;
+    getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
+
+
+    // Faces that do not move
+    PackedList<1> isZonedFace(mesh.nFaces(), 0);
+    {
+        // 1. All faces on zoned surfaces
+        const wordList& faceZoneNames = surfaces.faceZoneNames();
+        const faceZoneMesh& fZones = mesh.faceZones();
+
+        forAll(zonedSurfaces, i)
+        {
+            label zoneSurfI = zonedSurfaces[i];
+
+            label zoneI = fZones.findZoneID(faceZoneNames[zoneSurfI]);
+
+            const faceZone& fZone = fZones[zoneI];
+
+            forAll(fZone, i)
+            {
+                isZonedFace.set(fZone[i], 1);
+            }  
+        }
+    }
+
+
+    // Determine per pp face which patch it should be in
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Patch that face should be in
+    labelList closestPatch(pp.size(), -1);
+    {
+        // face snap distance as max of point snap distance
+        scalarField faceSnapDist(pp.size(), -GREAT);
+        {
+            // Distance to attract to nearest feature on surface
+            const scalarField snapDist(calcSnapDistance(snapParams, pp));
+
+            const faceList& localFaces = pp.localFaces();
+
+            forAll(localFaces, faceI)
+            {
+                const face& f = localFaces[faceI];
+
+                forAll(f, fp)
+                {
+                    faceSnapDist[faceI] = max
+                    (
+                        faceSnapDist[faceI],
+                        snapDist[f[fp]]
+                    );
+                }
+            }
+        }
+
+        pointField localFaceCentres(pp.size());
+        forAll(pp, i)
+        {
+            localFaceCentres[i] = mesh.faceCentres()[pp.addressing()[i]];
+        }
+
+        // Get nearest surface and region
+        labelList hitSurface;
+        labelList hitRegion;
+        surfaces.findNearestRegion
+        (
+            unzonedSurfaces,
+            localFaceCentres,
+            sqr(4*faceSnapDist),    // sqr of attract distance
+            hitSurface,
+            hitRegion
+        );
+
+        // Get patch
+        forAll(pp, i)
+        {
+            label faceI = pp.addressing()[i];
+
+            if (hitSurface[i] != -1 && (isZonedFace.get(faceI) == 0))
+            {
+                closestPatch[i] = globalToPatch_
+                [
+                    surfaces.globalRegion
+                    (
+                        hitSurface[i],
+                        hitRegion[i]
+                    )
+                ];
+            }
+        }
+    }
+
+
+    // Change those faces for which there is a different closest patch
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList ownPatch(mesh.nFaces(), -1);
+    labelList neiPatch(mesh.nFaces(), -1);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        forAll(pp, i)
+        {
+            ownPatch[pp.start()+i] = patchI;
+            neiPatch[pp.start()+i] = patchI;
+        }
+    }
+
+    label nChanged = 0;
+    forAll(closestPatch, i)
+    {
+        label faceI = pp.addressing()[i];
+
+        if (closestPatch[i] != -1 && closestPatch[i] != ownPatch[faceI])
+        {
+            ownPatch[faceI] = closestPatch[i];
+            neiPatch[faceI] = closestPatch[i];
+            nChanged++;
+        }
+    }
+
+    Info<< "Repatched " << returnReduce(nChanged, sumOp<label>())
+        << " faces in = " << mesh.time().cpuTimeIncrement() << " s\n" << nl
+        << endl;
+
+    return meshRefiner_.createBaffles(ownPatch, neiPatch);
+}
+
+
 void Foam::autoSnapDriver::doSnap
 (
     const dictionary& snapDict,
@@ -1329,7 +1492,7 @@ void Foam::autoSnapDriver::doSnap
         );
         indirectPrimitivePatch& pp = ppPtr();
 
-        // Distance to attact to nearest feature on surface
+        // Distance to attract to nearest feature on surface
         const scalarField snapDist(calcSnapDistance(snapParams, pp));
 
 
@@ -1383,6 +1546,9 @@ void Foam::autoSnapDriver::doSnap
 
     // Merge any introduced baffles.
     mergeZoneBaffles(baffles);
+
+    // Repatch faces according to nearest.
+    repatchToSurface(snapParams, adaptPatchIDs);
 }
 
 
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 82cf448b98..1604e35859 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -215,6 +215,14 @@ public:
                 motionSmoother&
             );
 
+            //- Repatch faces according to surface nearest the face centre
+            autoPtr<mapPolyMesh> repatchToSurface
+            (
+                const snapParameters& snapParams,
+                const labelList& adaptPatchIDs
+            );
+
+
             void doSnap
             (
                 const dictionary& snapDict,
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 4fe3005766..7f5bc29a79 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -50,6 +50,7 @@ License
 #include "globalIndex.H"
 #include "meshTools.H"
 #include "OFstream.H"
+#include "geomDecomp.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -426,15 +427,14 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
 
 // Determine for multi-processor regions the lowest numbered cell on the lowest
 // numbered processor.
-void Foam::meshRefinement::getRegionMaster
+void Foam::meshRefinement::getCoupledRegionMaster
 (
+    const globalIndex& globalCells,
     const boolList& blockedFace,
     const regionSplit& globalRegion,
     Map<label>& regionToMaster
 ) const
 {
-    globalIndex globalCells(mesh_.nCells());
-
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
     forAll(patches, patchI)
@@ -483,6 +483,274 @@ void Foam::meshRefinement::getRegionMaster
 }
 
 
+void Foam::meshRefinement::calcLocalRegions
+(
+    const globalIndex& globalCells,
+    const labelList& globalRegion,
+    const Map<label>& coupledRegionToMaster,
+
+    Map<label>& globalToLocalRegion,
+    pointField& localPoints
+) const
+{
+    globalToLocalRegion.resize(globalRegion.size());
+    DynamicList<point> localCc(globalRegion.size()/2);
+
+    forAll(globalRegion, cellI)
+    {
+        Map<label>::const_iterator fndMaster =
+            coupledRegionToMaster.find(globalRegion[cellI]);
+
+        if (fndMaster != coupledRegionToMaster.end())
+        {
+            // Multi-processor region.
+            if (globalCells.toGlobal(cellI) == fndMaster())
+            {
+                // I am master. Allocate region for me.
+                globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
+                localCc.append(mesh_.cellCentres()[cellI]);
+            }
+        }
+        else
+        {
+            // Single processor region.
+            if (globalToLocalRegion.insert(globalRegion[cellI], localCc.size()))
+            {
+                localCc.append(mesh_.cellCentres()[cellI]);
+            }
+        }
+    }
+    localCc.shrink();
+    localPoints.transfer(localCc);
+
+    if (localPoints.size() != globalToLocalRegion.size())
+    {
+        FatalErrorIn("calcLocalRegions(..)")
+            << "localPoints:" << localPoints.size()
+            << " globalToLocalRegion:" << globalToLocalRegion.size()
+            << exit(FatalError);
+    }
+}
+
+
+Foam::label Foam::meshRefinement::getShiftedRegion
+(
+    const globalIndex& indexer,
+    const Map<label>& globalToLocalRegion,
+    const Map<label>& coupledRegionToShifted,
+    const label globalRegion
+)
+{
+    Map<label>::const_iterator iter =
+        globalToLocalRegion.find(globalRegion);
+
+    if (iter != globalToLocalRegion.end())
+    {
+        // Region is 'owned' locally. Convert local region index into global.
+        return indexer.toGlobal(iter());
+    }
+    else
+    {
+        return coupledRegionToShifted[globalRegion];
+    }
+}
+
+
+// Add if not yet present
+void Foam::meshRefinement::addUnique(const label elem, labelList& lst)
+{
+    if (findIndex(lst, elem) == -1)
+    {
+        label sz = lst.size();
+        lst.setSize(sz+1);
+        lst[sz] = elem;
+    }
+}
+
+
+void Foam::meshRefinement::calcRegionRegions
+(
+    const labelList& globalRegion,
+    const Map<label>& globalToLocalRegion,
+    const Map<label>& coupledRegionToMaster,
+    labelListList& regionRegions
+) const
+{
+    // Global region indexing since we now know the shifted regions.
+    globalIndex shiftIndexer(globalToLocalRegion.size());
+
+    // Redo the coupledRegionToMaster to be in shifted region indexing.
+    Map<label> coupledRegionToShifted(coupledRegionToMaster.size());
+    forAllConstIter(Map<label>, coupledRegionToMaster, iter)
+    {
+        label region = iter.key();
+
+        Map<label>::const_iterator fndRegion = globalToLocalRegion.find(region);
+
+        if (fndRegion != globalToLocalRegion.end())
+        {
+            // A local cell is master of this region. Get its shifted region.
+            coupledRegionToShifted.insert
+            (
+                region,
+                shiftIndexer.toGlobal(fndRegion())
+            );
+        }
+        Pstream::mapCombineGather(coupledRegionToShifted, minEqOp<label>());
+        Pstream::mapCombineScatter(coupledRegionToShifted);
+    }
+
+
+    // Determine region-region connectivity.
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // This is for all my regions (so my local ones or the ones I am master
+    // of) the neighbouring region indices.
+
+
+    // Transfer lists.
+    PtrList<HashSet<edge, Hash<edge> > > regionConnectivity(Pstream::nProcs());
+    forAll(regionConnectivity, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            regionConnectivity.set
+            (
+                procI,
+                new HashSet<edge, Hash<edge> >
+                (
+                    coupledRegionToShifted.size()
+                  / Pstream::nProcs()
+                )
+            );
+        }
+    }
+
+
+    // Connectivity. For all my local regions the connected regions.
+    regionRegions.setSize(globalToLocalRegion.size());
+
+    // Add all local connectivity to regionRegions; add all non-local
+    // to the transferlists.
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        label ownRegion = globalRegion[mesh_.faceOwner()[faceI]];
+        label neiRegion = globalRegion[mesh_.faceNeighbour()[faceI]];
+
+        if (ownRegion != neiRegion)
+        {
+            label shiftOwnRegion = getShiftedRegion
+            (
+                shiftIndexer,
+                globalToLocalRegion,
+                coupledRegionToShifted,
+                ownRegion
+            );
+            label shiftNeiRegion = getShiftedRegion
+            (
+                shiftIndexer,
+                globalToLocalRegion,
+                coupledRegionToShifted,
+                neiRegion
+            );
+
+
+            // Connection between two regions. Send to owner of region.
+            // - is ownRegion 'owned' by me
+            // - is neiRegion 'owned' by me
+
+            if (shiftIndexer.isLocal(shiftOwnRegion))
+            {
+                label localI = shiftIndexer.toLocal(shiftOwnRegion);
+                addUnique(shiftNeiRegion, regionRegions[localI]);
+            }
+            else
+            {
+                label masterProc = shiftIndexer.whichProcID(shiftOwnRegion);
+                regionConnectivity[masterProc].insert
+                (
+                    edge(shiftOwnRegion, shiftNeiRegion)
+                );
+            }
+
+            if (shiftIndexer.isLocal(shiftNeiRegion))
+            {
+                label localI = shiftIndexer.toLocal(shiftNeiRegion);
+                addUnique(shiftOwnRegion, regionRegions[localI]);
+            }
+            else
+            {
+                label masterProc = shiftIndexer.whichProcID(shiftNeiRegion);
+                regionConnectivity[masterProc].insert
+                (
+                    edge(shiftOwnRegion, shiftNeiRegion)
+                );
+            }
+        }
+    }
+
+
+    // Send
+    forAll(regionConnectivity, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            OPstream str(Pstream::blocking, procI);
+            str << regionConnectivity[procI];
+        }
+    }
+    // Receive
+    forAll(regionConnectivity, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            IPstream str(Pstream::blocking, procI);
+            str >> regionConnectivity[procI];
+        }
+    }
+
+    // Add to addressing.
+    forAll(regionConnectivity, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            for
+            (
+                HashSet<edge, Hash<edge> >::const_iterator iter =
+                    regionConnectivity[procI].begin();
+                iter != regionConnectivity[procI].end();
+                ++iter
+            )
+            {
+                const edge& e = iter.key();
+
+                bool someLocal = false;
+                if (shiftIndexer.isLocal(e[0]))
+                {
+                    label localI = shiftIndexer.toLocal(e[0]);
+                    addUnique(e[1], regionRegions[localI]);
+                    someLocal = true;
+                }
+                if (shiftIndexer.isLocal(e[1]))
+                {
+                    label localI = shiftIndexer.toLocal(e[1]);
+                    addUnique(e[0], regionRegions[localI]);
+                    someLocal = true;
+                }
+
+                if (!someLocal)
+                {
+                    FatalErrorIn("calcRegionRegions(..)")
+                        << "Received from processor " << procI
+                        << " connection " << e
+                        << " where none of the elements is local to me."
+                        << abort(FatalError);
+                }
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
@@ -598,88 +866,88 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
     // the region might span multiple domains so we want to get
     // a "region master" per domain. Note that multi-processor
     // regions can only occur on cells on coupled patches.
-
-
-    // Determine per coupled region the master cell (lowest numbered cell
-    // on lowest numbered processor)
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    Map<label> regionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
-    getRegionMaster(blockedFace, globalRegion, regionToMaster);
+    // Note: since the number of regions does not change by this the
+    // process can be seen as just a shift of a region onto a single
+    // processor.
 
 
     // Global cell numbering engine
     globalIndex globalCells(mesh_.nCells());
 
 
-    // Determine cell centre per region
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Now we can divide regions into
-    // - single-processor regions (almost all)
-    //   - allocate local region + coordinate for it
-    // - multi-processor for which I am master
-    //   - allocate local region + coordinate for it
-    // - multi-processor for which I am not master
-    //   - do not allocate region.
-    //   - but inherit new distribution from master.
+    // Determine per coupled region the master cell (lowest numbered cell
+    // on lowest numbered processor)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // (does not determine master for non-coupled=fully-local regions)
 
-    Map<label> globalToLocalRegion(mesh_.nCells());
-    DynamicList<point> localCc(mesh_.nCells()/10);
+    Map<label> coupledRegionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
+    getCoupledRegionMaster
+    (
+        globalCells,
+        blockedFace,
+        globalRegion,
+        coupledRegionToMaster
+    );
 
-    forAll(globalRegion, cellI)
-    {
-        Map<label>::const_iterator fndMaster =
-            regionToMaster.find(globalRegion[cellI]);
+    // Determine my regions
+    // ~~~~~~~~~~~~~~~~~~~~
+    // These are the fully local ones or the coupled ones of which I am master.
 
-        if (fndMaster != regionToMaster.end())
-        {
-            // Multi-processor region.
-            if (globalCells.toGlobal(cellI) == fndMaster())
-            {
-                // I am master. Allocate region for me.
-                globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
-                localCc.append(mesh_.cellCentres()[cellI]);
-            }
-        }
-        else
-        {
-            // Single processor region.
-            Map<label>::iterator iter =
-                globalToLocalRegion.find(globalRegion[cellI]);
-
-            if (iter == globalToLocalRegion.end())
-            {
-                globalToLocalRegion.insert(globalRegion[cellI], localCc.size());
-                localCc.append(mesh_.cellCentres()[cellI]);
-            }
-        }
-    }
-    localCc.shrink();
+    Map<label> globalToLocalRegion;
     pointField localPoints;
-    localPoints.transfer(localCc);
+    calcLocalRegions
+    (
+        globalCells,
+        globalRegion,
+        coupledRegionToMaster,
+
+        globalToLocalRegion,
+        localPoints
+    );
 
 
-    // Call decomposer on localCc
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    labelList localDistribution = decomposer.decompose(localPoints);
+    // Find distribution for regions
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+    labelList regionDistribution;
 
-    // Distribute the destination processor for coupled regions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    Map<label> regionToDist(regionToMaster.size());
-
-    forAllConstIter(Map<label>, regionToMaster, iter)
+    if (isA<geomDecomp>(decomposer))
     {
-        if (globalCells.isLocal(iter()))
+        regionDistribution = decomposer.decompose(localPoints);
+    }
+    else
+    {
+        labelListList regionRegions;
+        calcRegionRegions
+        (
+            globalRegion,
+            globalToLocalRegion,
+            coupledRegionToMaster,
+
+            regionRegions
+        );
+
+        regionDistribution = decomposer.decompose(regionRegions, localPoints);
+    }
+
+
+
+    // Convert region-based decomposition back to cell-based one
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Transfer destination processor back to all. Use global reduce for now.
+    Map<label> regionToDist(coupledRegionToMaster.size());
+    forAllConstIter(Map<label>, coupledRegionToMaster, iter)
+    {
+        label region = iter.key();
+
+        Map<label>::const_iterator regionFnd = globalToLocalRegion.find(region);
+
+        if (regionFnd != globalToLocalRegion.end())
         {
-            // Master cell is local.
-            regionToDist.insert
-            (
-                iter.key(),
-                localDistribution[globalToLocalRegion[iter.key()]]
-            );
+            // Master cell is local. Store my distribution.
+            regionToDist.insert(iter.key(), regionDistribution[regionFnd()]);
         }
         else
         {
@@ -691,29 +959,27 @@ Foam::labelList Foam::meshRefinement::decomposeCombineRegions
     Pstream::mapCombineScatter(regionToDist);
 
 
-    // Convert region-based decomposition back to cell-based one
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
+    // Determine destination for all cells
     labelList distribution(mesh_.nCells());
 
     forAll(globalRegion, cellI)
     {
-        Map<label>::const_iterator fndMaster =
+        Map<label>::const_iterator fndRegion =
             regionToDist.find(globalRegion[cellI]);
 
-        if (fndMaster != regionToDist.end())
+        if (fndRegion != regionToDist.end())
         {
-            // Special handling
-            distribution[cellI] = fndMaster();
+            distribution[cellI] = fndRegion();
         }
         else
         {
             // region is local to the processor.
             label localRegionI = globalToLocalRegion[globalRegion[cellI]];
 
-            distribution[cellI] = localDistribution[localRegionI];
+            distribution[cellI] = regionDistribution[localRegionI];
         }
     }
+
     return distribution;
 }
 
@@ -839,9 +1105,18 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
         if (debug)
         {
-            Pout<< "Wanted distribution:"
-                << distributor.countCells(distribution)
-                << endl;
+            labelList nProcCells(distributor.countCells(distribution));
+            Pout<< "Wanted distribution:" << nProcCells << endl;
+
+            Pstream::listCombineGather(nProcCells, plusEqOp<label>());
+            Pstream::listCombineScatter(nProcCells);
+
+            Pout<< "Wanted resulting decomposition:" << endl;
+            forAll(nProcCells, procI)
+            {
+                Pout<< "    " << procI << '\t' << nProcCells[procI] << endl;
+            }
+            Pout<< endl;
         }
         // Do actual sending/receiving of mesh
         map = distributor.distribute(distribution);
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 012ef08d00..d3eb52b7c4 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -68,6 +68,7 @@ class featureEdgeMesh;
 class fvMeshDistribute;
 class searchableSurface;
 class regionSplit;
+class globalIndex;
 
 /*---------------------------------------------------------------------------*\
                            Class meshRefinement Declaration
@@ -157,15 +158,6 @@ private:
             labelList& neiLevel,
             pointField& neiCc
         ) const;
-        ////- Calculate coupled boundary end vector and refinement level. Sort
-        ////  so both sides of coupled face have data in same order.
-        //void calcCanonicalBoundaryData
-        //(
-        //    labelList& leftLevel,
-        //    pointField& leftCc,
-        //    labelList& rightLevel,
-        //    pointField& rightCc
-        //)  const;
 
         //- Find any intersection of surface. Store in surfaceIndex_.
         void updateIntersections(const labelList& changedFaces);
@@ -194,15 +186,51 @@ private:
             const label exposedPatchI
         );
 
-        //- Used in decomposeCombineRegions. Given global region per cell
-        //  determines master processor/cell for regions straddling
-        //  procboundaries.
-        void getRegionMaster
-        (
-            const boolList& blockedFace,
-            const regionSplit& globalRegion,
-            Map<label>& regionToMaster
-        ) const;
+        // For decomposeCombineRegions
+
+            //- Used in decomposeCombineRegions. Given global region per cell
+            //  determines master processor/cell for regions straddling
+            //  procboundaries.
+            void getCoupledRegionMaster
+            (
+                const globalIndex& globalCells,
+                const boolList& blockedFace,
+                const regionSplit& globalRegion,
+                Map<label>& regionToMaster
+            ) const;
+
+            //- Determine regions that are local to me or coupled ones that
+            //  are owned by me. Determine representative location.
+            void calcLocalRegions
+            (
+                const globalIndex& globalCells,
+                const labelList& globalRegion,
+                const Map<label>& coupledRegionToMaster,
+
+                Map<label>& globalToLocalRegion,
+                pointField& localPoints
+            ) const;
+
+            //- Convert region into global index.
+            static label getShiftedRegion
+            (
+                const globalIndex& indexer,
+                const Map<label>& globalToLocalRegion,
+                const Map<label>& coupledRegionToShifted,
+                const label globalRegion
+            );
+
+            //- helper: add element if not in list. Linear search.
+            static void addUnique(const label, labelList&);
+
+            //- Calculate region connectivity. Major communication.
+            void calcRegionRegions
+            (
+                const labelList& globalRegion,
+                const Map<label>& globalToLocalRegion,
+                const Map<label>& coupledRegionToMaster,
+                labelListList& regionRegions
+            ) const;
 
 
         // Refinement candidate selection
@@ -334,6 +362,14 @@ private:
                 const labelList& globalToPatch
             ) const;
 
+            //- Returns list with for every internal face -1 or the patch
+            //  they should be baffled into.
+            labelList markFacesOnProblemCellsGeometric
+            (
+                const dictionary& motionDict,
+                const labelList& globalToPatch
+            ) const;
+
 
         // Baffle merging
 
@@ -552,6 +588,7 @@ public:
             (
                 const bool handleSnapProblems,
                 const bool mergeFreeStanding,
+                const dictionary& motionDict,
                 Time& runTime,
                 const labelList& globalToPatch,
                 const point& keepPoint
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 8b0236a7e2..e6bd37f670 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -44,9 +44,9 @@ License
 #include "OFstream.H"
 #include "regionSplit.H"
 #include "removeCells.H"
-
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+#include "motionSmoother.H"
+#include "polyMeshGeometry.H"
+#include "IOmanip.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -555,28 +555,28 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     boolList isBoundaryFace(mesh_.nFaces(), false);
 
     // Fill boundary data. All elements on meshed patches get marked.
-    forAll(globalToPatch, i)
+    // Get the labels of added patches.
+    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+    forAll(adaptPatchIDs, i)
     {
-        label patchI = globalToPatch[i];
+        label patchI = adaptPatchIDs[i];
 
-        if (patchI != -1)
+        const polyPatch& pp = patches[patchI];
+
+        label faceI = pp.start();
+
+        forAll(pp, j)
         {
-            const polyPatch& pp = patches[patchI];
+            markBoundaryFace
+            (
+                faceI,
+                isBoundaryFace,
+                isBoundaryEdge,
+                isBoundaryPoint
+            );
 
-            label faceI = pp.start();
-
-            forAll(pp, j)
-            {
-                markBoundaryFace
-                (
-                    faceI,
-                    isBoundaryFace,
-                    isBoundaryEdge,
-                    isBoundaryPoint
-                );
-
-                faceI++;
-            }
+            faceI++;
         }
     }
 
@@ -872,6 +872,264 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
     return facePatch;
 }
+//XXXXXXXXXXXXXX
+// Mark faces to be baffled to prevent snapping problems. Does
+// test to find nearest surface and checks which faces would get squashed.
+Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
+(
+    const dictionary& motionDict,
+    const labelList& globalToPatch
+) const
+{
+    // Get the labels of added patches.
+    labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+    // Construct addressing engine.
+    autoPtr<indirectPrimitivePatch> ppPtr
+    (
+        meshRefinement::makePatch
+        (
+            mesh_,
+            adaptPatchIDs
+        )
+    );
+    const indirectPrimitivePatch& pp = ppPtr();
+    const pointField& localPoints = pp.localPoints();
+    const labelList& meshPoints = pp.meshPoints();
+
+    // Find nearest (non-baffle) surface
+    pointField newPoints(mesh_.points());
+    {
+        List<pointIndexHit> hitInfo;
+        labelList hitSurface;
+        surfaces_.findNearest
+        (
+            surfaces_.getUnnamedSurfaces(),
+            localPoints,
+            scalarField(localPoints.size(), sqr(GREAT)),    // sqr of attraction
+            hitSurface,
+            hitInfo
+        );
+    
+        forAll(hitInfo, i)
+        {
+            if (hitInfo[i].hit())
+            {
+                //label pointI = meshPoints[i];
+                //Pout<< "   " << pointI << " moved from "
+                //    << mesh_.points()[pointI] << " by "
+                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
+                //    << endl;
+                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
+            }
+        }
+    }
+
+    // Per face (internal or coupled!) the patch that the
+    // baffle should get (or -1).
+    labelList facePatch(mesh_.nFaces(), -1);
+    // Count of baffled faces
+    label nBaffleFaces = 0;
+
+
+//    // Sync position? Or not since same face on both side so just sync
+//    // result of baffle.
+//
+//    const scalar minArea(readScalar(motionDict.lookup("minArea")));
+//
+//    Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
+//        << minArea << endl;
+//
+//    pointField facePoints;
+//    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+//    {
+//        const face& f = mesh_.faces()[faceI];
+//
+//        bool usesPatchPoint = false;
+//
+//        facePoints.setSize(f.size());
+//        forAll(f, fp)
+//        {
+//            Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
+//
+//            if (iter != pp.meshPointMap().end())
+//            {
+//                facePoints[fp] = newPosition[iter()];
+//                usesPatchPoint = true;
+//            }
+//            else
+//            {
+//                facePoints[fp] = mesh_.points()[f[fp]];
+//            }
+//        }
+//
+//        if (usesPatchPoint)
+//        {
+//            // Check area of face wrt original area
+//            face identFace(identity(f.size()));
+//
+//            if (identFace.mag(facePoints) < minArea)
+//            {
+//                facePatch[faceI] = getBafflePatch(facePatch, faceI);
+//                nBaffleFaces++;
+//            }
+//        }
+//    }
+//
+//
+//    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+//    forAll(patches, patchI)
+//    {
+//        const polyPatch& pp = patches[patchI];
+//
+//        if (pp.coupled())
+//        {
+//            forAll(pp, i)
+//            {
+//                label faceI = pp.start()+i;
+//
+//                const face& f = mesh_.faces()[faceI];
+//
+//                bool usesPatchPoint = false;
+//
+//                facePoints.setSize(f.size());
+//                forAll(f, fp)
+//                {
+//                    Map<label>::const_iterator iter =
+//                        pp.meshPointMap().find(f[fp]);
+//
+//                    if (iter != pp.meshPointMap().end())
+//                    {
+//                        facePoints[fp] = newPosition[iter()];
+//                        usesPatchPoint = true;
+//                    }
+//                    else
+//                    {
+//                        facePoints[fp] = mesh_.points()[f[fp]];
+//                    }
+//                }
+//
+//                if (usesPatchPoint)
+//                {
+//                    // Check area of face wrt original area
+//                    face identFace(identity(f.size()));
+//
+//                    if (identFace.mag(facePoints) < minArea)
+//                    {
+//                        facePatch[faceI] = getBafflePatch(facePatch, faceI);
+//                        nBaffleFaces++;
+//                    }
+//                }
+//            }
+//        }
+//    }
+
+    {
+        pointField oldPoints(mesh_.points());
+        mesh_.movePoints(newPoints);
+        faceSet wrongFaces(mesh_, "wrongFaces", 100);
+        {
+            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
+
+            // Just check the errors from squashing
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            const labelList allFaces(identity(mesh_.nFaces()));
+            label nWrongFaces = 0;
+
+            scalar minArea(readScalar(motionDict.lookup("minArea")));
+            if (minArea > -SMALL)
+            {
+                polyMeshGeometry::checkFaceArea
+                (
+                    false,
+                    minArea,
+                    mesh_,
+                    mesh_.faceAreas(),
+                    allFaces,
+                    &wrongFaces
+                );
+
+                label nNewWrongFaces = returnReduce
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces with area < "
+                    << setw(5) << minArea
+                    << " m^2                            : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+
+//            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
+            scalar minDet = 0.01;
+            if (minDet > -1)
+            {
+                polyMeshGeometry::checkCellDeterminant
+                (
+                    false,
+                    minDet,
+                    mesh_,
+                    mesh_.faceAreas(),
+                    allFaces,
+                    polyMeshGeometry::affectedCells(mesh_, allFaces),
+                    &wrongFaces
+                );
+
+                label nNewWrongFaces = returnReduce
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces on cells with determinant < "
+                    << setw(5) << minDet << "                : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+        }
+
+
+        forAllConstIter(faceSet, wrongFaces, iter)
+        {
+            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
+
+            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
+            {
+                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
+                nBaffleFaces++;
+
+                //Pout<< "    " << iter.key()
+                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
+                //    << " is destined for patch " << facePatch[iter.key()]
+                //    << endl;
+            }
+        }
+        // Restore points.
+        mesh_.movePoints(oldPoints);
+    }
+
+
+    Info<< "markFacesOnProblemCellsGeometric : marked "
+        << returnReduce(nBaffleFaces, sumOp<label>())
+        << " additional internal and coupled faces"
+        << " to be converted into baffles." << endl;
+
+    syncTools::syncFaceList
+    (
+        mesh_,
+        facePatch,
+        maxEqOp<label>(),
+        false               // no separation
+    );
+
+    return facePatch;
+}
+//XXXXXXXX
 
 
 // Return a list of coupled face pairs, i.e. faces that use the same vertices.
@@ -1554,6 +1812,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
 (
     const bool handleSnapProblems,
     const bool mergeFreeStanding,
+    const dictionary& motionDict,
     Time& runTime,
     const labelList& globalToPatch,
     const point& keepPoint
@@ -1584,6 +1843,12 @@ void Foam::meshRefinement::baffleAndSplitMesh
         ownPatch,
         neiPatch
     );
+
+    if (debug)
+    {
+        runTime++;
+    }
+
     createBaffles(ownPatch, neiPatch);
 
     if (debug)
@@ -1621,14 +1886,65 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
         labelList facePatch
         (
-            markFacesOnProblemCells
+            //markFacesOnProblemCells
+            //(
+            //    globalToPatch
+            //)
+            markFacesOnProblemCellsGeometric
             (
+                motionDict,
                 globalToPatch
             )
         );
         Info<< "Analyzed problem cells in = "
             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
 
+        if (debug)
+        {
+            // Dump all these faces to a faceSet.
+            faceSet problemGeom(mesh_, "problemFacesGeom", 100);
+
+            const labelList facePatchGeom
+            (
+                markFacesOnProblemCellsGeometric
+                (
+                    motionDict,
+                    globalToPatch
+                )
+            );
+            forAll(facePatchGeom, faceI)
+            {
+                if (facePatchGeom[faceI] != -1)
+                {
+                    problemGeom.insert(faceI);
+                }
+            }
+            Pout<< "Dumping " << problemGeom.size()
+                << " problem faces to " << problemGeom.objectPath() << endl;
+            problemGeom.write();
+
+
+            faceSet problemTopo(mesh_, "problemFacesTopo", 100);
+
+            const labelList facePatchTopo
+            (    
+                markFacesOnProblemCells
+                (
+                    globalToPatch
+                )
+            );
+            forAll(facePatchTopo, faceI)
+            {
+                if (facePatchTopo[faceI] != -1)
+                {
+                    problemTopo.insert(faceI);
+                }
+            }
+            Pout<< "Dumping " << problemTopo.size()
+                << " problem faces to " << problemTopo.objectPath() << endl;
+            problemTopo.write();
+        }
+
         Info<< "Introducing baffles to delete problem cells." << nl << endl;
 
         if (debug)
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index b8ee73b8ce..e9eb2cb445 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -1277,22 +1277,13 @@ Foam::autoPtr<Foam::mapDistributePolyMesh>
 
     if (Pstream::nProcs() > 1)
     {
-        labelList distribution(decomposer.decompose(mesh_.cellCentres()));
-        // Get distribution such that baffle faces stay internal to the
-        // processor.
-        //labelList distribution(decomposePreserveBaffles(decomposer));
-
-        if (debug)
-        {
-            Pout<< "Wanted distribution:"
-                << distributor.countCells(distribution)
-                << endl;
-        }
-        // Do actual sending/receiving of mesh
-        distMap = distributor.distribute(distribution);
-
-        // Update numbering of meshRefiner
-        distribute(distMap);
+        distMap = balance
+        (
+            false,  //keepZoneFaces
+            false,  //keepBaffles
+            decomposer,
+            distributor
+        );
 
         Info<< "Balanced mesh in = "
             << mesh_.time().cpuTimeIncrement() << " s" << endl;
@@ -1302,7 +1293,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh>
 
         if (debug)
         {
-            Pout<< "Writing " << msg
+            Pout<< "Writing balanced " << msg
                 << " mesh to time " << mesh_.time().timeName() << endl;
             write
             (
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index ce09e597e2..1654e3cd7f 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -377,15 +377,34 @@ Foam::refinementSurfaces::refinementSurfaces
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// Get indices of named surfaces (surfaces with cellZoneName)
+// Get indices of unnamed surfaces (surfaces without faceZoneName)
+Foam::labelList Foam::refinementSurfaces::getUnnamedSurfaces() const
+{
+   labelList anonymousSurfaces(faceZoneNames_.size());
+
+    label i = 0;
+    forAll(faceZoneNames_, surfI)
+    {
+        if (faceZoneNames_[surfI].size() == 0)
+        {
+            anonymousSurfaces[i++] = surfI;
+        }
+    }
+    anonymousSurfaces.setSize(i);
+
+    return anonymousSurfaces;
+}
+
+
+// Get indices of named surfaces (surfaces with faceZoneName)
 Foam::labelList Foam::refinementSurfaces::getNamedSurfaces() const
 {
-   labelList namedSurfaces(cellZoneNames_.size());
+   labelList namedSurfaces(faceZoneNames_.size());
 
     label namedI = 0;
-    forAll(cellZoneNames_, surfI)
+    forAll(faceZoneNames_, surfI)
     {
-        if (cellZoneNames_[surfI].size() > 0)
+        if (faceZoneNames_[surfI].size() > 0)
         {
             namedSurfaces[namedI++] = surfI;
         }
@@ -846,6 +865,69 @@ void Foam::refinementSurfaces::findNearest
 }
 
 
+void Foam::refinementSurfaces::findNearestRegion
+(
+    const labelList& surfacesToTest,
+    const pointField& samples,
+    const  scalarField& nearestDistSqr,
+    labelList& hitSurface,
+    labelList& hitRegion
+) const
+{
+    labelList geometries(IndirectList<label>(surfaces_, surfacesToTest));
+
+    // Do the tests. Note that findNearest returns index in geometries.
+    List<pointIndexHit> hitInfo;
+    searchableSurfacesQueries::findNearest
+    (
+        allGeometry_,
+        geometries,
+        samples,
+        nearestDistSqr,
+        hitSurface,
+        hitInfo
+    );
+
+    // Rework the hitSurface to be surface (i.e. index into surfaces_)
+    forAll(hitSurface, pointI)
+    {
+        if (hitSurface[pointI] != -1)
+        {
+            hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
+        }
+    }
+
+    // Collect the region
+    hitRegion.setSize(hitSurface.size());
+    hitRegion = -1;
+
+    forAll(surfacesToTest, i)
+    {
+        label surfI = surfacesToTest[i];
+
+        // Collect hits for surfI
+        const labelList localIndices(findIndices(hitSurface, surfI));
+
+        List<pointIndexHit> localHits
+        (
+            IndirectList<pointIndexHit> 
+            (
+                hitInfo,
+                localIndices
+            )
+        );
+
+        labelList localRegion;
+        allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
+
+        forAll(localIndices, i)
+        {
+            hitRegion[localIndices[i]] = localRegion[i];
+        }
+    }
+}
+
+
 //// Find intersection with max of edge. Return -1 or the surface
 //// with the highest maxLevel above currentLevel
 //Foam::label Foam::refinementSurfaces::findHighestIntersection
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index edb613d966..ee5dd219a5 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -154,7 +154,10 @@ public:
                 return cellZoneNames_;
             }
 
-            //- Get indices of named surfaces (surfaces with cellZoneName)
+            //- Get indices of named surfaces (surfaces with faceZoneName)
+            labelList getUnnamedSurfaces() const;
+
+            //- Get indices of named surfaces (surfaces with faceZoneName)
             labelList getNamedSurfaces() const;
 
             //- Get indices of closed named surfaces
@@ -249,6 +252,8 @@ public:
 
             //- Find intersection nearest to the endpoints. surface1,2 are
             //  not indices into surfacesToTest but refinement surface indices.
+            //  Returns surface, region on surface (so not global surface)
+            //  and position on surface.
             void findNearestIntersection
             (
                 const labelList& surfacesToTest,
@@ -282,6 +287,17 @@ public:
                 List<pointIndexHit>&
             ) const;
 
+            //- Find nearest point on surfaces. Return surface and region on
+            //  surface (so not global surface)
+            void findNearestRegion
+            (
+                const labelList& surfacesToTest,
+                const pointField& samples,
+                const  scalarField& nearestDistSqr,
+                labelList& hitSurface,
+                labelList& hitRegion
+            ) const;
+
             //- Detect if a point is 'inside' (closed) surfaces.
             //  Returns -1 if not, returns first surface it is.
             void findInside
diff --git a/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C b/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
index c5927a6128..907eb097f1 100644
--- a/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
+++ b/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
@@ -301,7 +301,7 @@ void Foam::layerAdditionRemoval::removeCellLayer
 
         // Is any of the faces a boundary face?  If so, grab the patch
         // A boundary-to-boundary collapse is checked for in validCollapse()
-        // and cannnot happen here.  
+        // and cannot happen here.  
 
         if (!mesh.isInternalFace(mf[faceI]))
         {
diff --git a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C
index 1606e29b07..02471bdc2e 100644
--- a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C
+++ b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C
@@ -50,7 +50,8 @@ Foam::Analytical<Type>::~Analytical()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-Type Foam::Analytical<Type>::integrate
+typename Foam::IntegrationScheme<Type>::integrationResult
+Foam::Analytical<Type>::integrate
 (
     const Type phi,
     const scalar dt,
@@ -58,7 +59,11 @@ Type Foam::Analytical<Type>::integrate
     const scalar beta
 ) const
 {
-    return alpha + (phi - alpha)*exp(-beta*dt);
+    typename IntegrationScheme<Type>::integrationResult retValue;
+    retValue.average() = alpha + (phi - alpha)*(1 - exp(-beta*dt))/(beta*dt);
+    retValue.value() =  alpha + (phi - alpha)*exp(-beta*dt);
+
+    return retValue;
 }
 
 
diff --git a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H
index ebbb61522c..3c6bdd3b1f 100644
--- a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H
+++ b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H
@@ -74,7 +74,7 @@ public:
     // Member Functions
 
         //- Perform the integration
-        virtual Type integrate
+        virtual typename IntegrationScheme<Type>::integrationResult integrate
         (
             const Type phi,
             const scalar dt,
diff --git a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C
index 5b9f921819..2b39e9fd68 100644
--- a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C
+++ b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C
@@ -50,7 +50,8 @@ Foam::Euler<Type>::~Euler()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-Type Foam::Euler<Type>::integrate
+typename Foam::IntegrationScheme<Type>::integrationResult
+Foam::Euler<Type>::integrate
 (
     const Type phi,
     const scalar dt,
@@ -58,7 +59,11 @@ Type Foam::Euler<Type>::integrate
     const scalar beta
 ) const
 {
-    return (phi + dt*alpha)/(1.0 + dt/beta);
+    typename IntegrationScheme<Type>::integrationResult retValue;
+    retValue.value() = (phi + beta*dt*alpha)/(1.0 + beta*dt);
+    retValue.average() = 0.5*(phi + retValue.value());
+
+    return retValue;
 }
 
 
diff --git a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H
index 724320a09e..a3c925e9f0 100644
--- a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H
+++ b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H
@@ -74,7 +74,7 @@ public:
     // Member Functions
 
         //- Perform the integration
-        virtual Type integrate
+        virtual typename IntegrationScheme<Type>::integrationResult integrate
         (
             const Type phi,
             const scalar dt,
diff --git a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H
index b8b3859a0f..bac74a898c 100644
--- a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H
+++ b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H
@@ -53,6 +53,63 @@ namespace Foam
 template<class Type>
 class IntegrationScheme
 {
+
+public:
+
+    //- Helper class to supply results of integration
+    class integrationResult
+    {
+        //- Integration value
+        Type value_;
+
+        //- Average value across integration step
+        Type average_;
+
+
+    public:
+
+        //- Constructor
+        integrationResult()
+        :
+            value_(pTraits<Type>::zero),
+            average_(pTraits<Type>::zero)
+        {}
+
+
+        // Member functions
+
+            // Access
+
+                //- Return const access to the value
+                Type value() const
+                {
+                    return value_;
+                }
+
+                //- Return const access to the average
+                Type average() const
+                {
+                    return average_;
+                }
+
+
+            // Edit
+
+                //- Return access to the value for changing
+                Type& value()
+                {
+                    return value_;
+                }
+
+                //- Return access to the average for changing
+                Type& average()
+                {
+                    return average_;
+                }
+        };
+
+
+private:
     // Private data
 
         //- Name of the Integration variable
@@ -120,7 +177,7 @@ public:
     // Member Functions
 
         //- Perform the Integration
-        virtual Type integrate
+        virtual integrationResult integrate
         (
             const Type phi,
             const scalar dt,
@@ -151,7 +208,7 @@ public:
                                                                               \
     defineNamedTemplateTypeNameAndDebug(SS<Type>, 0);                         \
                                                                               \
-    IntegrationScheme<Type>::adddictionaryConstructorToTable<SS<Type> >          \
+    IntegrationScheme<Type>::adddictionaryConstructorToTable<SS<Type> >       \
         add##SS##Type##ConstructorToTable_;
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 05c5dc728f..6ee04f1235 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -332,21 +332,11 @@ void Foam::KinematicCloud<ParcelType>::evolve()
 
     inject(td);
 
-    move(td);
-}
-
-
-template<class ParcelType>
-template<class TrackingData>
-void Foam::KinematicCloud<ParcelType>::move
-(
-    TrackingData& td
-)
-{
     if (coupled_)
     {
         resetSourceTerms();
     }
+
     Cloud<ParcelType>::move(td);
 }
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index dff0859c53..d0925be0f5 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -257,10 +257,6 @@ protected:
             ParcelType* p
         );
 
-        //- Move the parcels
-        template<class TrackingData>
-        void move(TrackingData& td);
-
         //- Post-injection checks
         void postInjectCheck();
 
@@ -436,7 +432,7 @@ public:
             //- Reset the spray source terms
             void resetSourceTerms();
 
-            //- Evolve the spray (move, inject)
+            //- Evolve the spray (inject, inject)
             void evolve();
 };
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index f30f70bcb3..dcf4730f7a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -185,7 +185,12 @@ void Foam::ReactingCloud<ParcelType>::evolve()
 
     inject(td);
 
-    this->move(td);
+    if (this->coupled())
+    {
+        resetSourceTerms();
+    }
+
+    Cloud<ParcelType>::move(td);
 }
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
index 7afa096c63..6f7bb05e09 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
@@ -204,7 +204,7 @@ public:
             //- Reset the spray source terms
             void resetSourceTerms();
 
-            //- Evolve the spray (move, inject)
+            //- Evolve the spray (inject, move)
             void evolve();
 };
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index 0e07f89769..ad0d520b74 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -176,7 +176,12 @@ void Foam::ThermoCloud<ParcelType>::evolve()
 
     inject(td);
 
-    this->move(td);
+    if (this->coupled())
+    {
+        resetSourceTerms();
+    }
+
+    Cloud<ParcelType>::move(td);
 }
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
index da9c5e4fb5..072c86165b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
@@ -209,11 +209,11 @@ public:
 
         // Cloud evolution functions
 
-            //- Evolve the spray (move, inject)
-            void evolve();
-
             //- Reset the spray source terms
             void resetSourceTerms();
+
+            //- Evolve the spray (inject, move)
+            void evolve();
 };
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 9b0dbf45cd..07e20192c8 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -155,15 +155,16 @@ Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Update velocity - treat as 3-D
-    const scalar bp = 1.0/(Cud + VSMALL);
-    const vector ap = Uc_/bp + rhoc_/rho_*td.g();
+    const vector ap = Uc_ + (1 - rhoc_/rho_)/(Cud + VSMALL)*td.g();
+    const scalar bp = Cud;
 
-    vector Unew = td.cloud().UIntegrator().integrate(U_, dt, ap, bp);
+    vector Unew = td.cloud().UIntegrator().integrate(U_, dt, ap, bp).value();
 
 //    Info<< "U_, Unew = " << U_ << ", " << Unew << endl;
 
     // Calculate the momentum transfer to the continuous phase
-    dUTrans = -mass()*(Unew - U_);
+    // - do not include gravity impulse
+    dUTrans = -mass()*(Unew - U_ - dt*td.g());
 
     // Make corrections for 2-D cases
     if (meshInfo.caseIs2d())
@@ -233,13 +234,17 @@ bool Foam::KinematicParcel<ParcelType>::move
         // Update cell based properties
         p.updateCellQuantities(td, dt, celli);
 
-        if (td.cloud().coupled())
+        // Avoid problems with extremely small timesteps
+        if (dt > ROOTVSMALL)
         {
-            p.calcCoupled(td, dt, celli);
-        }
-        else
-        {
-            p.calcUncoupled(td, dt, celli);
+            if (td.cloud().coupled())
+            {
+                p.calcCoupled(td, dt, celli);
+            }
+            else
+            {
+                p.calcUncoupled(td, dt, celli);
+            }
         }
 
         if (p.onBoundary() && td.keepParticle)
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
index 98ecd22675..b000ed0237 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
@@ -203,7 +203,7 @@ Foam::Ostream& Foam::operator<<
         os  << static_cast<const Particle<ParcelType>& >(p);
         os.write
         (
-            reinterpret_cast<const char*>(p.typeId()),
+            reinterpret_cast<const char*>(&p.typeId_),
             sizeof(p.typeId())
           + sizeof(p.d())
           + sizeof(p.U())
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index d9d47894a5..7fe198e773 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -59,7 +59,6 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     const scalar mass0 = this->mass();
     const scalar np0 = this->nParticle_;
     const scalar T0 = this->T_;
-    const scalar cp0 = this->cp_;
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
     // Initialise transfer terms
@@ -79,9 +78,12 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // - components do not necessarily exist in particle already
     scalarList dMassSR(td.cloud().gases().size(), 0.0);
 
-    // Total mass lost from particle due to surface reactions
-    // - sub-model will adjust component mass fractions
-    scalar dMassMTSR = 0.0;
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Calculate velocity - update U
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    scalar Cud = 0.0;
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -91,13 +93,6 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate velocity - update U
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
-
-
     // ~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate mass transfer
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -107,27 +102,36 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate surface reactions
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMTSR, dMassSR);
+
+    // Initialise enthalpy retention to zero
+    scalar dhRet = 0.0;
+
+    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
 
     // New total mass
-    const scalar mass1 = mass0 - sum(dMassMT) - dMassMTSR;
+    const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
 
-    // Ratio of mass devolatilised to the total volatile mass of the particle
-    const scalar fVol = 1 -
-        (YMixture_[0]*mass1)
-       /(td.cloud().composition().YMixture0()[0]*mass0_);
+    // Correct particle temperature to account for latent heat of
+    // devolatilisation
+    T1 -=
+        td.constProps().Ldevol()
+       *sum(dMassMT)
+       /(0.5*(mass0 + mass1)*this->cp_);
 
-    // Specific heat capacity of non-volatile components
-    const scalar cpNonVolatile =
-        (
-            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_)
-          + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_)
-        )/(YMixture_[1] + YMixture_[2]);
+    // Add retained enthalpy from surface reaction to particle and remove
+    // from gas
+    T1 += dhRet/(0.5*(mass0 + mass1)*this->cp_);
+    dhTrans -= dhRet;
 
-    // New specific heat capacity - linear variation until volatiles
-    // have evolved
-    const scalar cp1 = (cpNonVolatile - td.constProps().cp0())*fVol
-       + td.constProps().cp0();
+    // Correct dhTrans to account for enthalpy of evolved volatiles
+    dhTrans +=
+        sum(dMassMT)
+       *td.cloud().composition().HGas(YGas_, 0.5*(T0 + T1));
+
+    // Correct dhTrans to account for enthalpy of consumed solids
+    dhTrans +=
+        sum(dMassSR)
+       *td.cloud().composition().HSolid(YSolid_, 0.5*(T0 + T1));
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -147,8 +151,8 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     td.cloud().UCoeff()[celli] += np0*mass0*Cud;
 
     // Update enthalpy transfer
-//    td.cloud().hTrans()[celli] += np0*(mass0*cp0*T0 - mass1*cp1*T1);
-    td.cloud().hTrans()[celli] += np0*((mass0*cp0 - mass1*cp1)*T0 + dhTrans);
+    // - enthalpy of lost solids already accounted for
+    td.cloud().hTrans()[celli] += np0*dhTrans;
 
     // Accumulate coefficient to be applied in carrier phase enthalpy coupling
     td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
@@ -166,7 +170,12 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
         {
             td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i];
         }
-        td.cloud().hTrans()[celli] += np0*mass1*cp1*T1;
+        td.cloud().hTrans()[celli] +=
+            np0*mass1
+           *(
+                YMixture_[0]*td.cloud().composition().HGas(YGas_, T1)
+              + YMixture_[2]*td.cloud().composition().HSolid(YSolid_, T1)
+            );
         td.cloud().UTrans()[celli] += np0*mass1*U1;
     }
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -176,7 +185,10 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     {
         this->U_ = U1;
         this->T_ = T1;
-        this->cp_ = cp1;
+        this->cp_ =
+            YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1)
+          + YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1)
+          + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_);
 
         // Update particle density or diameter
         if (td.cloud().massTransfer().changesVolume())
@@ -205,7 +217,7 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     const scalar T0 = this->T_;
     const scalar mass0 = this->mass();
-//    const scalar cp0 = this->cp();
+    const scalar cp0 = this->cp_;
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
     // Initialise transfer terms
@@ -225,9 +237,12 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // - components do not necessarily exist in particle already
     scalarList dMassSR(td.cloud().gases().size(), 0.0);
 
-    // Total mass lost from particle due to surface reactions
-    // - sub-model will adjust component mass fractions
-    scalar dMassMTSR = 0.0;
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Calculate velocity - update U
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    scalar Cud = 0.0;
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -236,16 +251,6 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     scalar htc = 0.0;
     scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
 
-    // Limit new temp max by vapourisarion temperature
-    T1 = min(td.constProps().Tvap(), T1);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate velocity - update U
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
-
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate mass transfer
@@ -256,37 +261,23 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate surface reactions
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    calcSurfaceReactions
-    (
-        td,
-        dt,
-        celli,
-        T0,
-        T1,
-        dMassMTSR,
-        dMassSR
-    );
+
+    // Initialise enthalpy retention to zero
+    scalar dhRet = 0.0;
+
+    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
 
     // New total mass
-    const scalar mass1 = mass0 - sum(dMassMT) - dMassMTSR;
+    const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
 
-    // Ratio of mass devolatilised to the total volatile mass of the particle
-    const scalar fVol = 1 -
-        (YMixture_[0]*mass1)
-       /(td.cloud().composition().YMixture0()[0]*mass0_);
-
-    // Specific heat capacity of non-volatile components
-    const scalar cpNonVolatile =
-        (
-            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_)
-          + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_)
-        )/(YMixture_[1] + YMixture_[2]);
-
-    // New specific heat capacity - linear variation until volatiles
-    // have evolved
-    const scalar cp1 = (cpNonVolatile - td.constProps().cp0())*fVol
-       + td.constProps().cp0();
+    // New specific heat capacity
+    const scalar cp1 =
+        YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1)
+      + YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1)
+      + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_);
 
+    // Add retained enthalpy to particle
+    T1 += dhRet/(mass0*0.5*(cp0 + cp1));
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Remove the particle when mass falls below minimum threshold
@@ -389,8 +380,9 @@ void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
     const label celli,
     const scalar T0,
     const scalar T1,
-    scalar& dMassMTSR,
-    scalarList& dMassMT
+    const scalarList& dMassMT,
+    scalarList& dMassSR,
+    scalar& dhRet
 )
 {
     // Check that model is active
@@ -409,21 +401,20 @@ void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
         T0,
         T1,
         this->Tc_,
+        pc_,
         this->rhoc_,
         this->mass(),
+        dMassMT,
         YGas_,
         YLiquid_,
         YSolid_,
         YMixture_,
-        dMassMTSR,
-        dMassMT
+        dMassSR,
+        dhRet
     );
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 
 #include "ReactingParcelIO.C"
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index 43e125f13b..a6b8eb5ad5 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -89,6 +89,9 @@ public:
             //- Boiling point [K]
             const scalar Tbp_;
 
+            //- Latent heat of devolatilisation [J/kg]
+            const scalar Ldevol_;
+
 
     public:
 
@@ -102,6 +105,9 @@ public:
 
             //- Return const access to the boiling point
             inline scalar Tbp() const;
+
+            //- Return const access to the latent heat of devolatilisation
+            inline scalar Ldevol() const;
     };
 
 
@@ -210,8 +216,9 @@ protected:
             const label celli,
             const scalar T0,
             const scalar T1,
-            scalar& dMassMTSR,
-            scalarList& dMassMT
+            const scalarList& dMassMT,
+            scalarList& dMassSR,
+            scalar& dhRet
         );
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
index 7138666b2e..0c50b0a853 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
@@ -34,7 +34,8 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
 :
     ThermoParcel<ParcelType>::constantProperties(dict),
     Tvap_(dimensionedScalar(dict.lookup("Tvap")).value()),
-    Tbp_(dimensionedScalar(dict.lookup("Tbp")).value())
+    Tbp_(dimensionedScalar(dict.lookup("Tbp")).value()),
+    Ldevol_(dimensionedScalar(dict.lookup("Ldevol")).value())
 {}
 
 
@@ -127,6 +128,14 @@ Foam::ReactingParcel<ParcelType>::constantProperties::Tbp() const
 }
 
 
+template<class ParcelType>
+inline Foam::scalar
+Foam::ReactingParcel<ParcelType>::constantProperties::Ldevol() const
+{
+    return Ldevol_;
+}
+
+
 // * * * * * * * * * * * trackData Member Functions  * * * * * * * * * * * * //
 
 template<class ParcelType>
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
index ae0dcd0573..b1797b97de 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
@@ -297,10 +297,7 @@ Foam::Ostream& Foam::operator<<
         os  << static_cast<const ThermoParcel<ParcelType>& >(p);
         os.write
         (
-            reinterpret_cast<const char*>
-            (
-                &const_cast<ReactingParcel<ParcelType>&>(p).mass0()
-            ),
+            reinterpret_cast<const char*>(&p.mass0_),
             sizeof(p.mass0())
         );
         os << p.YMixture() << YGasLoc << YLiquidLoc << YSolidLoc;
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 89c0e4f2e4..75ab727b77 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -203,17 +203,17 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     // Set new particle temperature
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    scalar Tnew = td.cloud().TIntegrator().integrate(T_, dt, ap, bp);
+    // Integrate to find the new parcel temperature
+    IntegrationScheme<scalar>::integrationResult Tres =
+        td.cloud().TIntegrator().integrate(T_, dt, ap, bp);
 
-    dhTrans = -this->mass()*cp_*(Tnew - T_);
+    // Using average parcel temperature for enthalpy transfer calculation
+    dhTrans = dt*this->areaS()*htc*(Tres.average() - Tc_);
 
-    return Tnew;
+    return Tres.value();
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 
 #include "ThermoParcelIO.C"
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
index 04d444e6bb..d0242e5e51 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
@@ -149,10 +149,7 @@ Foam::Ostream& Foam::operator<<
         os  << static_cast<const KinematicParcel<ParcelType>& >(p);
         os.write
         (
-            reinterpret_cast<const char*>
-            (
-                &const_cast<ThermoParcel<ParcelType>&>(p).T()
-            ),
+            reinterpret_cast<const char*>(&p.T_),
             sizeof(p.T()) + sizeof(p.cp())
         );
     }
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
index ce85b19a5c..9caea94a34 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
@@ -209,13 +209,20 @@ public:
             //- Return the gas constant for the gas mixture
             virtual scalar RGas(const scalarField& YGas) const = 0;
 
-            //- Return enthalpy for the gas mixture
+            //- Return enthalpy for the gas mixture [energy per unit mass]
             virtual scalar HGas
             (
                  const scalarField& YGas,
                  const scalar T
             ) const = 0;
 
+            //- Return enthalpy for the solid mixture [energy per unit mass]
+            virtual scalar HSolid
+            (
+                const scalarField& YSolid,
+                const scalar T
+            ) const = 0;
+
             //- Return specific heat caparcity for the gas mixture
             virtual scalar cpGas
             (
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.C b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.C
index 350546b5fe..74fb92338c 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.C
@@ -519,6 +519,29 @@ const
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::SingleMixtureFraction<CloudType>::HSolid
+(
+    const scalarField& YSolid,
+    const scalar T
+)
+const
+{
+    scalar HMixture = 0.0;
+    forAll(YSolid, i)
+    {
+        label id = solidGlobalIds_[i];
+        HMixture +=
+            YSolid[i]
+           *(
+                this->solids().properties()[id].Hf()
+              + this->solids().properties()[id].cp()*T
+            );
+    }
+    return HMixture;
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::SingleMixtureFraction<CloudType>::cpGas
 (
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.H b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.H
index bd03d69456..de1fb8daf7 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.H
@@ -188,9 +188,12 @@ public:
             //- Return the gas constant for the gas mixture
             scalar RGas(const scalarField& YGas) const;
 
-            //- Return enthalpy for the gas mixture
+            //- Return enthalpy for the gas mixture [energy per unit mass]
             scalar HGas(const scalarField& YGas, const scalar T) const;
 
+            //- Return enthalpy for the solid mixture [energy per unit mass]
+            scalar HSolid(const scalarField& YSolid, const scalar T) const;
+
             //- Return specific heat caparcity for the gas mixture
             scalar cpGas(const scalarField& YGas, const scalar T) const;
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
index 8ea8816f04..d164a500c2 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
@@ -64,16 +64,20 @@ void Foam::NoSurfaceReaction<CloudType>::calculate
     const scalar T0,
     const scalar T1,
     const scalar Tc,
+    const scalar pc,
     const scalar rhoc,
     const scalar massp,
+    const scalarList& dMassMT,
     scalarField& YGas,
     scalarField& YLiquid,
     scalarField& YSolid,
     scalarField& YMixture,
-    scalar& dMassMTSR,
-    scalarList& dMassSR
+    scalarList& dMassSR,
+    scalar& dhRet
 ) const
-{}
+{
+    // do nothing
+}
 
 
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
index f3072df598..8a58352b5d 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
@@ -84,14 +84,16 @@ public:
             const scalar T0,
             const scalar T1,
             const scalar Tc,
+            const scalar pc,
             const scalar rhoc,
             const scalar massp,
+            const scalarList& dMassMT,
             scalarField& YGas,
             scalarField& YLiquid,
             scalarField& YSolid,
             scalarField& YMixture,
-            scalar& dMassMTSR,
-            scalarList& dMassSR
+            scalarList& dMassSR,
+            scalar& dhRet
         ) const;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
index b50b5a90b4..a3b3e4e00f 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
@@ -141,14 +141,16 @@ public:
             const scalar T0,
             const scalar T1,
             const scalar Tc,
+            const scalar pc,
             const scalar rhoc,
             const scalar massp,
+            const scalarList& dMassMT,
             scalarField& YGas,
             scalarField& YLiquid,
             scalarField& YSolid,
             scalarField& YMixture,
-            scalar& dMassMTSR,
-            scalarList& dMassSR
+            scalarList& dMassSR,
+            scalar& dhRet
         ) const = 0;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/NoHeatTransfer/NoHeatTransfer.C b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/NoHeatTransfer/NoHeatTransfer.C
index edeb799fc0..92ddbd7663 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/NoHeatTransfer/NoHeatTransfer.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/NoHeatTransfer/NoHeatTransfer.C
@@ -62,11 +62,6 @@ Foam::scalar Foam::NoHeatTransfer<CloudType>::Nu
     const scalar
 ) const
 {
-    notImplemented
-    (
-        "Foam::scalar Foam::NoHeatTransfer<CloudType>::Nu"
-        "(const scalar, const scalar)"
-    );
     return 0.0;
 }
 
@@ -74,8 +69,7 @@ Foam::scalar Foam::NoHeatTransfer<CloudType>::Nu
 template <class CloudType>
 Foam::scalar Foam::NoHeatTransfer<CloudType>::Pr() const
 {
-    notImplemented("Foam::scalar Foam::NoHeatTransfer<CloudType>::Pr()");
-    return 0.0;
+    return 1.0;
 }
 
 
diff --git a/src/meshTools/searchableSurface/searchableBox.C b/src/meshTools/searchableSurface/searchableBox.C
index b1776b3c9e..1858905d49 100644
--- a/src/meshTools/searchableSurface/searchableBox.C
+++ b/src/meshTools/searchableSurface/searchableBox.C
@@ -436,7 +436,6 @@ void Foam::searchableBox::findLineAll
     // Work array
     DynamicList<pointIndexHit, 1, 1> hits;
 
-//XXX
     // Tolerances:
     // To find all intersections we add a small vector to the last intersection
     // This is chosen such that
diff --git a/tutorials/pimpleFoam/t-junction/0/U b/tutorials/pimpleFoam/t-junction/0/U
new file mode 100644
index 0000000000..8851e66583
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/0/U
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet      
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+    }
+
+    outlet1    
+    {
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
+    }
+
+    outlet2   
+    {
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
+    }
+
+    defaultFaces      
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/0/epsilon b/tutorials/pimpleFoam/t-junction/0/epsilon
new file mode 100644
index 0000000000..8e5b3e553d
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/0/epsilon
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    inlet
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.01;         // 1cm - half channel height
+        value           uniform 1;
+    }
+
+    outlet1
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+    }
+
+    outlet2
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+    }
+
+    defaultFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/0/k b/tutorials/pimpleFoam/t-junction/0/k
new file mode 100644
index 0000000000..35b2bd3ff9
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/0/k
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    inlet
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.05;           // 5% turbulence
+        value           uniform 1;
+    }
+
+    outlet1
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+    }
+
+    outlet2
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+    }
+
+    defaultFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/0/nuTilda b/tutorials/pimpleFoam/t-junction/0/nuTilda
new file mode 100644
index 0000000000..6ea76a6b63
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/0/nuTilda
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      nuTilda;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet1
+    {
+        type            zeroGradient;
+    }
+
+    outlet2
+    {
+        type            zeroGradient;
+    }
+
+    defaultFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/0/p b/tutorials/pimpleFoam/t-junction/0/p
new file mode 100644
index 0000000000..b3d4646058
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/0/p
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    inlet      
+    {
+        type            totalPressure;
+        p0              uniform 100040;
+        U               U;
+        phi             phi;
+        rho             none;
+        psi             none;
+        gamma           1;
+        value           uniform 100040;
+    }
+
+    outlet1     
+    {
+        type            fixedValue;
+        value           uniform 100010;
+    }
+
+    outlet2    
+    {
+        type            fixedValue;
+        value           uniform 100000;
+    }
+
+    defaultFaces
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/README.txt b/tutorials/pimpleFoam/t-junction/README.txt
new file mode 100644
index 0000000000..0905f8784d
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/README.txt
@@ -0,0 +1,2 @@
+15/8/8 Simple T-junction. Inlet on left, one outlet at bottom, one at top.
+To test multiple outlets.
diff --git a/tutorials/pimpleFoam/t-junction/constant/RASProperties b/tutorials/pimpleFoam/t-junction/constant/RASProperties
new file mode 100644
index 0000000000..e46e0a858b
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/constant/RASProperties
@@ -0,0 +1,200 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel            kEpsilon;
+
+turbulence          on;
+
+printCoeffs         on;
+
+laminarCoeffs
+{
+}
+
+kEpsilonCoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphaEps         0.76923;
+}
+
+RNGkEpsilonCoeffs
+{
+    Cmu              0.0845;
+    C1               1.42;
+    C2               1.68;
+    alphak           1.39;
+    alphaEps         1.39;
+    eta0             4.38;
+    beta             0.012;
+}
+
+realizableKECoeffs
+{
+    Cmu              0.09;
+    A0               4.0;
+    C2               1.9;
+    alphak           1;
+    alphaEps         0.833333;
+}
+
+kOmegaSSTCoeffs
+{
+    alphaK1          0.85034;
+    alphaK2          1.0;
+    alphaOmega1      0.5;
+    alphaOmega2      0.85616;
+    gamma1           0.5532;
+    gamma2           0.4403;
+    beta1            0.0750;
+    beta2            0.0828;
+    betaStar         0.09;
+    a1               0.31;
+    c1               10;
+
+    Cmu              0.09;
+}
+
+NonlinearKEShihCoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphak           1;
+    alphaEps         0.76932;
+    A1               1.25;
+    A2               1000;
+    Ctau1            -4;
+    Ctau2            13;
+    Ctau3            -2;
+    alphaKsi         0.9;
+}
+
+LienCubicKECoeffs
+{
+    C1               1.44;
+    C2               1.92;
+    alphak           1;
+    alphaEps         0.76923;
+    A1               1.25;
+    A2               1000;
+    Ctau1            -4;
+    Ctau2            13;
+    Ctau3            -2;
+    alphaKsi         0.9;
+}
+
+QZetaCoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphaZeta        0.76923;
+    anisotropic     no;
+}
+
+LaunderSharmaKECoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphaEps         0.76923;
+}
+
+LamBremhorstKECoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphaEps         0.76923;
+}
+
+LienCubicKELowReCoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphak           1;
+    alphaEps         0.76923;
+    A1               1.25;
+    A2               1000;
+    Ctau1            -4;
+    Ctau2            13;
+    Ctau3            -2;
+    alphaKsi         0.9;
+    Am               0.016;
+    Aepsilon         0.263;
+    Amu              0.00222;
+}
+
+LienLeschzinerLowReCoeffs
+{
+    Cmu              0.09;
+    C1               1.44;
+    C2               1.92;
+    alphak           1;
+    alphaEps         0.76923;
+    Am               0.016;
+    Aepsilon         0.263;
+    Amu              0.00222;
+}
+
+LRRCoeffs
+{
+    Cmu              0.09;
+    Clrr1            1.8;
+    Clrr2            0.6;
+    C1               1.44;
+    C2               1.92;
+    Cs               0.25;
+    Ceps             0.15;
+    alphaEps         0.76923;
+}
+
+LaunderGibsonRSTMCoeffs
+{
+    Cmu              0.09;
+    Clg1             1.8;
+    Clg2             0.6;
+    C1               1.44;
+    C2               1.92;
+    C1Ref            0.5;
+    C2Ref            0.3;
+    Cs               0.25;
+    Ceps             0.15;
+    alphaEps         0.76923;
+    alphaR           1.22;
+}
+
+SpalartAllmarasCoeffs
+{
+    alphaNut         1.5;
+    Cb1              0.1355;
+    Cb2              0.622;
+    Cw2              0.3;
+    Cw3              2;
+    Cv1              7.1;
+    Cv2              5.0;
+}
+
+wallFunctionCoeffs
+{
+    kappa            0.4187;
+    E                9;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/constant/polyMesh/blockMeshDict b/tutorials/pimpleFoam/t-junction/constant/polyMesh/blockMeshDict
new file mode 100644
index 0000000000..d8f6e3bb6e
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/constant/polyMesh/blockMeshDict
@@ -0,0 +1,109 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//           outlet1
+//             +-+
+//             | |
+//             | |
+//             | |
+//             | |
+// +-----------+ |
+// |inlet        |
+// +-----------+ |
+//             | |
+//             | |
+//             | |
+//             | |
+//             +-+
+//           outlet2
+
+convertToMeters 1;
+
+vertices        
+(
+    (0.0  -0.01 0)   //0
+    (0.2  -0.01 0)
+    (0.2   0.01 0)   //2
+    (0.0   0.01 0)
+
+    (0.22 -0.01 0)  //4
+    (0.22  0.01 0)
+
+    (0.2  -0.21 0)  //6
+    (0.22 -0.21 0)
+
+    (0.2   0.21 0)  //8
+    (0.22  0.21 0)
+
+    // Z
+    (0.0  -0.01 0.02)   //0
+    (0.2  -0.01 0.02)
+    (0.2   0.01 0.02)   //2
+    (0.0   0.01 0.02)
+
+    (0.22 -0.01 0.02)  //4
+    (0.22  0.01 0.02)
+
+    (0.2  -0.21 0.02)  //6
+    (0.22 -0.21 0.02)
+
+    (0.2   0.21 0.02)  //8
+    (0.22  0.21 0.02)
+
+);
+
+blocks          
+(
+    // inlet block
+    hex (0 1 2 3  10 11 12 13) (50 5 5) simpleGrading (1 1 1)
+
+    // central block
+    hex (1 4 5 2  11 14 15 12) (5 5 5) simpleGrading (1 1 1)
+
+    // bottom block
+    hex (6 7 4 1  16 17 14 11) (5 50 5) simpleGrading (1 1 1)
+
+    // top block
+    hex (2 5 9 8  12 15 19 18) (5 50 5) simpleGrading (1 1 1)
+);
+
+edges           
+(
+);
+
+patches         
+(
+    patch inlet
+    (
+        (0 10 13 3)
+    )
+
+    patch outlet1
+    (
+        (6 7 17 16)
+    )
+
+    patch outlet2
+    (
+        (8 18 19 9)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/constant/polyMesh/boundary b/tutorials/pimpleFoam/t-junction/constant/polyMesh/boundary
new file mode 100644
index 0000000000..4ea7f95a3e
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/constant/polyMesh/boundary
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+4
+(
+    inlet
+    {
+        type            patch;
+        nFaces          25;
+        startFace       10050;
+    }
+    outlet1
+    {
+        type            patch;
+        nFaces          25;
+        startFace       10075;
+    }
+    outlet2
+    {
+        type            patch;
+        nFaces          25;
+        startFace       10100;
+    }
+    defaultFaces
+    {
+        type            wall;
+        nFaces          3075;
+        startFace       10125;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/constant/transportProperties b/tutorials/pimpleFoam/t-junction/constant/transportProperties
new file mode 100644
index 0000000000..fd5fc94660
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/constant/transportProperties
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              nu [0 2 -1 0 0 0 0] 1e-05;
+
+CrossPowerLawCoeffs
+{
+    nu0             nu0 [0 2 -1 0 0 0 0] 1e-06;
+    nuInf           nuInf [0 2 -1 0 0 0 0] 1e-06;
+    m               m [0 0 1 0 0 0 0] 1;
+    n               n [0 0 0 0 0 0 0] 1;
+}
+
+BirdCarreauCoeffs
+{
+    nu0             nu0 [0 2 -1 0 0 0 0] 1e-06;
+    nuInf           nuInf [0 2 -1 0 0 0 0] 1e-06;
+    k               k [0 0 1 0 0 0 0] 0;
+    n               n [0 0 0 0 0 0 0] 1;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/system/controlDict b/tutorials/pimpleFoam/t-junction/system/controlDict
new file mode 100644
index 0000000000..b5f0d20f7e
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/system/controlDict
@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application turbFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+maxCo           5;
+
+functions
+(
+    probes
+    {
+        // Type of functionObject
+        type probes;
+
+        // Where to load it from (if not already in solver)
+        functionObjectLibs ("libsampling.so");
+
+        // Name of the directory for the probe data
+        name probes;
+
+        // Locations to be probed. runTime modifiable!
+        probeLocations
+        (
+            (1E-6 0 0.01)           // at inlet
+            (0.21 -0.20999 0.01)    // at outlet1
+            (0.21  0.20999 0.01)    // at outlet2
+            (0.21 0 0.01)           // at central block
+        );
+
+        // Fields to be probed. runTime modifiable!
+        fields
+        (
+            p
+            U
+        );
+    }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/system/fvSchemes b/tutorials/pimpleFoam/t-junction/system/fvSchemes
new file mode 100644
index 0000000000..bbe7fab0d2
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/system/fvSchemes
@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+    grad(U)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss limitedLinearV 1;
+    div(phi,k)      Gauss limitedLinear 1;
+    div(phi,epsilon) Gauss limitedLinear 1;
+    div(phi,R)      Gauss limitedLinear 1;
+    div(R)          Gauss linear;
+    div(phi,nuTilda) Gauss limitedLinear 1;
+    div((nuEff*dev(grad(U).T()))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(nuEff,U) Gauss linear corrected;
+    laplacian((1|A(U)),p) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DREff,R) Gauss linear corrected;
+    laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+    interpolate(U)  linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/pimpleFoam/t-junction/system/fvSolution b/tutorials/pimpleFoam/t-junction/system/fvSolution
new file mode 100644
index 0000000000..995edf623c
--- /dev/null
+++ b/tutorials/pimpleFoam/t-junction/system/fvSolution
@@ -0,0 +1,88 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p GAMG
+    {
+        tolerance        1e-6;
+        relTol           0.01;
+
+        smoother         GaussSeidel;
+
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator     faceAreaPair;
+        mergeLevels      1;
+    };
+    pFinal GAMG
+    {
+        tolerance        1e-6;
+        relTol           0.0;
+
+        smoother         GaussSeidel;
+
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator     faceAreaPair;
+        mergeLevels      1;
+    };
+
+    U PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0.1;
+    };
+    UFinal PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0;
+    };
+    k PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0;
+    };
+    epsilon PBiCG
+    {
+        preconditioner   DILU;
+        tolerance        1e-05;
+        relTol           0;
+    };
+}
+
+PIMPLE
+{
+    nOuterCorrectors 2;
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       0;
+}
+
+relaxationFactors
+{
+    //p               0.3;
+    U               1.0;
+    k               1.0;
+    epsilon         1.0;
+}
+
+
+// ************************************************************************* //