diff --git a/applications/solvers/basic/potentialFoam/createFields.H b/applications/solvers/basic/potentialFoam/createFields.H
index faad23831295bde396f81282e8f7efb89534cb54..af0adfd8f9f7ab2099d744f79d93ef1ac6a82b12 100644
--- a/applications/solvers/basic/potentialFoam/createFields.H
+++ b/applications/solvers/basic/potentialFoam/createFields.H
@@ -93,21 +93,48 @@ if (args.optionFound("writep"))
 
 
 Info<< "Constructing velocity potential field Phi\n" << endl;
-volScalarField Phi
+autoPtr<volScalarField> PhiPtr;
+
+IOobject io
 (
-    IOobject
-    (
-        "Phi",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::NO_WRITE
-    ),
+    "Phi",
+    runTime.timeName(),
     mesh,
-    dimensionedScalar("Phi", dimLength*dimVelocity, 0),
-    p.boundaryField().types()
+    IOobject::MUST_READ,
+    IOobject::NO_WRITE
 );
 
+if (io.typeHeaderOk<volScalarField>())
+{
+    PhiPtr.reset(new volScalarField(io, mesh));
+}
+else
+{
+    // Cannot just use p.boundaryField().types() since does not initialise
+    // complex boundary types. Instead re-clone them from p.
+
+    io.readOpt() = IOobject::NO_READ;
+    PhiPtr.reset
+    (
+        new volScalarField
+        (
+            io,
+            mesh,
+            dimensionedScalar("Phi", dimLength*dimVelocity, 0),
+            p.boundaryField().types()
+        )
+    );
+
+    const volScalarField::GeometricBoundaryField& bp = p.boundaryField();
+    volScalarField::GeometricBoundaryField& bPhi = PhiPtr().boundaryField();
+
+    forAll(bp, patchI)
+    {
+        bPhi.set(patchI, bp[patchI].clone(PhiPtr().dimensionedInternalField()));
+    }
+}
+volScalarField& Phi = PhiPtr();
+
 label PhiRefCell = 0;
 scalar PhiRefValue = 0;
 setRefCell
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriterTemplates.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriterTemplates.C
index da660b779ca2b4c53d21472789eef70e1a7e7891..01518d3444148f8513c0bf35e2937276c3a66963 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriterTemplates.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriterTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,7 +49,8 @@ void Foam::lagrangianWriter::writeIOField(const wordList& objects)
 
         IOField<Type> fld(header);
 
