diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
index db59fd4fc3..81f81dc508 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
@@ -268,8 +268,6 @@ void Foam::meshRefinement::calcCellCellRays
 // Returns first intersection if there are more than one.
 void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 {
-    const pointField& cellCentres = mesh_.cellCentres();
-
     // Stats on edges to test. Count proc faces only once.
     PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
 
@@ -309,28 +307,17 @@ void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
     // Collect segments we want to test for
     pointField start(changedFaces.size());
     pointField end(changedFaces.size());
-
-    forAll(changedFaces, i)
     {
-        label faceI = changedFaces[i];
-        label own = mesh_.faceOwner()[faceI];
-
-        start[i] = cellCentres[own];
-        if (mesh_.isInternalFace(faceI))
-        {
-            end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
-        }
-        else
-        {
-            end[i] = neiCc[faceI-mesh_.nInternalFaces()];
-        }
-    }
-
-    // Extend segments a bit
-    {
-        const vectorField smallVec(ROOTSMALL*(end-start));
-        start -= smallVec;
-        end += smallVec;
+        labelList minLevel;
+        calcCellCellRays
+        (
+            neiCc,
+            neiLevel,
+            changedFaces,
+            start,
+            end,
+            minLevel
+        );
     }
 
 
@@ -2964,8 +2951,6 @@ void Foam::meshRefinement::dumpRefinementLevel() const
 void Foam::meshRefinement::dumpIntersections(const fileName& prefix) const
 {
     {
-        const pointField& cellCentres = mesh_.cellCentres();
-
         OFstream str(prefix + "_edges.obj");
         label vertI = 0;
         Pout<< "meshRefinement::dumpIntersections :"
@@ -2986,27 +2971,17 @@ void Foam::meshRefinement::dumpIntersections(const fileName& prefix) const
         // Collect segments we want to test for
         pointField start(intersectionFaces.size());
         pointField end(intersectionFaces.size());
-
-        forAll(intersectionFaces, i)
         {
-            label faceI = intersectionFaces[i];
-            start[i] = cellCentres[mesh_.faceOwner()[faceI]];
-
-            if (mesh_.isInternalFace(faceI))
-            {
-                end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
-            }
-            else
-            {
-                end[i] = neiCc[faceI-mesh_.nInternalFaces()];
-            }
-        }
-
-        // Extend segments a bit
-        {
-            const vectorField smallVec(ROOTSMALL*(end-start));
-            start -= smallVec;
-            end += smallVec;
+            labelList minLevel;
+            calcCellCellRays
+            (
+                neiCc,
+                labelList(neiCc.size(), -1),
+                intersectionFaces,
+                start,
+                end,
+                minLevel
+            );
         }
 
 
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
index c0b2724da1..50523060f9 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
@@ -669,7 +669,6 @@ private:
             //- Determines cell zone from cell region information.
             bool calcRegionToZone
             (
-                const label backgroundZoneID,
                 const label surfZoneI,
                 const label ownRegion,
                 const label neiRegion,
@@ -681,7 +680,6 @@ private:
             //  marked in namedSurfaceIndex regarded as blocked.
             void findCellZoneTopo
             (
-                const label backgroundZoneID,
                 const pointField& locationsInMesh,
                 const labelList& allSurfaceIndex,
                 const labelList& namedSurfaceIndex,
@@ -702,7 +700,6 @@ 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 9af9efa7d8..8f6986d682 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -311,7 +311,6 @@ void Foam::meshRefinement::getBafflePatches
         zonify
         (
             true,               // allowFreeStandingZoneFaces
-            -2,                 // zone to put unreached cells into
             locationsInMesh,
             zonesInMesh,
 
@@ -1811,75 +1810,8 @@ 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,
@@ -1893,17 +1825,12 @@ 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.
-    // 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
+    // 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.
 
 
     // Assumes:
@@ -1993,95 +1920,107 @@ void Foam::meshRefinement::findCellZoneTopo
         }
     }
 
+    // Synchronise any changes due to locationInMesh 
+    Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
+    Pstream::listCombineScatter(regionToCellZone);
 
-    // Find correspondence between cell zone and surface
-    // by changing cell zone every time we cross a surface.
-    while (true)
+    // Get coupled neighbour cellRegion
+    labelList neiCellRegion;
+    syncTools::swapBoundaryCellList(mesh_, cellRegion, neiCellRegion);
+
+
+
+
+    // 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++)
     {
-        // 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);
+        label own = mesh_.faceOwner()[faceI];
+        label ownRegion = cellRegion[own];
+        label nei = mesh_.faceNeighbour()[faceI];
+        label neiRegion = cellRegion[nei];
 
-
-        bool changed = false;
-
-        // Internal faces
-
-        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+        // Check whether inbetween different regions
+        if (ownRegion != neiRegion)
         {
-            label surfI = namedSurfaceIndex[faceI];
-
-            // Connected even if no cellZone defined for surface
-            if (surfI != -1)
+            if (regionToCellZone[ownRegion] == -2)
             {
-                // 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;
+                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)
+            {
+                if (growCellToZone[nei] == -2)
+                {
+                    // First visit of cell
+                    growCellToZone[nei] = regionToCellZone[ownRegion];
+                }
+                else if (regionToCellZone[ownRegion] != growCellToZone[nei])
+                {
+                    cellToZone[nei] = growCellToZone[nei];
+                }
             }
         }
+    }
 
-        // Boundary faces
 
-        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-        // Get coupled neighbour cellRegion
-        labelList neiCellRegion;
-        syncTools::swapBoundaryCellList(mesh_, cellRegion, neiCellRegion);
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
 
-        // Calculate region to zone from cellRegions on either side of coupled
-        // face.
-        forAll(patches, patchI)
+        if (pp.coupled())
         {
-            const polyPatch& pp = patches[patchI];
-
-            if (pp.coupled())
+            forAll(pp, i)
             {
-                forAll(pp, i)
+                label faceI = pp.start()+i;
+
+                label own = mesh_.faceOwner()[faceI];
+                label ownRegion = cellRegion[own];
+                label neiRegion = neiCellRegion[faceI-mesh_.nInternalFaces()];
+
+                // Check whether inbetween different regions
+                if (ownRegion != neiRegion)
                 {
-                    label faceI = pp.start()+i;
-
-                    label surfI = namedSurfaceIndex[faceI];
-
-                    // Connected even if no cellZone defined for surface
-                    if (surfI != -1)
+                    if (regionToCellZone[ownRegion] == -2)
                     {
-                        bool changedCell = calcRegionToZone
-                        (
-                            backgroundZoneID,
-                            surfaceToCellZone[surfI],
-                            cellRegion[mesh_.faceOwner()[faceI]],
-                            neiCellRegion[faceI-mesh_.nInternalFaces()],
-                            regionToCellZone
-                        );
-
-                        changed = changed | changedCell;
+                        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];
+                            }
+                        }
                     }
                 }
             }
         }
-
-        if (!returnReduce(changed, orOp<bool>()))
-        {
-            break;
-        }
     }
 
