diff --git a/applications/test/Distribution/DistributionTest.C b/applications/test/Distribution/DistributionTest.C
index e59cb1af1e..3d181b22fb 100644
--- a/applications/test/Distribution/DistributionTest.C
+++ b/applications/test/Distribution/DistributionTest.C
@@ -45,7 +45,8 @@ Description
 #include "Distribution.H"
 #include "Random.H"
 #include "dimensionedTypes.H"
-
+#include "argList.H"
+#include "PstreamReduceOps.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,6 +54,8 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
+    #   include "setRootCase.H"
+
     Random R(918273);
 
     {
@@ -76,7 +79,67 @@ int main(int argc, char *argv[])
             << "Median " << dS.median()
             << endl;
 
-        dS.write("Distribution_scalar_test", dS.normalised());
+        dS.write("Distribution_scalar_test_1");
+
+        Distribution<scalar> dS2(scalar(1e-2));
+
+        Info<< nl << "Distribution<scalar>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from GaussNormal distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dS2.add(1.5*R.GaussNormal() -6.0);
+        }
+
+        Info<< "Mean " << dS2.mean() << nl
+            << "Median " << dS2.median()
+            << endl;
+
+        dS2.write("Distribution_scalar_test_2");
+
+        Info<< nl << "Adding previous two Distribution<scalar>" << endl;
+
+        dS = dS + dS2;
+
+        dS.write("Distribution_scalar_test_1+2");
+    }
+
+    if (Pstream::parRun())
+    {
+        // scalar in parallel
+        label randomDistributionTestSize = 100000000;
+
+        Distribution<scalar> dS(scalar(1e-1));
+
+        Pout<< "Distribution<scalar>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from uniform distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dS.add(R.scalar01() + 10*Pstream::myProcNo());
+        }
+
+        Pout<< "Mean " << dS.mean() << nl
+            << "Median " << dS.median()
+            << endl;
+
+        reduce(dS, sumOp< Distribution<scalar> >());
+
+        if (Pstream::master())
+        {
+            Info<< "Reducing parallel Distribution<scalar>" << nl
+                << "Mean " << dS.mean() << nl
+                << "Median " << dS.median()
+                << endl;
+
+            dS.write("Distribution_scalar_test_parallel_reduced");
+        }
     }
 
     {
@@ -114,40 +177,40 @@ int main(int argc, char *argv[])
             << "Median " << dV.median()
             << endl;
 
-        dV.write("Distribution_vector_test", dV.normalised());
+        dV.write("Distribution_vector_test");
     }
 
-    {
-        // labelVector
-        Distribution<labelVector> dLV(labelVector::one*10);
+    // {
+    //     // labelVector
+    //     Distribution<labelVector> dLV(labelVector::one*10);
 
-        label randomDistributionTestSize = 2000000;
+    //     label randomDistributionTestSize = 2000000;
 
-        Info<< nl << "Distribution<labelVector>" << nl
-            << "Sampling "
-            << randomDistributionTestSize
-            << " times from uniform distribution."
-            << endl;
+    //     Info<< nl << "Distribution<labelVector>" << nl
+    //         << "Sampling "
+    //         << randomDistributionTestSize
+    //         << " times from uniform distribution."
+    //         << endl;
 
-        for (label i = 0; i < randomDistributionTestSize; i++)
-        {
-            dLV.add
-            (
-                labelVector
-                (
-                    R.integer(-1000, 1000),
-                    R.integer(-5000, 5000),
-                    R.integer(-2000, 7000)
-                )
-            );
-        }
+    //     for (label i = 0; i < randomDistributionTestSize; i++)
+    //     {
+    //         dLV.add
+    //         (
+    //             labelVector
+    //             (
+    //                 R.integer(-1000, 1000),
+    //                 R.integer(-5000, 5000),
+    //                 R.integer(-2000, 7000)
+    //             )
+    //         );
+    //     }
 
-        Info<< "Mean " << dLV.mean() << nl
-            << "Median " << dLV.median()
-            << endl;
+    //     Info<< "Mean " << dLV.mean() << nl
+    //         << "Median " << dLV.median()
+    //         << endl;
 
-        dLV.write("Distribution_labelVector_test", dLV.normalised());
-    }
+    //     dLV.write("Distribution_labelVector_test");
+    // }
 
     {
         // tensor
@@ -170,7 +233,7 @@ int main(int argc, char *argv[])
             << "Median " << dT.median()
             << endl;
 
-        dT.write("Distribution_tensor_test", dT.normalised());
+        dT.write("Distribution_tensor_test");
     }
 
     {
@@ -194,7 +257,7 @@ int main(int argc, char *argv[])
             << "Median " << dSyT.median()
             << endl;
 
-        dSyT.write("Distribution_symmTensor_test", dSyT.normalised());
+        dSyT.write("Distribution_symmTensor_test");
     }
 
     {
@@ -218,7 +281,7 @@ int main(int argc, char *argv[])
             << "Median " << dSpT.median()
             << endl;
 
-        dSpT.write("Distribution_sphericalTensor_test", dSpT.normalised());
+        dSpT.write("Distribution_sphericalTensor_test");
     }
 
     Info<< nl << "End" << nl << endl;
