diff --git a/src/parallel/decompose/Allwclean b/src/parallel/decompose/Allwclean
index 402ddedc7a..bb8efaecaa 100755
--- a/src/parallel/decompose/Allwclean
+++ b/src/parallel/decompose/Allwclean
@@ -16,30 +16,39 @@ fi
 # Define how to clean an mpi-versioned library
 wcleanMpiLib()
 {
-    for libName
-    do
     (
+        version="$1"
+        libName="$2"
+
         WM_OPTIONS="$WM_OPTIONS$WM_MPLIB"
         whichmpi="$WM_PROJECT_DIR/platforms/$WM_OPTIONS/src/parallel/decompose/$libName/using:$FOAM_MPI"
-        whichscotch="$WM_PROJECT_DIR/platforms/$WM_OPTIONS/src/parallel/decompose/$libName/using:$SCOTCH_VERSION"
+        whichscotch="$WM_PROJECT_DIR/platforms/$WM_OPTIONS/src/parallel/decompose/$libName/using:$version"
         wclean $libName
     )
-    done
 }
 
+
 if [ -n "$SCOTCH_ARCH_PATH" ]
 then
     wclean scotchDecomp
 
     if [ -d "$FOAM_LIBBIN/$FOAM_MPI" ]
     then
-        wcleanMpiLib ptscotchDecomp
+        wcleanMpiLib $SCOTCH_VERSION ptscotchDecomp
     fi
 else
     echo "    skipping scotchDecomp (ptscotchDecomp)"
 fi
 
+
+if [ -n "$ZOLTAN_ARCH_PATH" ]
+then
+   wcleanMpiLib 3.90 zoltanDecomp
+fi
+
+
 wclean metisDecomp
 wclean decompositionMethods
 
+
 #------------------------------------------------------------------------------
diff --git a/src/parallel/decompose/Allwmake b/src/parallel/decompose/Allwmake
index 0c213e59c7..287a0922fd 100755
--- a/src/parallel/decompose/Allwmake
+++ b/src/parallel/decompose/Allwmake
@@ -22,33 +22,52 @@ fi
 # use sentinel file to handle version changes
 wmakeMpiLib()
 {
-    for libName
-    do
     (
+        version="$1"
+        libName="$2"
+
         WM_OPTIONS="$WM_OPTIONS$WM_MPLIB"
         whichmpi="$WM_PROJECT_DIR/platforms/$WM_OPTIONS/src/parallel/decompose/$libName/using:$FOAM_MPI"
-        whichscotch="$WM_PROJECT_DIR/platforms/$WM_OPTIONS/src/parallel/decompose/$libName/using:$SCOTCH_VERSION"
-        [ -e "$whichmpi" -a -e "$whichscotch" ] || wclean $libName
+        whichtarget="$WM_PROJECT_DIR/platforms/$WM_OPTIONS/src/parallel/decompose/$libName/using:$version"
+        [ -e "$whichmpi" -a -e "$whichtarget" ] || wclean $libName
         echo "wmake $targetType $libName"
         wmake $targetType $libName
-        touch "$whichmpi" "$whichscotch"
+        touch "$whichmpi" "$whichtarget"
     )
-    done
 }
 
+
 if [ -n "$SCOTCH_ARCH_PATH" ]
 then
     wmake $targetType scotchDecomp
 
     if [ -n "$WM_MPLIB" ] && [ -d "$FOAM_LIBBIN/$FOAM_MPI" ]
     then
-        wmakeMpiLib ptscotchDecomp
+        wmakeMpiLib $SCOTCH_VERSION ptscotchDecomp
     fi
 else
     echo "    skipping scotchDecomp (ptscotchDecomp)"
 fi
 
 