-
     if (debug)
     {
         forAll(regionToCellZone, regionI)
@@ -2200,10 +2139,6 @@ void Foam::meshRefinement::getIntersections
     PackedBoolList& posOrientation
 ) const
 {
-    const pointField& cellCentres = mesh_.cellCentres();
-    const labelList& faceOwner = mesh_.faceOwner();
-    const labelList& faceNeighbour = mesh_.faceNeighbour();
-
     namedSurfaceIndex.setSize(mesh_.nFaces());
     namedSurfaceIndex = -1;
 
@@ -2219,30 +2154,19 @@ void Foam::meshRefinement::getIntersections
     pointField start(testFaces.size());
     pointField end(testFaces.size());
 
-    forAll(testFaces, i)
     {
-        label faceI = testFaces[i];
-
-        if (mesh_.isInternalFace(faceI))
-        {
-            start[i] = cellCentres[faceOwner[faceI]];
-            end[i] = cellCentres[faceNeighbour[faceI]];
-        }
-        else
-        {
-            start[i] = cellCentres[faceOwner[faceI]];
-            end[i] = neiCc[faceI-mesh_.nInternalFaces()];
-        }
+        labelList minLevel;
+        calcCellCellRays
+        (
+            neiCc,
+            labelList(neiCc.size(), -1),
+            testFaces,
+            start,
+            end,
+            minLevel
+        );
     }
 
-    // Extend segments a bit
-    {
-        const vectorField smallVec(ROOTSMALL*(end-start));
-        start -= smallVec;
-        end += smallVec;
-    }
-
-
     // Do test for intersections
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
     // Note that we intersect all intersected faces again. Could reuse
@@ -2340,7 +2264,6 @@ void Foam::meshRefinement::getIntersections
 void Foam::meshRefinement::zonify
 (
     const bool allowFreeStandingZoneFaces,
-    const label backgroundZoneID,
     const pointField& locationsInMesh,
     const wordList& zonesInMesh,
 
@@ -2535,24 +2458,53 @@ void Foam::meshRefinement::zonify
         );
     }
 
+    if (debug&MESH)
+    {
+        Pout<< "Writing cell zone allocation on mesh to time "
+            << timeName() << endl;
+        volScalarField volCellToZone
+        (
+            IOobject
+            (
+                "cellToZoneBeforeTopo",
+                timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimless, 0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+
+        forAll(cellToZone, cellI)
+        {
+            volCellToZone[cellI] = cellToZone[cellI];
+        }
+        volCellToZone.write();
+    }
+
+
+    // 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)
         {
@@ -3621,6 +3573,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
     if (debug&MESH)
     {
+        runTime++;
         Pout<< "Writing subsetted mesh to time " << timeName()
             << endl;
         write
@@ -4316,7 +4269,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     zonify
     (
         allowFreeStandingZoneFaces,
-        -1,                 // Set all cells with cellToZone -2 to -1
         locationsInMesh,
         zonesInMesh,
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
index 484229bf93..d5e45f3d47 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
@@ -1316,6 +1316,8 @@ void Foam::snappyRefineDriver::removeInsideCells
 
     if (debug&meshRefinement::MESH)
     {
+        const_cast<Time&>(mesh.time())++;
+
         Pout<< "Writing subsetted mesh to time "
             << meshRefiner_.timeName() << endl;
         meshRefiner_.write