diff --git a/src/OpenFOAM/containers/Lists/Distribution/Distribution.C b/src/OpenFOAM/containers/Lists/Distribution/Distribution.C
index 263853428d..d242feb11a 100644
--- a/src/OpenFOAM/containers/Lists/Distribution/Distribution.C
+++ b/src/OpenFOAM/containers/Lists/Distribution/Distribution.C
@@ -28,7 +28,73 @@ License
 #include "OFstream.H"
 #include "ListOps.H"
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Distribution<Type>::Distribution()
+:
+    List< List<scalar> >(pTraits<Type>::nComponents),
+    binWidth_(pTraits<Type>::one),
+    listStarts_(pTraits<Type>::nComponents, 0)
+{}
+
+
+template<class Type>
+Foam::Distribution<Type>::Distribution(const Type& binWidth)
+:
+    List< List<scalar> >(pTraits<Type>::nComponents),
+    binWidth_(binWidth),
+    listStarts_(pTraits<Type>::nComponents, 0)
+{}
+
+
+template<class Type>
+Foam::Distribution<Type>::Distribution(const Distribution<Type>& d)
+:
+    List< List<scalar> >(static_cast< const List< List<scalar> >& >(d)),
+    binWidth_(d.binWidth()),
+    listStarts_(d.listStarts())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Distribution<Type>::~Distribution()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::scalar Foam::Distribution<Type>::totalWeight(direction cmpt) const
+{
+    const List<scalar>& cmptDistribution = (*this)[cmpt];
+
+    scalar sumOfWeights = 0.0;
+
+    forAll(cmptDistribution, i)
+    {
+        sumOfWeights += cmptDistribution[i];
+    }
+
+    return sumOfWeights;
+}
+
+
+template<class Type>
+Foam::List<Foam::label> Foam::Distribution<Type>::keys(direction cmpt) const
+{
+    List<label> keys = identity((*this)[cmpt].size());
+
+    forAll(keys, k)
+    {
+        keys[k] += listStarts_[cmpt];
+    }
+
+    return keys;
+}
+
 
 template<class Type>
 Foam::label Foam::Distribution<Type>::index