+if settings=`$WM_PROJECT_DIR/bin/foamEtcFile config.sh/zoltan`
+then
+    . $settings
+    echo "    using ZOLTAN_ARCH_PATH=$ZOLTAN_ARCH_PATH"
+
+    if [ -n "$ZOLTAN_ARCH_PATH" ] && [ -n "$WM_MPLIB" ] && [ -d "$FOAM_LIBBIN/$FOAM_MPI" ]
+    then
+        wmakeMpiLib $ZOLTAN_VERSION zoltanDecomp
+    else
+        echo "    skipping zoltanDecomp"
+    fi
+else
+    echo
+    echo "    Error: no config.sh/zoltan settings"
+    echo
+fi
+
+
 # Try and build metisDecomp (has own logic)
 (cd metisDecomp && ./Allwmake $targetType)
 
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
index 4023e2a736..4ecbbb0761 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -108,6 +108,13 @@ Foam::decompositionMethod::NewDecomposer
 
     Info<< "Selecting decomposer " << methodType << endl;
 
+    libs.open
+    (
+        decompositionDict,
+        "libs",
+        decomposerConstructorTablePtr_
+    );
+
     decomposerConstructorTable::iterator cstrIter =
         decomposerConstructorTablePtr_->find(methodType);
 
@@ -141,6 +148,13 @@ Foam::decompositionMethod::NewDistributor
 
     Info<< "Selecting distributor " << methodType << endl;
 
+    libs.open
+    (
+        distributionDict,
+        "libs",
+        distributorConstructorTablePtr_
+    );
+
     distributorConstructorTable::iterator cstrIter =
         distributorConstructorTablePtr_->find(methodType);
 
