diff --git a/src/fvOptions/sources/derived/meanVelocityForce/patchMeanVelocityForce/patchMeanVelocityForce.C b/src/fvOptions/sources/derived/meanVelocityForce/patchMeanVelocityForce/patchMeanVelocityForce.C
index 9b1acd02b3..cccc4db367 100644
--- a/src/fvOptions/sources/derived/meanVelocityForce/patchMeanVelocityForce/patchMeanVelocityForce.C
+++ b/src/fvOptions/sources/derived/meanVelocityForce/patchMeanVelocityForce/patchMeanVelocityForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,12 +74,52 @@ Foam::scalar Foam::fv::patchMeanVelocityForce::magUbarAve
     const volVectorField& U
 ) const
 {
-    return
-        gSum
+    vector2D sumAmagUsumA
+    (
+        sum
         (
             (flowDir_ & U.boundaryField()[patchi_])
            *mesh_.boundary()[patchi_].magSf()
-        )/gSum(mesh_.boundary()[patchi_].magSf());
+        ),
+        sum(mesh_.boundary()[patchi_].magSf())
+    );
+
+
+    // If the mean velocity force is applied to a cyclic patch
+    // for parallel runs include contributions from processorCyclic patches
+    // generated from the decomposition of the cyclic patch
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    if (Pstream::parRun() && isA<cyclicPolyPatch>(patches[patchi_]))
+    {
+        const keyType processorCyclicPatchNames
+        (
+            string("procBoundary.*to.*through" + patch_)
+        );
+
+        labelList processorCyclicPatches
+        (
+            patches.findIndices(processorCyclicPatchNames)
+        );
+
+        forAll(processorCyclicPatches, pcpi)
+        {
+            const label patchi = processorCyclicPatches[pcpi];
+
+            sumAmagUsumA.x() +=
+                sum
+                (
+                    (flowDir_ & U.boundaryField()[patchi])
+                   *mesh_.boundary()[patchi].magSf()
+                );
+
+            sumAmagUsumA.y() += sum(mesh_.boundary()[patchi].magSf());
+        }
+    }
+
+    mesh_.reduce(sumAmagUsumA, sumOp<vector2D>());
+
+    return sumAmagUsumA.x()/sumAmagUsumA.y();
 }