diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomerate.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomerate.C
index b0f5e757da..508c0c3f94 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomerate.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomerate.C
@@ -26,177 +26,7 @@ License
 #include "pairGAMGAgglomeration.H"
 #include "lduAddressing.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
-(
-    label& nCoarseCells,
-    const lduAddressing& fineMatrixAddressing,
-    const scalarField& faceWeights
-)
-{
-    const label nFineCells = fineMatrixAddressing.size();
-
-    const labelUList& upperAddr = fineMatrixAddressing.upperAddr();
-    const labelUList& lowerAddr = fineMatrixAddressing.lowerAddr();
-
-    // For each cell calculate faces
-    labelList cellFaces(upperAddr.size() + lowerAddr.size());
-    labelList cellFaceOffsets(nFineCells + 1);
-
-    // memory management
-    {
-        labelList nNbrs(nFineCells, 0);
-
-        forAll(upperAddr, facei)
-        {
-            nNbrs[upperAddr[facei]]++;
-        }
-
-        forAll(lowerAddr, facei)
-        {
-            nNbrs[lowerAddr[facei]]++;
-        }
-
-        cellFaceOffsets[0] = 0;
-        forAll(nNbrs, celli)
-        {
-            cellFaceOffsets[celli+1] = cellFaceOffsets[celli] + nNbrs[celli];
-        }
-
-        // reset the whole list to use as counter
-        nNbrs = 0;
-
-        forAll(upperAddr, facei)
-        {
-            cellFaces
-            [
-                cellFaceOffsets[upperAddr[facei]] + nNbrs[upperAddr[facei]]
-            ] = facei;
-
-            nNbrs[upperAddr[facei]]++;
-        }
-
-        forAll(lowerAddr, facei)
-        {
-            cellFaces
-            [
-                cellFaceOffsets[lowerAddr[facei]] + nNbrs[lowerAddr[facei]]
-            ] = facei;
-
-            nNbrs[lowerAddr[facei]]++;
-        }
-    }
-
-
-    // go through the faces and create clusters
-
-    tmp<labelField> tcoarseCellMap(new labelField(nFineCells, -1));
-    labelField& coarseCellMap = tcoarseCellMap();
-
-    nCoarseCells = 0;
-
-    for (label celli=0; celli<nFineCells; celli++)
-    {
-        if (coarseCellMap[celli] < 0)
-        {
-            label matchFaceNo = -1;
-            scalar maxFaceWeight = -GREAT;
-
-            // check faces to find ungrouped neighbour with largest face weight
-            for
-            (
-                label faceOs=cellFaceOffsets[celli];
-                faceOs<cellFaceOffsets[celli+1];
-                faceOs++
-            )
-            {
-                label facei = cellFaces[faceOs];
-
-                // I don't know whether the current cell is owner or neighbour.
-                // Therefore I'll check both sides
-                if
-                (
-                    coarseCellMap[upperAddr[facei]] < 0
-                 && coarseCellMap[lowerAddr[facei]] < 0
-                 && faceWeights[facei] > maxFaceWeight
-                )
-                {
-                    // Match found. Pick up all the necessary data
-                    matchFaceNo = facei;
-                    maxFaceWeight = faceWeights[facei];
-                }
-            }
-
-            if (matchFaceNo >= 0)
-            {
-                // Make a new group
-                coarseCellMap[upperAddr[matchFaceNo]] = nCoarseCells;
-                coarseCellMap[lowerAddr[matchFaceNo]] = nCoarseCells;
-                nCoarseCells++;
-            }
-            else
-            {
-                // No match. Find the best neighbouring cluster and
-                // put the cell there
-                label clusterMatchFaceNo = -1;
-                scalar clusterMaxFaceCoeff = -GREAT;
-
-                for
-                (
-                    label faceOs=cellFaceOffsets[celli];
-                    faceOs<cellFaceOffsets[celli+1];
-                    faceOs++
-                )
-                {
-                    label facei = cellFaces[faceOs];
-
-                    if (faceWeights[facei] > clusterMaxFaceCoeff)
-                    {
-                        clusterMatchFaceNo = facei;
-                        clusterMaxFaceCoeff = faceWeights[facei];
-                    }
-                }
-
-                if (clusterMatchFaceNo >= 0)
-                {
-                    // Add the cell to the best cluster
-                    coarseCellMap[celli] = max
-                    (
-                        coarseCellMap[upperAddr[clusterMatchFaceNo]],
-                        coarseCellMap[lowerAddr[clusterMatchFaceNo]]
-                    );
-                }
-            }
-        }
-    }
-
-
-    // Check that all cells are part of clusters,
-    // if not create single-cell "clusters" for each
-    for (label celli=0; celli<nFineCells; celli++)
-    {
-        if (coarseCellMap[celli] < 0)
-        {
-            coarseCellMap[celli] = nCoarseCells;
-            nCoarseCells++;
-        }
-    }
-
-    // Reverse the map ordering to improve the next level of agglomeration
-    // (doesn't always help and is sometimes detrimental)
-    nCoarseCells--;
-
-    forAll(coarseCellMap, celli)
-    {
-        coarseCellMap[celli] = nCoarseCells - coarseCellMap[celli];
-    }
-
-    nCoarseCells++;
-
-    return tcoarseCellMap;
-}
-
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
 void Foam::pairGAMGAgglomeration::agglomerate
 (
@@ -285,4 +115,187 @@ void Foam::pairGAMGAgglomeration::agglomerate
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
+(
+    label& nCoarseCells,
+    const lduAddressing& fineMatrixAddressing,
+    const scalarField& faceWeights
+)
+{
+    const label nFineCells = fineMatrixAddressing.size();
+
+    const labelUList& upperAddr = fineMatrixAddressing.upperAddr();
+    const labelUList& lowerAddr = fineMatrixAddressing.lowerAddr();
+
+    // For each cell calculate faces
+    labelList cellFaces(upperAddr.size() + lowerAddr.size());
+    labelList cellFaceOffsets(nFineCells + 1);
+
+    // memory management
+    {
+        labelList nNbrs(nFineCells, 0);
+
+        forAll(upperAddr, facei)
+        {
+            nNbrs[upperAddr[facei]]++;
+        }
+
+        forAll(lowerAddr, facei)
+        {
+            nNbrs[lowerAddr[facei]]++;
+        }
+
+        cellFaceOffsets[0] = 0;
+        forAll(nNbrs, celli)
+        {
+            cellFaceOffsets[celli+1] = cellFaceOffsets[celli] + nNbrs[celli];
+        }
+
+        // reset the whole list to use as counter
+        nNbrs = 0;
+
+        forAll(upperAddr, facei)
+        {
+            cellFaces
+            [
+                cellFaceOffsets[upperAddr[facei]] + nNbrs[upperAddr[facei]]
+            ] = facei;
+
+            nNbrs[upperAddr[facei]]++;
+        }
+
+        forAll(lowerAddr, facei)
+        {
+            cellFaces
+            [
+                cellFaceOffsets[lowerAddr[facei]] + nNbrs[lowerAddr[facei]]
+            ] = facei;
+
+            nNbrs[lowerAddr[facei]]++;
+        }
+    }
+
+
+    // go through the faces and create clusters
+
+    tmp<labelField> tcoarseCellMap(new labelField(nFineCells, -1));
+    labelField& coarseCellMap = tcoarseCellMap();
+
+    nCoarseCells = 0;
+    label celli;
+
+    for (label cellfi=0; cellfi<nFineCells; cellfi++)
+    {
+        // Change cell ordering depending on direction for this level
+        celli = forward_ ? cellfi : nFineCells - cellfi - 1;
+
+        if (coarseCellMap[celli] < 0)
+        {
+            label matchFaceNo = -1;
+            scalar maxFaceWeight = -GREAT;
+
+            // check faces to find ungrouped neighbour with largest face weight
+            for
+            (
+                label faceOs=cellFaceOffsets[celli];
+                faceOs<cellFaceOffsets[celli+1];
+                faceOs++
+            )
+            {
+                label facei = cellFaces[faceOs];
+
+                // I don't know whether the current cell is owner or neighbour.
+                // Therefore I'll check both sides
+                if
+                (
+                    coarseCellMap[upperAddr[facei]] < 0
+                 && coarseCellMap[lowerAddr[facei]] < 0
+                 && faceWeights[facei] > maxFaceWeight
+                )
+                {
+                    // Match found. Pick up all the necessary data
+                    matchFaceNo = facei;
+                    maxFaceWeight = faceWeights[facei];
+                }
+            }
+
+            if (matchFaceNo >= 0)
+            {
+                // Make a new group
+                coarseCellMap[upperAddr[matchFaceNo]] = nCoarseCells;
+                coarseCellMap[lowerAddr[matchFaceNo]] = nCoarseCells;
+                nCoarseCells++;
+            }
+            else
+            {
+                // No match. Find the best neighbouring cluster and
+                // put the cell there
+                label clusterMatchFaceNo = -1;
+                scalar clusterMaxFaceCoeff = -GREAT;
+
+                for
+                (
+                    label faceOs=cellFaceOffsets[celli];
+                    faceOs<cellFaceOffsets[celli+1];
+                    faceOs++
+                )
+                {
+                    label facei = cellFaces[faceOs];
+
+                    if (faceWeights[facei] > clusterMaxFaceCoeff)
+                    {
+                        clusterMatchFaceNo = facei;
+                        clusterMaxFaceCoeff = faceWeights[facei];
+                    }
+                }
+
+                if (clusterMatchFaceNo >= 0)
+                {
+                    // Add the cell to the best cluster
+                    coarseCellMap[celli] = max
+                    (
+                        coarseCellMap[upperAddr[clusterMatchFaceNo]],
+                        coarseCellMap[lowerAddr[clusterMatchFaceNo]]
+                    );
+                }
+            }
+        }
+    }
+
+    // Check that all cells are part of clusters,
+    // if not create single-cell "clusters" for each
+    for (label cellfi=0; cellfi<nFineCells; cellfi++)
+    {
+        // Change cell ordering depending on direction for this level
+        celli = forward_ ? cellfi : nFineCells - cellfi - 1;
+
+        if (coarseCellMap[celli] < 0)
+        {
+            coarseCellMap[celli] = nCoarseCells;
+            nCoarseCells++;
+        }
+    }
+
+    if (!forward_)
+    {
+        nCoarseCells--;
+
+        forAll(coarseCellMap, celli)
+        {
+            coarseCellMap[celli] = nCoarseCells - coarseCellMap[celli];
+        }
+
+        nCoarseCells++;
+    }
+
+    // Reverse the map ordering for the next level
+    // to improve the next level of agglomeration
+    forward_ = !forward_;
+
+    return tcoarseCellMap;
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.C
index 1c4cbfcd7f..31359e692c 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,7 +29,8 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(pairGAMGAgglomeration, 0);
+    defineTypeNameAndDebug(pairGAMGAgglomeration, 0);
+    bool pairGAMGAgglomeration::forward_(true);
 }
 
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.H
index 59de4b1e4d..02d6c5f1c4 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.H
@@ -56,6 +56,9 @@ class pairGAMGAgglomeration
         //- Number of levels to merge, 1 = don't merge, 2 = merge pairs etc.
         label mergeLevels_;
 
+        //- Direction of cell loop for the current level
+        static bool forward_;
+
 
 protected:
 
@@ -97,7 +100,6 @@ public:
             const lduAddressing& fineMatrixAddressing,
             const scalarField& faceWeights
         );
-
 };