@@ -288,7 +302,7 @@ void Foam::decompositionMethod::calcCellCells
     // Get agglomerate owner on other side of coupled faces
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    labelList globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
+    labelList globalNeighbour(mesh.nFaces() - mesh.nInternalFaces());
 
     forAll(patches, patchi)
     {
@@ -338,7 +352,7 @@ void Foam::decompositionMethod::calcCellCells
         if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
         {
             label facei = pp.start();
-            label bFacei = pp.start()-mesh.nInternalFaces();
+            label bFacei = pp.start() - mesh.nInternalFaces();
 
             forAll(pp, i)
             {
@@ -389,7 +403,7 @@ void Foam::decompositionMethod::calcCellCells
         if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
         {
             label facei = pp.start();
-            label bFacei = pp.start()-mesh.nInternalFaces();
+            label bFacei = pp.start() - mesh.nInternalFaces();
 
             forAll(pp, i)
             {
@@ -479,7 +493,7 @@ void Foam::decompositionMethod::calcCellCells
     // Get agglomerate owner on other side of coupled faces
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    labelList globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
+    labelList globalNeighbour(mesh.nFaces() - mesh.nInternalFaces());
 
     forAll(patches, patchi)
     {
diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
index 53f1e9cd6f..e4c0c80cf8 100644
--- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
+++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
@@ -796,6 +796,7 @@ Foam::labelList Foam::ptscotchDecomp::decompose
     {
         decomp[i] = finalDecomp[i];
     }
+
     return decomp;
 }
 
diff --git a/src/parallel/decompose/zoltanDecomp/Make/files b/src/parallel/decompose/zoltanDecomp/Make/files
new file mode 100644
index 0000000000..761e86cd80
--- /dev/null
+++ b/src/parallel/decompose/zoltanDecomp/Make/files
@@ -0,0 +1,3 @@
+zoltanDecomp.C
+
+LIB = $(FOAM_LIBBIN)/$(FOAM_MPI)/libzoltanDecomp
diff --git a/src/parallel/decompose/zoltanDecomp/Make/options b/src/parallel/decompose/zoltanDecomp/Make/options
new file mode 100644
index 0000000000..a9e798f03d
--- /dev/null
+++ b/src/parallel/decompose/zoltanDecomp/Make/options
@@ -0,0 +1,10 @@
+-include $(GENERAL_RULES)/mplibType
+
+EXE_INC = \
+    $(PFLAGS) $(PINC) \
+    -I$(FOAM_SRC)/Pstream/mpi/lnInclude \
+    -I$(ZOLTAN_ARCH_PATH)/include/ \
+    -I../decompositionMethods/lnInclude
+
+LIB_LIBS = \
+    -L$(FOAM_EXT_LIBBIN)/$(FOAM_MPI) -lzoltan
diff --git a/src/parallel/decompose/zoltanDecomp/zoltanDecomp.C b/src/parallel/decompose/zoltanDecomp/zoltanDecomp.C
new file mode 100644
index 0000000000..aa8bda2b64
--- /dev/null
+++ b/src/parallel/decompose/zoltanDecomp/zoltanDecomp.C
@@ -0,0 +1,537 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "zoltanDecomp.H"
+#include "globalIndex.H"
+#include "PstreamGlobals.H"
+#include "addToRunTimeSelectionTable.H"
+
+#include "sigFpe.H"
+#ifdef LINUX_GNUC
+    #ifndef __USE_GNU
+        #define __USE_GNU
+    #endif
+    #include <fenv.h>
+#endif
+
+#include "zoltan.h"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(zoltanDecomp, 0);
+
+    addToRunTimeSelectionTable
+    (
+        decompositionMethod,
+        zoltanDecomp,
+        distributor
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+static int get_number_of_vertices(void *data, int *ierr)
+{
+    const Foam::pointField& points =
+        *static_cast<const Foam::pointField*>(data);
+
+    *ierr = ZOLTAN_OK;
+
+    return points.size();
+}
+
+
+static void get_vertex_list
+(
+    void* data,
+    int nGID,
+    int nLID,
+    ZOLTAN_ID_PTR globalIDs,
+    ZOLTAN_ID_PTR localIDs,
+    int wgt_dim,
+    float* obj_wgts,
+    int* ierr
+)
+{
+    const Foam::Tuple2<const Foam::pointField&, const Foam::scalarField&>&
+        vertexData =
+        *static_cast
+        <
+            const Foam::Tuple2
+            <
+                const Foam::pointField&,
+                const Foam::scalarField&
+            >*
+        >(data);
+
+    const Foam::pointField& points = vertexData.first();
+    const Foam::scalarField& weights = vertexData.second();
+
+    Foam::globalIndex globalMap(points.size());
+
+    for (Foam::label i=0; i<points.size(); i++)
+    {
+        localIDs[i] = i;
+        globalIDs[i] = globalMap.toGlobal(i);
+        obj_wgts[i] = weights[i];
+    }
+
+    *ierr = ZOLTAN_OK;
+}
+
+
+static int get_mesh_dim(void* data, int* ierr)
+{
+    return 3;
+}
+
+
+static void get_geom_list
+(
+    void* data,
+    int nGID,
+    int nLID,
+    int nPoints,
+    ZOLTAN_ID_PTR globalIDs,
+    ZOLTAN_ID_PTR localIDs,
+    int nDim,
+    double* vertices,
+    int* ierr
+)
+{
+    const Foam::pointField& points =
+        *static_cast<const Foam::pointField*>(data);
+
+    if
+    (
+        (nGID != 1)
+     || (nLID != 1)
+     || (nPoints != points.size())
+     || (nDim != 3)
+    )
+    {
+        *ierr = ZOLTAN_FATAL;
+        return;
+    }
+
+    double* p = vertices;
+
+    for (Foam::label celli=0; celli<nPoints; celli++)
+    {
+        const Foam::point& pt = points[celli];
+
+        for (Foam::direction cmpt=0; cmpt<Foam::vector::nComponents; cmpt++)
+        {
+            *p++ = pt[cmpt];
+        }
+    }
+
+    *ierr = ZOLTAN_OK;
+}
+
+
+static void get_num_edges_list
+(
+    void* data,
+    int nGID,
+    int nLID,
+    int nPoints,
+    ZOLTAN_ID_PTR globalIDs,
+    ZOLTAN_ID_PTR localIDs,
+    int* numEdges,
+    int* ierr
+)
+{
+    const Foam::CompactListList<Foam::label>& adjacency =
+        *static_cast<const Foam::CompactListList<Foam::label>*>(data);
+
+    Foam::labelList sizes(adjacency.sizes());
+
+    if ((nGID != 1) || (nLID != 1) || (nPoints != sizes.size()))
+    {
+        *ierr = ZOLTAN_FATAL;
+        return;
+    }
+
+    for (Foam::label i=0; i<nPoints; i++)
+    {
+        numEdges[i] = sizes[i];
+    }
+
+    *ierr = ZOLTAN_OK;
+}
+
+
+static void get_edge_list
+(
+    void* data,
+    int nGID,
+    int nLID,
+    int nPoints,
+    ZOLTAN_ID_PTR globalIDs,
+    ZOLTAN_ID_PTR localIDs,
+    int* num_edges,
+    ZOLTAN_ID_PTR nborGID,
+    int* nborProc,
+    int wgt_dim,
+    float* ewgts,
+    int* ierr
+)
+{
+    const Foam::CompactListList<Foam::label>& adjacency =
+        *static_cast<const Foam::CompactListList<Foam::label>*>(data);
+
+    const Foam::labelList& m(adjacency.m());
+
+    Foam::globalIndex globalMap(nPoints);
+
+    if
+    (
+        (nGID != 1)
+     || (nLID != 1)
+     || (wgt_dim != 0)
+    )
+    {
+        *ierr = ZOLTAN_FATAL;
+        return;
+    }
+
+    forAll(m, i)
+    {
+        nborGID[i] = m[i];
+        nborProc[i] = globalMap.whichProcID(m[i]);
+    }
+
+    *ierr = ZOLTAN_OK;
+}
+
+
+Foam::label Foam::zoltanDecomp::decompose
+(
+    const CompactListList<label>& adjacency,
+    const pointField& points,
+    const scalarField& pWeights,
+    List<label>& finalDecomp
+) const
+{
+    stringList args(1);
+    args[0] = "zoltanDecomp";
+
+    int argc = args.size();
+    char* argv[argc];
+    for (label i = 0; i < argc; i++)
+    {
+        argv[i] = strdup(args[i].c_str());
+    }
+
+    float ver;
+    int rc = Zoltan_Initialize(argc, argv, &ver);
+
+    if (rc != ZOLTAN_OK)
+    {
+        FatalErrorInFunction
+            << "Failed initialising Zoltan" << exit(FatalError);
+    }
+
+    struct Zoltan_Struct *zz = Zoltan_Create(PstreamGlobals::MPI_COMM_FOAM);
+
+    // Set internal parameters
+    Zoltan_Set_Param(zz, "return_lists", "export");
+    Zoltan_Set_Param(zz, "obj_weight_dim", "1");
+    Zoltan_Set_Param(zz, "edge_weight_dim", "0");
+
+    // General default paramaters
+    Zoltan_Set_Param(zz, "debug_level", "0");
+    Zoltan_Set_Param(zz, "imbalance_tol", "1.05");
+
+    // Set default method parameters
+    Zoltan_Set_Param(zz, "lb_method", "graph");
+    Zoltan_Set_Param(zz, "lb_approach", "repartition");
+
+    if (decompositionDict_.found("zoltanCoeffs"))
+    {
+        const dictionary& coeffsDict_ =
+            decompositionDict_.subDict("zoltanCoeffs");
+
+        forAllConstIter(IDLList<entry>, coeffsDict_, iter)
+        {
+            if (!iter().isDict())
+            {
+                const word& key = iter().keyword();
+                const word value(iter().stream());
+
+                if (debug)
+                {
+                    Info<< typeName << " : setting parameter " << key
+                        << " to " << value << endl;
+                }
+
+                Zoltan_Set_Param(zz, key.c_str(), value.c_str());
+            }
+        }
+    }
+
+    void* pointsPtr = &const_cast<pointField&>(points);
+    Zoltan_Set_Num_Obj_Fn(zz, get_number_of_vertices, pointsPtr);
+
+    Tuple2<const pointField&, const scalarField&> vertexData(points, pWeights);
+    Zoltan_Set_Obj_List_Fn(zz, get_vertex_list, &vertexData);
+
+    // Callbacks for geometry
+    Zoltan_Set_Num_Geom_Fn(zz, get_mesh_dim, pointsPtr);
+    Zoltan_Set_Geom_Multi_Fn(zz, get_geom_list, pointsPtr);
+
+    // Callbacks for connectivity
+    void* adjacencyPtr = &const_cast<CompactListList<label>&>(adjacency);
+    Zoltan_Set_Num_Edges_Multi_Fn(zz, get_num_edges_list, adjacencyPtr);
+    Zoltan_Set_Edge_List_Multi_Fn(zz, get_edge_list, adjacencyPtr);
+
+    int changed, num_imp, num_exp, *imp_procs, *exp_procs;
+    int *imp_to_part, *exp_to_part;
+    int num_gid_entries, num_lid_entries;
+    ZOLTAN_ID_PTR imp_global_ids, exp_global_ids;
+    ZOLTAN_ID_PTR imp_local_ids, exp_local_ids;
+
+    // Switch off FPU error trapping to work around divide-by-zero bug
+    // in the Zoltan graph and hypergraph methods
+    #ifdef FE_NOMASK_ENV
+    int oldExcepts = fedisableexcept
+    (
+        FE_DIVBYZERO
+      | FE_INVALID
+      | FE_OVERFLOW
+    );
+    #endif
+
+    // Perform load balancing
+    Zoltan_LB_Partition
+    (
+        zz,
+        &changed,
+        &num_gid_entries,
+        &num_lid_entries,
+        &num_imp,
+        &imp_global_ids,
+        &imp_local_ids,
+        &imp_procs,
+        &imp_to_part,
+        &num_exp,
+        &exp_global_ids,
+        &exp_local_ids,
+        &exp_procs,
+        &exp_to_part
+    );
+
+    // Reinstate FPU error trapping
+    #ifdef FE_NOMASK_ENV
+    feenableexcept(oldExcepts);
+    #endif
+
+    if (debug && changed)
+    {
+        Pout << "Zoltan number to move " << changed << " " << num_exp << endl;
+    }
+
+    finalDecomp.setSize(points.size());
+    finalDecomp = Pstream::myProcNo();
+
+    if (changed)
+    {
+        for(int i=0; i<num_exp; i++)
+        {
+            finalDecomp[exp_local_ids[i]] = exp_procs[i];
+        }
+    }
+
+    // Free memory allocated for load-balancing results by Zoltan library
+
+    Zoltan_LB_Free_Part
+    (
+        &imp_global_ids,
+        &imp_local_ids,
+        &imp_procs,
+        &imp_to_part
+    );
+
+    Zoltan_LB_Free_Part
+    (
+        &exp_global_ids,
+        &exp_local_ids,
+        &exp_procs,
+        &exp_to_part
+    );
+
+    Zoltan_Destroy(&zz);
+
+    return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::zoltanDecomp::zoltanDecomp(const dictionary& decompositionDict)
+:
+    decompositionMethod(decompositionDict)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::zoltanDecomp::decompose
+(
+    const polyMesh& mesh,
+    const pointField& cellCentres,
+    const scalarField& cellWeights
+)
+{
+    if (cellCentres.size() != mesh.nCells())
+    {
+        FatalErrorInFunction
+            << "Can use this decomposition method only for the whole mesh"
+            << endl
+            << "and supply one coordinate (cellCentre) for every cell." << endl
+            << "The number of coordinates " << cellCentres.size() << endl
+            << "The number of cells in the mesh " << mesh.nCells()
+            << exit(FatalError);
+    }
+
+
+    // Convert mesh.cellCells() into compressed format
+    CompactListList<label> cellCells;
+    calcCellCells
+    (
+        mesh,
+        identity(mesh.nCells()),
+        mesh.nCells(),
+        true,
+        cellCells
+    );
+
+    // Decompose using default weights
+    List<label> finalDecomp;
+    decompose
+    (
+        cellCells,
+        cellCentres,
+        cellWeights,
+        finalDecomp
+    );
+
+    // Copy back to labelList
+    labelList decomp(cellCentres.size());
+    forAll(decomp, i)
+    {
+        decomp[i] = finalDecomp[i];
+    }
+    return decomp;
+}
+
+
+Foam::labelList Foam::zoltanDecomp::decompose
+(
+    const polyMesh& mesh,
+    const labelList& agglom,
+    const pointField& agglomPoints,
+    const scalarField& pointWeights
+)
+{
+    if (agglom.size() != mesh.nCells())
+    {
+        FatalErrorInFunction
+            << "Size of cell-to-coarse map " << agglom.size()
+            << " differs from number of cells in mesh " << mesh.nCells()
+            << exit(FatalError);
+    }
+
+    // Convert agglom into compressed format
+    CompactListList<label> cellCells;
+    calcCellCells
+    (
+        mesh,
+        agglom,
+        agglomPoints.size(),
+        true,
+        cellCells
+    );
+
+    // Decompose using weights
+    List<label> finalDecomp;
+    decompose
+    (
+        cellCells,
+        agglomPoints,
+        pointWeights,
+        finalDecomp
+    );
+
+    // Rework back into decomposition for original mesh
+    labelList fineDistribution(agglom.size());
+
+    forAll(fineDistribution, i)
+    {
+        fineDistribution[i] = finalDecomp[agglom[i]];
+    }
+
+    return fineDistribution;
+}
+
+
+Foam::labelList Foam::zoltanDecomp::decompose
+(
+    const labelListList& globalPointPoints,
+    const pointField& points,
+    const scalarField& pWeights
+)
+{
+    if (points.size() != globalPointPoints.size())
+    {
+        FatalErrorInFunction
+            << "Inconsistent number of cells (" << globalPointPoints.size()
+            << ") and number of cell centres (" << points.size()
+            << ")." << exit(FatalError);
+    }
+
+
+    // Convert globalPointPoints into compressed format
+    CompactListList<label> pointPoints(globalPointPoints);
+
+    // Decompose using weights
+    List<label> finalDecomp;
+    decompose
+    (
+        pointPoints,
+        points,
+        pWeights,
+        finalDecomp
+    );
+
+    return finalDecomp;
+}
+
+
+// ************************************************************************* //
diff --git a/src/parallel/decompose/zoltanDecomp/zoltanDecomp.H b/src/parallel/decompose/zoltanDecomp/zoltanDecomp.H
new file mode 100644
index 0000000000..ef2b5216bc
--- /dev/null
+++ b/src/parallel/decompose/zoltanDecomp/zoltanDecomp.H
@@ -0,0 +1,186 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::zoltanDecomp
+
+Description
+    Zoltan redistribution in parallel
+
+    Note: Zoltan methods do not support serial operation.
+
+    Parameters
+    - lb_method : The load-balancing algorithm
+      - block : block partitioning
+      - random : random partitioning
+      - rcb : recursive coordinate bisection
+      - rib : ecursive inertial bisection
+      - hsfc : Hilbert space-filling curve partitioning
+      - reftree : refinement tree based partitioning
+      - graph : choose from collection of methods for graphs
+      - hypergraph : choose from a collection of methods for hypergraphs
+
+    - lb_approach The desired load balancing approach. Only lb_method =
+      hypergraph or graph uses the lb_approach parameter. Valid values are
+
+      - partition : Partition without reference to the current distribution,
+        recommended for static load balancing.
+
+      - repartition : Partition starting from the current data distribution
+        to keep data migration low, recommended for dynamic load balancing.
+
+      - refine : Quickly improve the current data distribution
+
+      Default values
+      - debug_level     0
+      - imbalance_tol   1.05
+      - lb_method       graph
+      - lb_approach     repartition
+
+Usage
+    To select the Zoltan graph repartition method add the following entries to
+    decomposeParDict:
+
+        distributor     zoltan;
+        libs            ("libzoltanRenumber.so");
+
+    The Zoltan lb_method and lb_approach can be changed by adding the
+    corresponding entries to the optional zoltanCeoffs sub-dictionary, e.g.:
+
+    zoltanCoeffs
+    {
+        lb_method       hypergraph;
+        lb_approach     partition;
+    }
+
+SourceFiles
+    zoltanDecomp.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef zoltanDecomp_H
+#define zoltanDecomp_H
+
+#include "decompositionMethod.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class zoltanDecomp Declaration
+\*---------------------------------------------------------------------------*/
+
+class zoltanDecomp
+:
+    public decompositionMethod
+{
+    // Private Member Functions
+
+        //- Low level decompose
+        label decompose
+        (
+            const CompactListList<label>& adjacency,
+            const pointField& points,
+            const scalarField& cWeights,
+            List<label>& finalDecomp
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("zoltan");
+
+
+    // Constructors
+
+        //- Construct given the decomposition dictionary and mesh
+        zoltanDecomp(const dictionary& decompositionDict);
+
+        //- Disallow default bitwise copy construction
+        zoltanDecomp(const zoltanDecomp&) = delete;
+
+
+    //- Destructor
+    virtual ~zoltanDecomp()
+    {}
+
+
+    // Member Functions
+
+        //- Inherit decompose from decompositionMethod
+        using decompositionMethod::decompose;
+
+        //- Return for every coordinate the wanted processor number. Use the
+        //  mesh connectivity (if needed). See note on weights in scotchDecomp.H
+        virtual labelList decompose
+        (
+            const polyMesh& mesh,
+            const pointField& points,
+            const scalarField& pointWeights
+        );
+
+        //- Return for every coordinate the wanted processor number. Gets
+        //  passed agglomeration map (from fine to coarse cells) and coarse cell
+        //  location. Can be overridden by decomposers that provide this
+        //  functionality natively. See note on weights in scotchDecomp.H
+        virtual labelList decompose
+        (
+            const polyMesh& mesh,
+            const labelList& agglom,
+            const pointField& regionPoints,
+            const scalarField& regionWeights
+        );
+
+        //- Return for every coordinate the wanted processor number. Explicitly
+        //  provided mesh connectivity.
+        //  The connectivity is equal to mesh.cellCells() except for
+        //  - in parallel the cell numbers are global cell numbers (starting
+        //    from 0 at processor0 and then incrementing all through the
+        //    processors)
+        //  - the connections are across coupled patches
+        //  See note on weights in scotchDecomp.H
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc,
+            const scalarField& cWeights
+        );
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const zoltanDecomp&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/renumber/zoltanRenumber/zoltanRenumber.H b/src/renumber/zoltanRenumber/zoltanRenumber.H
index ef28b19942..b75c3e3544 100644
--- a/src/renumber/zoltanRenumber/zoltanRenumber.H
+++ b/src/renumber/zoltanRenumber/zoltanRenumber.H
@@ -30,7 +30,6 @@ Description
     Note: Zoltan methods do not support serial operation.
 
 Usage
-
     To select the Zoltan Hilbert space filling curves renumber method add the
     following entries to renumberMeshDict:
 
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/decomposeParDict b/tutorials/incompressible/simpleFoam/motorBike/system/decomposeParDict
index 6578c253d4..ea52eba579 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/decomposeParDict
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/decomposeParDict
@@ -18,6 +18,8 @@ numberOfSubdomains  6;
 
 decomposer      hierarchical;
 distributor     ptscotch;
+// distributor     zoltan;
+// libs            ("libzoltanDecomp.so");
 
 hierarchicalCoeffs
 {
@@ -25,5 +27,11 @@ hierarchicalCoeffs
     order           xyz;
 }
 
+zoltanCoeffs
+{
+    lb_method       hypergraph;
+    lb_approach     partition;
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/RAS/floatingObject/system/decomposeParDict b/tutorials/multiphase/interFoam/RAS/floatingObject/system/decomposeParDict
index 00e1de4cab..a828a111e8 100644
--- a/tutorials/multiphase/interFoam/RAS/floatingObject/system/decomposeParDict
+++ b/tutorials/multiphase/interFoam/RAS/floatingObject/system/decomposeParDict
@@ -18,6 +18,8 @@ numberOfSubdomains 8;
 
 decomposer      hierarchical;
 distributor     hierarchical;
+// distributor     zoltan;
+// libs            ("libzoltanDecomp.so");
 
 simpleCoeffs
 {