-        os_ << object << ' ' << pTraits<Type>::nComponents << ' '
+        os_ << object << ' '
+            << int(pTraits<Type>::nComponents) << ' '
             << fld.size() << " float" << std::endl;
 
         DynamicList<floatScalar> fField(pTraits<Type>::nComponents*fld.size());
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/patchWriterTemplates.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/patchWriterTemplates.C
index f6ea96d7b8bfe279efc6451e5e5e8986e6740df4..35915339ccf3ed57048ddcce63188e034ffac2bb 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/patchWriterTemplates.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/patchWriterTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,7 +38,8 @@ void Foam::patchWriter::write
     {
         const GeometricField<Type, fvPatchField, volMesh>& fld = flds[fieldI];
 
-        os_ << fld.name() << ' ' << pTraits<Type>::nComponents << ' '
+        os_ << fld.name() << ' '
+            << int(pTraits<Type>::nComponents) << ' '
             << nFaces_ << " float" << std::endl;
 
         DynamicList<floatScalar> fField(pTraits<Type>::nComponents*nFaces_);
@@ -74,7 +75,8 @@ void Foam::patchWriter::write
         const GeometricField<Type, pointPatchField, pointMesh>& fld =
             flds[fieldI];
 
-        os_ << fld.name() << ' ' << pTraits<Type>::nComponents << ' '
+        os_ << fld.name() << ' '
+            << int(pTraits<Type>::nComponents) << ' '
             << nPoints_ << " float" << std::endl;
 
         DynamicList<floatScalar> fField(pTraits<Type>::nComponents*nPoints_);
@@ -103,7 +105,8 @@ void Foam::patchWriter::write
     {
         const GeometricField<Type, fvPatchField, volMesh>& fld = flds[fieldI];
 
-        os_ << fld.name() << ' ' << pTraits<Type>::nComponents << ' '
+        os_ << fld.name() << ' '
+            << int(pTraits<Type>::nComponents) << ' '
             << nPoints_ << " float" << std::endl;
 
         DynamicList<floatScalar> fField(pTraits<Type>::nComponents*nPoints_);
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/surfaceMeshWriterTemplates.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/surfaceMeshWriterTemplates.C
index c3412a1faa4d6a6738f803b07398abf546379a5d..dc41ae03c364df1bb44f1361adba925c723d8975 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/surfaceMeshWriterTemplates.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/surfaceMeshWriterTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -71,7 +71,8 @@ void Foam::surfaceMeshWriter::write
         const GeometricField<Type, fvsPatchField, surfaceMesh>& fld =
             sflds[fieldI];
 
-        os_ << fld.name() << ' ' << pTraits<Type>::nComponents << ' '
+        os_ << fld.name() << ' '
+            << int(pTraits<Type>::nComponents) << ' '
             << pp_.size() << " float" << std::endl;
 
         DynamicList<floatScalar> fField(pTraits<Type>::nComponents*pp_.size());
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFunsTemplates.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFunsTemplates.C
index 6f0c368bec2e99de675fb21d8bd87c41b6a472bf..842ae84e4f324d42b450edb06c54dd224bdc0fe8 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFunsTemplates.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFunsTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -74,7 +74,8 @@ void Foam::writeFuns::write
 
     label nValues = mesh.nCells() + superCells.size();
 
-    os  << vvf.name() << ' ' << pTraits<Type>::nComponents << ' '
+    os  << vvf.name() << ' '
+        << int(pTraits<Type>::nComponents) << ' '
         << nValues << " float" << std::endl;
 
     DynamicList<floatScalar> fField(pTraits<Type>::nComponents*nValues);
@@ -106,7 +107,8 @@ void Foam::writeFuns::write
     const labelList& addPointCellLabels = topo.addPointCellLabels();
     const label nTotPoints = mesh.nPoints() + addPointCellLabels.size();
 
-    os  << pvf.name() << ' ' << pTraits<Type>::nComponents << ' '
+    os  << pvf.name() << ' '
+        << int(pTraits<Type>::nComponents) << ' '
         << nTotPoints << " float" << std::endl;
 
     DynamicList<floatScalar> fField(pTraits<Type>::nComponents*nTotPoints);
@@ -139,7 +141,8 @@ void Foam::writeFuns::write
     const labelList& addPointCellLabels = topo.addPointCellLabels();
     const label nTotPoints = mesh.nPoints() + addPointCellLabels.size();
 
-    os  << vvf.name() << ' ' << pTraits<Type>::nComponents << ' '
+    os  << vvf.name() << ' '
+        << int(pTraits<Type>::nComponents) << ' '
         << nTotPoints << " float" << std::endl;
 
     DynamicList<floatScalar> fField(pTraits<Type>::nComponents*nTotPoints);
diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C
index 8499dc509f176b50813dafbd9076da8d355545ce..6cf2c56d207ae3a326960a19fb1f96abf7cab630 100644
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -119,8 +119,9 @@ void writeVTKFields
     forAll(values, fieldI)
     {
         Info<< "        writing field " << fieldNames[fieldI] << endl;
-        os  << nl << fieldNames[fieldI] << ' ' << pTraits<Type>::nComponents
-            << ' ' << values[fieldI].size() << " float" << nl;
+        os  << nl << fieldNames[fieldI] << ' '
+            << int(pTraits<Type>::nComponents) << ' '
+            << values[fieldI].size() << " float" << nl;
         label offset = 0;
         forAll(addr, trackI)
         {
diff --git a/src/fileFormats/sampledSetWriters/vtk/vtkSetWriter.C b/src/fileFormats/sampledSetWriters/vtk/vtkSetWriter.C
index 4663b0878f46cf4eba3af28d12062661b9c4ab2f..5a111061826de06cb1c780c29d9cdf63972c7555 100644
--- a/src/fileFormats/sampledSetWriters/vtk/vtkSetWriter.C
+++ b/src/fileFormats/sampledSetWriters/vtk/vtkSetWriter.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -85,7 +85,8 @@ void Foam::vtkSetWriter<Type>::write
 
     forAll(valueSetNames, setI)
     {
-        os  << valueSetNames[setI] << ' ' << pTraits<Type>::nComponents << ' '
+        os  << valueSetNames[setI] << ' '
+            << int(pTraits<Type>::nComponents) << ' '
             << points.size() << " float" << nl;
 
         const Field<Type>& fld = *valueSets[setI];
@@ -169,7 +170,8 @@ void Foam::vtkSetWriter<Type>::write
 
     forAll(valueSetNames, setI)
     {
-        os  << valueSetNames[setI] << ' ' << pTraits<Type>::nComponents << ' '
+        os  << valueSetNames[setI] << ' '
+            << int(pTraits<Type>::nComponents) << ' '
             << nPoints << " float" << nl;
 
         const List<Field<Type>>& fieldVals = valueSets[setI];
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H
index a24b8d0ac478a389ad725772d02f8c00834541cf..b21a63b2ae5a1693387ab3e68e8435863f7528b1 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H
@@ -66,7 +66,6 @@ Description
         Property     | Description             | Required    | Default value
         p            | pressure field name     | no          | p
         cyclicPatch  | cyclic patch name       | yes         |
-        orientation  | 1 to open or -1 to close | yes|
         openFraction | current open fraction [0-1] | yes |
         openingTime  | time taken to open or close the baffle | yes |
         maxOpenFractionDelta | max fraction change per timestep | yes |
@@ -83,7 +82,6 @@ Description
         type            activePressureForceBaffleVelocity;
         p               p;
         cyclicPatch     cyclic1;
-        orientation     1;
         openFraction    0.2;
         openingTime     5.0;
         maxOpenFractionDelta 0.1;
diff --git a/src/fvAgglomerationMethods/Allwmake b/src/fvAgglomerationMethods/Allwmake
index f5cfa2cc5b88ecd698f30d349e0cc06c8892b5b7..ebc2bd18e39f076ad767d4b60c073b5c6bbfbf43 100755
--- a/src/fvAgglomerationMethods/Allwmake
+++ b/src/fvAgglomerationMethods/Allwmake
@@ -15,17 +15,4 @@ fi
 
 wmake $targetType pairPatchAgglomeration
 
-
-## get SCOTCH_VERSION, SCOTCH_ARCH_PATH
-if settings=`$WM_PROJECT_DIR/bin/foamEtcFile config.sh/scotch`
-then
-    . $settings
-    echo "using SCOTCH_ARCH_PATH=$SCOTCH_ARCH_PATH"
-else
-    echo
-    echo "Error: no config.sh/scotch settings"
-    echo
-fi
-
-
 #------------------------------------------------------------------------------
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
index 5103acc332fdad8d74d62f029b325b4541172afd..b86dd760cde7e87a325bab02f08a51ac8e4919f8 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
@@ -66,7 +66,7 @@ Foam::cloudSolution::cloudSolution(const fvMesh& mesh, const dictionary& dict)
 
         // transient default to false asks for extra massFlowRate
         // in transient lagrangian
-        dict_.lookup("transient") >> transient_;
+        transient_ = true;
     }
 }
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudToVTK/CloudToVTK.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudToVTK/CloudToVTK.C
index 68c71b7f7de9adc4abd91b2557b410553c69c475..3c0f9a93cfa5c38e137e06871d84f9b0299fd3e7 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudToVTK/CloudToVTK.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudToVTK/CloudToVTK.C
@@ -65,7 +65,8 @@ void Foam::CloudToVTK<CloudType>::writeFieldData
 ) const
 {
     vtkOs
-        << title << ' ' << pTraits<Type>::nComponents << ' '
+        << title << ' '
+        << int(pTraits<Type>::nComponents) << ' '
         << nParcels << " float" << std::endl;
     writeData(vtkOs, binary, data);
 }
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
index 50523060f95c2f450c7da2c8fc2257a812060a25..c0b2724da1d11e67dae850b20a2e4916b1ce5482 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
@@ -669,6 +669,7 @@ private:
             //- Determines cell zone from cell region information.
             bool calcRegionToZone
             (
+                const label backgroundZoneID,
                 const label surfZoneI,
                 const label ownRegion,
                 const label neiRegion,
@@ -680,6 +681,7 @@ private:
             //  marked in namedSurfaceIndex regarded as blocked.
             void findCellZoneTopo
             (
+                const label backgroundZoneID,
                 const pointField& locationsInMesh,
                 const labelList& allSurfaceIndex,
                 const labelList& namedSurfaceIndex,
@@ -700,6 +702,7 @@ private:
             void zonify
             (
                 const bool allowFreeStandingZoneFaces,
+                const label backgroundZoneID,
                 const pointField& locationsInMesh,
                 const wordList& zonesInMesh,
 
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
index 8f6986d682f51ba655765f4c1b8e35233a185c3d..ebc759119fea5d2f686c31b9524e748bc8e0fcda 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -311,6 +311,7 @@ void Foam::meshRefinement::getBafflePatches
         zonify
         (
             true,               // allowFreeStandingZoneFaces
+            -2,                 // zone to put unreached cells into
             locationsInMesh,
             zonesInMesh,
 
@@ -1810,8 +1811,75 @@ void Foam::meshRefinement::findCellZoneInsideWalk
 }
 
 
+bool Foam::meshRefinement::calcRegionToZone
+(
+    const label backgroundZoneID,
+    const label surfZoneI,
+    const label ownRegion,
+    const label neiRegion,
+
+    labelList& regionToCellZone
+) const
+{
+    bool changed = false;
+
+    // Check whether inbetween different regions
+    if (ownRegion != neiRegion)
+    {
+        // Jump. Change one of the sides to my type.
+
+        // 1. Interface between my type and unset region.
+        // Set region to keepRegion
+
+        if (regionToCellZone[ownRegion] == -2)
+        {
+            if (regionToCellZone[neiRegion] == surfZoneI)
+            {
+                // Face between unset and my region. Put unset
+                // region into keepRegion
+                //MEJ: see comment in findCellZoneTopo
+                if (backgroundZoneID != -2)
+                {
+                    regionToCellZone[ownRegion] = backgroundZoneID;
+                    changed = true;
+                }
+            }
+            else if (regionToCellZone[neiRegion] != -2)
+            {
+                // Face between unset and other region.
+                // Put unset region into my region
+                regionToCellZone[ownRegion] = surfZoneI;
+                changed = true;
+            }
+        }
+        else if (regionToCellZone[neiRegion] == -2)
+        {
+            if (regionToCellZone[ownRegion] == surfZoneI)
+            {
+                // Face between unset and my region. Put unset
+                // region into keepRegion
+                if (backgroundZoneID != -2)
+                {
+                    regionToCellZone[neiRegion] = backgroundZoneID;
+                    changed = true;
+                }
+            }
+            else if (regionToCellZone[ownRegion] != -2)
+            {
+                // Face between unset and other region.
+                // Put unset region into my region
+                regionToCellZone[neiRegion] = surfZoneI;
+                changed = true;
+            }
+        }
+    }
+    return changed;
+}
+
+
 void Foam::meshRefinement::findCellZoneTopo
 (
+    const label backgroundZoneID,
     const pointField& locationsInMesh,
     const labelList& allSurfaceIndex,
     const labelList& namedSurfaceIndex,
@@ -1825,12 +1893,17 @@ void Foam::meshRefinement::findCellZoneTopo
     // know which side of the face it relates to. So all we are doing here
     // is get the correspondence between surface/cellZone and regionSplit
     // region.
-    // This can occasionally go wrong when surfaces pass through a
-    // cell centre and a cell completely gets blocked off so becomes a
-    // separate region. The problem is to distinguish between a cell which
-    // properly should be deleted (so in the 'background' mesh)
-    // or just a single-cell region from incorrect baffling. Currently the
-    // logic is just to detect cells that are on different regions.
+    // See the logic in calcRegionToZone. The problem is what to do
+    // with unreachable parts of the mesh (cellToZone = -2).
+    // In 23x this routine was only called for actually zoneing an existing
+    // mesh so we had to do something with these cells and they
+    // would get set to -1 (background). However in this version this routine
+    // also gets used much earlier on when after the surface refinement it
+    // removes unreachable cells ('Removing mesh beyond surface intersections')
+    // and this is when we keep -2 so it gets removed.
+    // So the zone to set these unmarked cells to is provided as argument:
+    // - backgroundZoneID = -2 : do not change so remove cells
+    // - backgroundZoneID = -1 : put into background
 
 
     // Assumes:
@@ -1920,107 +1993,95 @@ void Foam::meshRefinement::findCellZoneTopo
         }
     }
 
-    // Synchronise any changes due to locationInMesh 
-    Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
-    Pstream::listCombineScatter(regionToCellZone);
-
-    // Get coupled neighbour cellRegion
-    labelList neiCellRegion;
-    syncTools::swapBoundaryCellList(mesh_, cellRegion, neiCellRegion);
 
+    // Find correspondence between cell zone and surface
+    // by changing cell zone every time we cross a surface.
+    while (true)
+    {
+        // Synchronise regionToCellZone.
+        // Note:
+        // - region numbers are identical on all processors
+        // - keepRegion is identical ,,
+        // - cellZones are identical ,,
+        // This done at top of loop to account for geometric matching
+        // not being synchronised.
+        Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
+        Pstream::listCombineScatter(regionToCellZone);
 
 
+        bool changed = false;
 
-    // Detect unassigned cell (region) bordering two cellZones. Put these
-    // cells into either neighbouring cellZone.
-    labelList growCellToZone(mesh_.nCells(), -2);
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-    {
-        label own = mesh_.faceOwner()[faceI];
-        label ownRegion = cellRegion[own];
-        label nei = mesh_.faceNeighbour()[faceI];
-        label neiRegion = cellRegion[nei];
+        // Internal faces
 
-        // Check whether inbetween different regions
-        if (ownRegion != neiRegion)
+        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
         {
-            if (regionToCellZone[ownRegion] == -2)
-            {
-                if (regionToCellZone[neiRegion] != -2)
-                {
-                    if (growCellToZone[own] == -2)
-                    {
-                        // First visit of cell
-                        growCellToZone[own] = regionToCellZone[neiRegion];
-                    }
-                    else if (regionToCellZone[neiRegion] != growCellToZone[own])
-                    {
-                        // Found a cell bordering multiple cellZones. Grow.
-                        cellToZone[own] = growCellToZone[own];
-                    }
-                }
-            }
-            else if (regionToCellZone[neiRegion] == -2)
+            label surfI = namedSurfaceIndex[faceI];
+
+            // Connected even if no cellZone defined for surface
+            if (surfI != -1)
             {
-                if (growCellToZone[nei] == -2)
-                {
-                    // First visit of cell
-                    growCellToZone[nei] = regionToCellZone[ownRegion];
-                }
-                else if (regionToCellZone[ownRegion] != growCellToZone[nei])
-                {
-                    cellToZone[nei] = growCellToZone[nei];
-                }
+                // Calculate region to zone from cellRegions on either side
+                // of internal face.
+                bool changedCell = calcRegionToZone
+                (
+                    backgroundZoneID,
+                    surfaceToCellZone[surfI],
+                    cellRegion[mesh_.faceOwner()[faceI]],
+                    cellRegion[mesh_.faceNeighbour()[faceI]],
+                    regionToCellZone
+                );
+
+                changed = changed | changedCell;
             }
         }
-    }
 
+        // Boundary faces
 
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
+        // Get coupled neighbour cellRegion
+        labelList neiCellRegion;
+        syncTools::swapBoundaryCellList(mesh_, cellRegion, neiCellRegion);
 
-        if (pp.coupled())
+        // Calculate region to zone from cellRegions on either side of coupled
+        // face.
+        forAll(patches, patchI)
         {
-            forAll(pp, i)
+            const polyPatch& pp = patches[patchI];
+
+            if (pp.coupled())
             {
-                label faceI = pp.start()+i;
+                forAll(pp, i)
+                {
+                    label faceI = pp.start()+i;
 
-                label own = mesh_.faceOwner()[faceI];
-                label ownRegion = cellRegion[own];
-                label neiRegion = neiCellRegion[faceI-mesh_.nInternalFaces()];
+                    label surfI = namedSurfaceIndex[faceI];
 
-                // Check whether inbetween different regions
-                if (ownRegion != neiRegion)
-                {
-                    if (regionToCellZone[ownRegion] == -2)
+                    // Connected even if no cellZone defined for surface
+                    if (surfI != -1)
                     {
-                        if (regionToCellZone[neiRegion] != -2)
-                        {
-                            if (growCellToZone[own] == -2)
-                            {
-                                // First visit of cell
-                                growCellToZone[own] =
-                                    regionToCellZone[neiRegion];
-                            }
-                            else if
-                            (
-                                regionToCellZone[neiRegion]
-                             != growCellToZone[own]
-                            )
-                            {
-                                // Found a cell bordering multiple cellZones.
-                                cellToZone[own] = growCellToZone[own];
-                            }
-                        }
+                        bool changedCell = calcRegionToZone
+                        (
+                            backgroundZoneID,
+                            surfaceToCellZone[surfI],
+                            cellRegion[mesh_.faceOwner()[faceI]],
+                            neiCellRegion[faceI-mesh_.nInternalFaces()],
+                            regionToCellZone
+                        );
+
+                        changed = changed | changedCell;
                     }
                 }
             }
         }
+
+        if (!returnReduce(changed, orOp<bool>()))
+        {
+            break;
+        }
     }
 
+
     if (debug)
     {
         forAll(regionToCellZone, regionI)
@@ -2264,6 +2325,7 @@ void Foam::meshRefinement::getIntersections
 void Foam::meshRefinement::zonify
 (
     const bool allowFreeStandingZoneFaces,
+    const label backgroundZoneID,
     const pointField& locationsInMesh,
     const wordList& zonesInMesh,
 
@@ -2486,25 +2548,23 @@ void Foam::meshRefinement::zonify
     }
 
 
-    // 5. Find any unassigned regions (from regionSplit). This will
-    //    just walk slightly out to cover those single cells that
-    //    are baffled off into a separate region.
-
-    Info<< "Growing from known cellZones to fix cells inbetween cellZones"
-        << nl << endl;
-
-    findCellZoneTopo
-    (
-        pointField(0),
-        globalRegion1,      // To split up cells
-        namedSurfaceIndex,  // Step across named surfaces to propagate
-        surfaceToCellZone,
-        cellToZone
-    );
-
+    // 5. Find any unassigned regions (from regionSplit)
 
     if (namedSurfaces.size())
     {
+        Info<< "Walking from known cellZones; crossing a faceZone "
+            << "face changes cellZone" << nl << endl;
+
+        findCellZoneTopo
+        (
+            backgroundZoneID,
+            pointField(0),
+            globalRegion1,      // To split up cells
+            namedSurfaceIndex,  // Step across named surfaces to propagate
+            surfaceToCellZone,
+            cellToZone
+        );
+
         // Make sure namedSurfaceIndex is unset inbetween same cell zones.
         if (!allowFreeStandingZoneFaces)
         {
@@ -4269,6 +4329,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     zonify
     (
         allowFreeStandingZoneFaces,
+        -1,                 // Set all cells with cellToZone -2 to -1
         locationsInMesh,
         zonesInMesh,
 
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
index 3dcef1dea1c29800a9fad6fc0994d99e70a1534c..a8ca96b2912d9a86df2aa05ef38772f5f8c3894d 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
@@ -119,7 +119,7 @@ Foam::solidChemistryModel<CompType, SolidThermo>::RRsHs() const
 
     if (this->chemistry_)
     {
-        DimensionedField<scalar, volMesh>& RRs = tRRsHs();
+        DimensionedField<scalar, volMesh>& RRs = tRRsHs.ref();
 
         const volScalarField& T =  this->solidThermo().T();
         const volScalarField& p =  this->solidThermo().p();
diff --git a/tutorials/incompressible/simpleFoam/rotorDisk/Allrun b/tutorials/incompressible/simpleFoam/rotorDisk/Allrun
index 14d69c2a9a5ba99af346359cf50848ff1b8b47dc..f6f1258231271bbe7ef6e65f55c27603b3c0c01e 100755
--- a/tutorials/incompressible/simpleFoam/rotorDisk/Allrun
+++ b/tutorials/incompressible/simpleFoam/rotorDisk/Allrun
@@ -8,7 +8,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 runApplication blockMesh
 runApplication surfaceFeatureExtract
 runApplication snappyHexMesh -overwrite
-runApplication createPatch -overwrite
 
 runApplication $(getApplication)
 
diff --git a/tutorials/incompressible/simpleFoam/rotorDisk/system/snappyHexMeshDict b/tutorials/incompressible/simpleFoam/rotorDisk/system/snappyHexMeshDict
index 81afb27c643c7a2eda365f519cbd52c5dc761c46..35b0bb756f2bd3335c0eb2e998b26354a1fd5d3c 100644
--- a/tutorials/incompressible/simpleFoam/rotorDisk/system/snappyHexMeshDict
+++ b/tutorials/incompressible/simpleFoam/rotorDisk/system/snappyHexMeshDict
@@ -101,8 +101,9 @@ castellatedMeshControls
         }
     }
 
-    locationInMesh (1e-5 1e-5 1e-5); // Offset from (0 0 0) to avoid
-                                     // coinciding with face or edge
+    locationInMesh (1e-5 -1e-2 1e-5);// Offset from (0 0 0) to avoid
+                                     // coinciding with face or edge and keep
+                                     // away from disk itself
 }
 
 snapControls
@@ -120,9 +121,6 @@ addLayersControls
     expansionRatio      1.2;
     finalLayerThickness 0.5;
     minThickness        1e-3;
-//  firstLayerThickness 0.01;
-
-//  maxThicknessToMedialRatio 0.6;
 }
 
 meshQualityControls
@@ -130,12 +128,6 @@ meshQualityControls
 //    minTetQuality -1e+30;
 }
 
-writeFlags
-(
-    scalarLevels
-    layerSets
-    layerFields
-);
 
 mergeTolerance 1e-6;
 
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoiling/system/fvSolution b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoiling/system/fvSolution
index ff61cafc76f65fcbed5b3ddad02985fe94822998..0a396b99cabad8a3c2cf82861fcbf76b7007fbc0 100644
--- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoiling/system/fvSolution
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoiling/system/fvSolution
@@ -35,7 +35,7 @@ solvers
         agglomerator    faceAreaPair;
         mergeLevels     1;
         tolerance       1e-8;
-        relTol          0.1;
+        relTol          0.01;
         maxIter         100;
         minIter         2;
     }
@@ -98,7 +98,7 @@ relaxationFactors
 {
     fields
     {
-        iDmdt 0.4;
+        iDmdt           0.1;
     }
 
     equations