diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTools.C b/applications/utilities/mesh/manipulation/checkMesh/checkTools.C
index a77ca272cd..00c1f4b0f2 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTools.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTools.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2021 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -423,8 +423,8 @@ void Foam::mergeAndWrite
 {
     const polyMesh& mesh = refCast<const polyMesh>(set.db());
 
-    pointField mergedPts;
-    labelList mergedIDs;
+    labelField mergedIDs(set.sortedToc());
+    pointField mergedPts(mesh.points(), mergedIDs);
 
     if (Pstream::parRun())
     {
@@ -432,55 +432,12 @@ void Foam::mergeAndWrite
         // (mesh.globalData().mergePoints etc) since this might
         // hide any synchronisation problem
 
-        globalIndex globalNumbering(mesh.nPoints());
+        // Renumber local ids -> global ids
+        globalIndex(mesh.nPoints()).inplaceToGlobal(mergedIDs);
 
-        mergedPts.setSize(returnReduce(set.size(), sumOp<label>()));
-        mergedIDs.setSize(mergedPts.size());
-
-        labelList setPointIDs(set.sortedToc());
-
-        // Get renumbered local data
-        pointField myPoints(mesh.points(), setPointIDs);
-        labelList myIDs(globalNumbering.toGlobal(setPointIDs));
-
-        if (Pstream::master())
-        {
-            // Insert master data first
-            label pOffset = 0;
-            SubList<point>(mergedPts, myPoints.size(), pOffset) = myPoints;
-            SubList<label>(mergedIDs, myIDs.size(), pOffset) = myIDs;
-            pOffset += myPoints.size();
-
-            // Receive slave ones
-            for (const int slave : Pstream::subProcs())
-            {
-                IPstream fromSlave(Pstream::commsTypes::scheduled, slave);
-
-                pointField slavePts(fromSlave);
-                labelList slaveIDs(fromSlave);
-
-                SubList<point>(mergedPts, slavePts.size(), pOffset) = slavePts;
-                SubList<label>(mergedIDs, slaveIDs.size(), pOffset) = slaveIDs;
-                pOffset += slaveIDs.size();
-            }
-        }
-        else
-        {
-            // Construct processor stream with estimate of size. Could
-            // be improved.
-            OPstream toMaster
-            (
-                Pstream::commsTypes::scheduled,
-                Pstream::masterNo(),
-                myPoints.byteSize() + myIDs.byteSize()
-            );
-            toMaster << myPoints << myIDs;
-        }
-    }
-    else
-    {
-        mergedIDs = set.sortedToc();
-        mergedPts = pointField(mesh.points(), mergedIDs);
+        globalIndex gatherer(mergedIDs.size(), globalIndex::gatherOnly{});
+        gatherer.gatherInplace(mergedIDs);
+        gatherer.gatherInplace(mergedPts);
     }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
index 86c48920ad..6e22d5140a 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -136,7 +136,7 @@ void Foam::PatchTools::gatherAndMerge
 
     if (Pstream::parRun())
     {
-        // Renumber the setPatch points/faces into unique points
+        // Renumber the points/faces into unique points
         globalPointsPtr = mesh.globalData().mergePoints
         (
             meshPoints,
@@ -147,65 +147,22 @@ void Foam::PatchTools::gatherAndMerge
 
         globalFacesPtr.reset(new globalIndex(localFaces.size()));
 
-        if (Pstream::master())
+        // Renumber faces locally
+        List<FaceType> myFaces(localFaces);
+        for (auto& f : myFaces)
         {
-            // Get renumbered local data
-            pointField myPoints(mesh.points(), uniqueMeshPointLabels);
-            List<FaceType> myFaces(localFaces);
-            for (auto& f : myFaces)
-            {
-                inplaceRenumber(pointToGlobal, f);
-            }
-
-
-            mergedFaces.setSize(globalFacesPtr().totalSize());
-            mergedPoints.setSize(globalPointsPtr().totalSize());
-
-            // Insert master data first
-            label pOffset = globalPointsPtr().localStart(Pstream::masterNo());
-            SubList<point>(mergedPoints, myPoints.size(), pOffset) = myPoints;
-
-            label fOffset = globalFacesPtr().localStart(Pstream::masterNo());
-            SubList<FaceType>(mergedFaces, myFaces.size(), fOffset) = myFaces;
-
-
-            // Receive slave ones
-            for (const int slave : Pstream::subProcs())
-            {
-                IPstream fromSlave(Pstream::commsTypes::scheduled, slave);
-
-                pointField slavePoints(fromSlave);
-                List<FaceType> slaveFaces(fromSlave);
-
-                label pOffset = globalPointsPtr().localStart(slave);
-                SubList<point>(mergedPoints, slavePoints.size(), pOffset) =
-                    slavePoints;
-
-                label fOffset = globalFacesPtr().localStart(slave);
-                SubList<FaceType>(mergedFaces, slaveFaces.size(), fOffset) =
-                    slaveFaces;
-            }
+            inplaceRenumber(pointToGlobal, f);
         }
-        else
-        {
-            // Get renumbered local data
-            pointField myPoints(mesh.points(), uniqueMeshPointLabels);
-            List<FaceType> myFaces(localFaces);
-            for (auto& f : myFaces)
-            {
-                inplaceRenumber(pointToGlobal, f);
-            }
 
-            // Construct processor stream with estimate of size. Could
-            // be improved.
-            OPstream toMaster
-            (
-                Pstream::commsTypes::scheduled,
-                Pstream::masterNo(),
-                myPoints.byteSize() + 4*sizeof(label)*myFaces.size()
-            );
-            toMaster << myPoints << myFaces;
-        }
+        // Can also use
+        //     UIndirectList<point>(mesh.points(), uniqueMeshPointLabels)
+        // but favour communication over local memory use
+        globalPointsPtr().gather
+        (
+            pointField(mesh.points(), uniqueMeshPointLabels),
+            mergedPoints
+        );
+        globalFacesPtr().gather(myFaces, mergedFaces);
     }
     else
     {
diff --git a/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.H b/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.H
index bb1f014163..bf78c77516 100644
--- a/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.H
+++ b/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.H
@@ -35,11 +35,11 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SprayCloud_H
-#define SprayCloud_H
+#ifndef Foam_SprayCloud_H
+#define Foam_SprayCloud_H
 
 #include "sprayCloud.H"
-#include "SortableList.H"
+#include "SortList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H b/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H
index 4c1cbffffb..4ce6c576bc 100644
--- a/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H
+++ b/src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H
@@ -88,95 +88,73 @@ inline Foam::scalar Foam::SprayCloud<CloudType>::penetration
             << exit(FatalError);
     }
 
-    scalar distance = 0.0;
-
     const label nParcel = this->size();
-    globalIndex globalParcels(nParcel);
-    const label nParcelSum = globalParcels.totalSize();
+    const globalIndex globalParcels(nParcel);
+    const label nTotParcel = globalParcels.totalSize();
 
-    if (nParcelSum == 0)
+    if (nTotParcel == 0)
     {
-        return distance;
+        return 0;
     }
 
     // lists of parcels mass and distance from initial injection point
-    List<List<scalar>> procMass(Pstream::nProcs());
-    List<List<scalar>> procDist(Pstream::nProcs());
+    List<scalar> mass(nParcel);
+    List<scalar> dist(nParcel);
 
-    List<scalar>& mass = procMass[Pstream::myProcNo()];
-    List<scalar>& dist = procDist[Pstream::myProcNo()];
-
-    mass.setSize(nParcel);
-    dist.setSize(nParcel);
-
-    label i = 0;
-    scalar mSum = 0.0;
-    for (const parcelType& p : *this)
+    scalar mTotal = 0;
     {
-        scalar m = p.nParticle()*p.mass();
-        scalar d = mag(p.position() - p.position0());
-        mSum += m;
+        label i = 0;
+        for (const parcelType& p : *this)
+        {
+            scalar m = p.nParticle()*p.mass();
+            scalar d = mag(p.position() - p.position0());
+            mTotal += m;
 
-        mass[i] = m;
-        dist[i] = d;
-
-        ++i;
+            mass[i] = m;
+            dist[i] = d;
+            ++i;
+        }
     }
+    // Total mass across all processors
+    reduce(mTotal, sumOp<scalar>());
 
-    // calculate total mass across all processors
-    reduce(mSum, sumOp<scalar>());
-    Pstream::gatherList(procMass);
-    Pstream::gatherList(procDist);
+    scalar distance = 0;
+    globalParcels.gatherInplace(mass);
+    globalParcels.gatherInplace(dist);
 
     if (Pstream::master())
     {
-        // flatten the mass lists
-        List<scalar> allMass(nParcelSum, Zero);
-        SortableList<scalar> allDist(nParcelSum, Zero);
-        for (const int proci : Pstream::allProcs())
+        if (nTotParcel == 1)
         {
-            SubList<scalar>
-            (
-                allMass,
-                globalParcels.localSize(proci),
-                globalParcels.localStart(proci)
-            ) = procMass[proci];
-
-            // flatten the distance list
-            SubList<scalar>
-            (
-                allDist,
-                globalParcels.localSize(proci),
-                globalParcels.localStart(proci)
-            ) = procDist[proci];
+            distance = dist[0];
         }
-
-        // sort allDist distances into ascending order
-        // note: allMass masses are left unsorted
-        allDist.sort();
-
-        if (nParcelSum > 1)
+        else
         {
-            const scalar mLimit = fraction*mSum;
-            const labelList& indices = allDist.indices();
+            // Distances - sored into ascending order
+            // Masses - leave unsorted
 
-            if (mLimit > (mSum - allMass[indices.last()]))
+            const SortList<scalar> sortedDist(dist);
+
+            const scalar mLimit = fraction*mTotal;
+            const labelList& indices = sortedDist.indices();
+
+            if (mLimit > (mTotal - mass[indices.last()]))
             {
-                distance = allDist.last();
+                distance = sortedDist.last();
             }
             else
             {
                 // assuming that 'fraction' is generally closer to 1 than 0,
                 // loop through in reverse distance order
-                const scalar mThreshold = (1.0 - fraction)*mSum;
-                scalar mCurrent = 0.0;
+                const scalar mThreshold = (1.0 - fraction)*mTotal;
+                scalar mCurrent = 0;
                 label i0 = 0;
 
                 forAllReverse(indices, i)
                 {
                     label indI = indices[i];
 
-                    mCurrent += allMass[indI];
+                    mCurrent += mass[indI];
 
                     if (mCurrent > mThreshold)
                     {
@@ -187,21 +165,20 @@ inline Foam::scalar Foam::SprayCloud<CloudType>::penetration
 
                 if (i0 == indices.size() - 1)
                 {
-                    distance = allDist.last();
+                    distance = sortedDist.last();
                 }
                 else
                 {
                     // linearly interpolate to determine distance
-                    scalar alpha = (mCurrent - mThreshold)/allMass[indices[i0]];
+                    scalar alpha = (mCurrent - mThreshold)/mass[indices[i0]];
                     distance =
-                        allDist[i0] + alpha*(allDist[i0+1] - allDist[i0]);
+                    (
+                        sortedDist[i0]
+                      + alpha*(sortedDist[i0+1] - sortedDist[i0])
+                    );
                 }
             }
         }
-        else
-        {
-            distance = allDist.first();
-        }
     }
 
     Pstream::scatter(distance);