@@ -128,73 +194,6 @@ Foam::Pair<Foam::label> Foam::Distribution<Type>::validLimits
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::Distribution<Type>::Distribution()
-:
-    List< List<scalar> >(pTraits<Type>::nComponents),
-    binWidth_(pTraits<Type>::one),
-    listStarts_(pTraits<Type>::nComponents, 0)
-{}
-
-
-template<class Type>
-Foam::Distribution<Type>::Distribution(const Type& binWidth)
-:
-    List< List<scalar> >(pTraits<Type>::nComponents),
-    binWidth_(binWidth),
-    listStarts_(pTraits<Type>::nComponents, 0)
-{}
-
-
-template<class Type>
-Foam::Distribution<Type>::Distribution(const Distribution<Type>& d)
-:
-    List< List<scalar> >(static_cast< const List< List<scalar> >& >(d)),
-    binWidth_(d.binWidth()),
-    listStarts_(d.listStarts())
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::Distribution<Type>::~Distribution()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::scalar Foam::Distribution<Type>::totalWeight(direction cmpt) const
-{
-    const List<scalar>& cmptDistribution = (*this)[cmpt];
-
-    scalar sumOfWeights = 0.0;
-
-    forAll(cmptDistribution, i)
-    {
-        sumOfWeights += cmptDistribution[i];
-    }
-
-    return sumOfWeights;
-}
-
-template<class Type>
-Foam::List<Foam::label> Foam::Distribution<Type>::keys(direction cmpt) const
-{
-    List<label> keys = identity((*this)[cmpt].size());
-
-    forAll(keys, k)
-    {
-        keys[k] += listStarts_[cmpt];
-    }
-
-    return keys;
-}
-
-
 template<class Type>
 Type Foam::Distribution<Type>::mean() const
 {
@@ -225,7 +224,7 @@ Type Foam::Distribution<Type>::mean() const
 
 
 template<class Type>
-Type Foam::Distribution<Type>::median()
+Type Foam::Distribution<Type>::median() const
 {
     Type medianValue(pTraits<Type>::zero);
 
@@ -247,18 +246,17 @@ Type Foam::Distribution<Type>::median()
              && normDist[0].second()*component(binWidth_, cmpt) > 0.5
             )
             {
-                scalar xk = normDist[1].first();
+                scalar xk =
+                    normDist[0].first()
+                  + 0.5*component(binWidth_, cmpt);
 
-                scalar xkm1 = normDist[0].first();
+                scalar xkm1 =
+                    normDist[0].first()
+                  - 0.5*component(binWidth_, cmpt);
 
-                scalar Sk =
-                    (normDist[0].second() + normDist[1].second())
-                   *component(binWidth_, cmpt);
+                scalar Sk = (normDist[0].second())*component(binWidth_, cmpt);
 
-                scalar Skm1 = normDist[0].second()*component(binWidth_, cmpt);
-
-                setComponent(medianValue, cmpt) =
-                    (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
+                setComponent(medianValue, cmpt) = 0.5*(xk - xkm1)/(Sk) + xkm1;
             }
             else
             {
@@ -266,7 +264,7 @@ Type Foam::Distribution<Type>::median()
 
                 scalar cumulative = 0.0;
 
-                forAll(normDist,nD)
+                forAll(normDist, nD)
                 {
                     if
                     (
@@ -275,9 +273,13 @@ Type Foam::Distribution<Type>::median()
                       > 0.5
                     )
                     {
-                        scalar xk = normDist[nD].first();
+                        scalar xk =
+                            normDist[nD].first()
+                          + 0.5*component(binWidth_, cmpt);
 
-                        scalar xkm1 = normDist[previousNonZeroIndex].first();
+                        scalar xkm1 =
+                            normDist[previousNonZeroIndex].first()
+                          + 0.5*component(binWidth_, cmpt);
 
                         scalar Sk =
                             cumulative
@@ -314,7 +316,6 @@ void Foam::Distribution<Type>::add
     const Type& weight
 )
 {
-
     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
         List<scalar>& cmptDistribution = (*this)[cmpt];
@@ -332,13 +333,13 @@ void Foam::Distribution<Type>::add
 
 template<class Type>
 Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
-Distribution<Type>::normalised()
+Distribution<Type>::normalised() const
 {
     List< List < Pair<scalar> > > normDistribution(pTraits<Type>::nComponents);
 
     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
-        List<scalar>& cmptDistribution = (*this)[cmpt];
+        const List<scalar>& cmptDistribution = (*this)[cmpt];
 
         if (cmptDistribution.empty())
         {
@@ -380,13 +381,13 @@ Distribution<Type>::normalised()
 
 template<class Type>
 Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
-Distribution<Type>::raw()
+Distribution<Type>::raw() const
 {
     List< List < Pair<scalar> > > rawDistribution(pTraits<Type>::nComponents);
 
     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
-        List<scalar>& cmptDistribution = (*this)[cmpt];
+        const List<scalar>& cmptDistribution = (*this)[cmpt];
 
         if (cmptDistribution.empty())
         {
@@ -433,37 +434,28 @@ void Foam::Distribution<Type>::clear()
 
 
 template<class Type>
-void Foam::Distribution<Type>::write
-(
-    const fileName& filePrefix,
-    const List< List< Pair<scalar> > >& pairs
-) const
+void Foam::Distribution<Type>::write(const fileName& filePrefix) const
 {
-    if (pairs.size() != pTraits<Type>::nComponents)
-    {
-        FatalErrorIn
-        (
-            "Distribution::write"
-            "("
-                "const fileName& filePrefix,"
-                "const List< List< Pair<scalar> > >& pairs"
-            ")"
-        )
-            << "List of pairs (" << pairs.size()
-            << ") is not the same size as the number of components ("
-            << pTraits<Type>::nComponents << ")." << nl
-            << abort(FatalError);
-    }
+    List< List<Pair<scalar> > > rawDistribution = raw();
+
+    List< List < Pair<scalar> > > normDistribution = normalised();
 
     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
-        const List< Pair<scalar> >& cmptPairs = pairs[cmpt];
+        const List< Pair<scalar> >& rawPairs = rawDistribution[cmpt];
+
+        const List< Pair<scalar> >& normPairs = normDistribution[cmpt];
 
         OFstream os(filePrefix + '_' + pTraits<Type>::componentNames[cmpt]);
 
-        forAll(cmptPairs, i)
+        os  << "# key normalised raw" << endl;
+
+        forAll(normPairs, i)
         {
-            os  << cmptPairs[i].first() << ' ' << cmptPairs[i].second() << nl;
+            os  << normPairs[i].first()
+                << ' ' << normPairs[i].second()
+                << ' ' << rawPairs[i].second()
+                << nl;
         }
     }
 }
@@ -491,11 +483,31 @@ void Foam::Distribution<Type>::operator=
     List< List<scalar> >::operator=(rhs);
 
     binWidth_ = rhs.binWidth();
+
+    listStarts_ = rhs.listStarts();
 }
 
 
 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
 
+template <class Type>
+Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    Distribution<Type>& d
+)
+{
+    is  >> static_cast<List< List<scalar> >&>(d)
+        >> d.binWidth_
+        >> d.listStarts_;
+
+    // Check state of Istream
+    is.check("Istream& operator>>(Istream&, Distribution<Type>&)");
+
+    return is;
+}
+
+
 template<class Type>
 Foam::Ostream& Foam::operator<<
 (
@@ -503,18 +515,68 @@ Foam::Ostream& Foam::operator<<
     const Distribution<Type>& d
 )
 {
-    os  << d.binWidth_
-        << static_cast<const List< List<scalar> >& >(d);
+    os  <<  static_cast<const List< List<scalar> >& >(d)
+        << d.binWidth_ << token::SPACE
+        << d.listStarts_;
 
     // Check state of Ostream
-    os.check
-    (
-        "Ostream& operator<<(Ostream&, "
-        "const Distribution&)"
-    );
+    os.check("Ostream& operator<<(Ostream&, " "const Distribution&)");
 
     return os;
 }
 
 
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
+
+template <class Type>
+Foam::Distribution<Type> Foam::operator+
+(
+    const Distribution<Type>& d1,
+    const Distribution<Type>& d2
+)
+{
+    // The coarsest binWidth is the sensible choice
+    Distribution<Type> d(max(d1.binWidth(), d2.binWidth()));
+
+    List< List< List < Pair<scalar> > > > rawDists(2);
+
+    rawDists[0] = d1.raw();
+    rawDists[1] = d2.raw();
+
+    forAll(rawDists, rDI)
+    {
+        for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+        {
+            List<scalar>& cmptDistribution = d[cmpt];
+
+            const List < Pair<scalar> >& cmptRaw = rawDists[rDI][cmpt];
+
+            forAll(cmptRaw, rI)
+            {
+                scalar valueToAdd = cmptRaw[rI].first();
+                scalar cmptWeight = cmptRaw[rI].second();
+
+                label n =
+                label
+                (
+                    component(valueToAdd, cmpt)
+                   /component(d.binWidth(), cmpt)
+                )
+                - label
+                (
+                    neg(component(valueToAdd, cmpt)
+                   /component(d.binWidth(), cmpt))
+                );
+
+                label listIndex = d.index(cmpt, n);
+
+                cmptDistribution[listIndex] += cmptWeight;
+            }
+        }
+    }
+
+    return Distribution<Type>(d);
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/containers/Lists/Distribution/Distribution.H b/src/OpenFOAM/containers/Lists/Distribution/Distribution.H
index 5dbdb444c4..ef6ef25676 100644
--- a/src/OpenFOAM/containers/Lists/Distribution/Distribution.H
+++ b/src/OpenFOAM/containers/Lists/Distribution/Distribution.H
@@ -50,6 +50,9 @@ namespace Foam
 template<class Type>
 class Distribution;
 
+template<class Type>
+Istream& operator>>(Istream&, Distribution<Type>&);
+
 template<class Type>
 Ostream& operator<<(Ostream&, const Distribution<Type>&);
 
@@ -70,15 +73,6 @@ class Distribution
         //- The start bin index of each component
         List<label> listStarts_;
 
-    // Private Member Functions
-
-        //- Return the appropriate List index for the given bin index.
-        //  Resizes the List if required
-        label index(direction cmpt, label n);
-
-        //- Returns the indices of the first and last non-zero entries
-        Pair<label> validLimits(direction cmpt) const;
-
 
 public:
 
@@ -109,12 +103,19 @@ public:
 
         List<label> keys(direction cmpt) const;
 
+        //- Return the appropriate List index for the given bin index.
+        //  Resizes the List if required
+        label index(direction cmpt, label n);
+
+        //- Returns the indices of the first and last non-zero entries
+        Pair<label> validLimits(direction cmpt) const;
+
         Type mean() const;
 
         // From http://mathworld.wolfram.com/StatisticalMedian.html
         // The statistical median is the value of the Distribution
         // variable where the cumulative Distribution = 0.5.
-        Type median();
+        Type median() const;
 
         //- Add a value to the distribution, optionally specifying a weight
         void add
@@ -124,10 +125,10 @@ public:
         );
 
         //- Return the normalised distribution and bins
-        List< List<Pair<scalar> > > normalised();
+        List< List<Pair<scalar> > > normalised() const;
 
         //- Return the distribution of the total bin weights
-        List< List < Pair<scalar> > > raw();
+        List< List < Pair<scalar> > > raw() const;
 
         //- Resets the Distribution by clearing the stored lists.
         //  Leaves the same number of them and the same binWidth.
@@ -144,13 +145,9 @@ public:
 
     // Write
 
-        //- Write the distribution to file.  Produces a separate file
-        //  for each component.
-        void write
-        (
-            const fileName& filePrefix,
-            const List< List<Pair<scalar> > >& pairs
-        ) const;
+        //- Write the distribution to file: key normalised raw.
+        //  Produces a separate file for each component.
+        void write(const fileName& filePrefix) const;
 
 
     // Member Operators
@@ -159,6 +156,12 @@ public:
 
     // IOstream Operators
 
+        friend Istream& operator>> <Type>
+        (
+            Istream&,
+            Distribution<Type>&
+        );
+
         friend Ostream& operator<< <Type>
         (
             Ostream&,
@@ -167,7 +170,15 @@ public:
 };
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
+
+template<class Type>
+Distribution<Type> operator+
+(
+    const Distribution<Type>&,
+    const Distribution<Type>&
+);
+
 
 } // End namespace Foam