diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
index 6ce7b3350ddb4e5eac90b3bc8338ece4e6504e78..600d3d9ecdc505d573b6a863de5cfb624b8d77ec 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
@@ -29,6 +29,8 @@ License
 #include "emptyPolyPatch.H"
 #include "coupledPolyPatch.H"
 #include "sampledSurface.H"
+#include "mergePoints.H"
+#include "indirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -222,6 +224,127 @@ void Foam::fieldValues::faceSource::sampledSurfaceFaces(const dictionary& dict)
 }
 
 
+void Foam::fieldValues::faceSource::combineSurfaceGeometry
+(
+    faceList& faces,
+    pointField& points
+) const
+{
+    List<faceList> allFaces(Pstream::nProcs());
+    List<pointField> allPoints(Pstream::nProcs());
+
+    labelList globalFacesIs(faceId_);
+    forAll(globalFacesIs, i)
+    {
+        if (facePatchId_[i] != -1)
+        {
+            label patchI = facePatchId_[i];
+            globalFacesIs[i] += mesh().boundaryMesh()[patchI].start();
+        }
+    }
+
+    // Add local faces and points to the all* lists
+    indirectPrimitivePatch pp
+    (
+        IndirectList<face>(mesh().faces(), globalFacesIs),
+        mesh().points()
+    );
+    allFaces[Pstream::myProcNo()] = pp.localFaces();
+    allPoints[Pstream::myProcNo()] = pp.localPoints();
+
+    Pstream::gatherList(allFaces);
+    Pstream::gatherList(allPoints);
+
+    // Renumber and flatten
+    label nFaces = 0;
+    label nPoints = 0;
+    forAll(allFaces, procI)
+    {
+        nFaces += allFaces[procI].size();
+        nPoints += allPoints[procI].size();
+    }
+
+    faces.setSize(nFaces);
+    points.setSize(nPoints);
+
+    nFaces = 0;
+    nPoints = 0;
+
+    // My own data first
+    {
+        const faceList& fcs = allFaces[Pstream::myProcNo()];
+        forAll(fcs, i)
+        {
+            const face& f = fcs[i];
+            face& newF = faces[nFaces++];
+            newF.setSize(f.size());
+            forAll(f, fp)
+            {
+                newF[fp] = f[fp] + nPoints;
+            }
+        }
+
+        const pointField& pts = allPoints[Pstream::myProcNo()];
+        forAll(pts, i)
+        {
+            points[nPoints++] = pts[i];
+        }
+    }
+
+    // Other proc data follows
+    forAll(allFaces, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            const faceList& fcs = allFaces[procI];
+            forAll(fcs, i)
+            {
+                const face& f = fcs[i];
+                face& newF = faces[nFaces++];
+                newF.setSize(f.size());
+                forAll(f, fp)
+                {
+                    newF[fp] = f[fp] + nPoints;
+                }
+            }
+
+            const pointField& pts = allPoints[procI];
+            forAll(pts, i)
+            {
+                points[nPoints++] = pts[i];
+            }
+        }
+    }
+
+    // Merge
+    labelList oldToNew;
+    pointField newPoints;
+    bool hasMerged = mergePoints
+    (
+        points,
+        SMALL,
+        false,
+        oldToNew,
+        newPoints
+    );
+
+    if (hasMerged)
+    {
+        if (debug)
+        {
+            Pout<< "Merged from " << points.size()
+                << " down to " << newPoints.size() << " points" << endl;
+        }
+
+        points.transfer(newPoints);
+        forAll(faces, i)
+        {
+            inplaceRenumber(oldToNew, faces[i]);
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 void Foam::fieldValues::faceSource::initialise(const dictionary& dict)
@@ -291,6 +414,14 @@ void Foam::fieldValues::faceSource::initialise(const dictionary& dict)
     }
 
     Info<< nl << endl;
+
+    if (valueOutput_)
+    {
+        surfaceWriterPtr_.reset
+        (
+            surfaceWriter::New(dict.lookup("surfaceFormat")).ptr()
+        );
+    }
 }
 
 
@@ -355,6 +486,7 @@ Foam::fieldValues::faceSource::faceSource
 )
 :
     fieldValue(name, obr, dict, loadFromFiles),
+    surfaceWriterPtr_(NULL),
     source_(sourceTypeNames_.read(dict.lookup("source"))),
     operation_(operationTypeNames_.read(dict.lookup("operation"))),
     weightFieldName_("none"),
diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
index fd95788894aa3e9d63ebc083045e422f66d6164c..c4473b352a8209921d1381805fe17c005e554db1 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
@@ -35,9 +35,10 @@ Description
         enabled         true;
         outputControl   outputTime;
         log             true;       // log to screen?
-        valueOutput     true;       // Write values at run-time output times?
-        source          faceZone;   // Type of face source:
-                                    // faceZone,patch,sampledSurface
+        valueOutput     true;       // write the output values
+        surfaceFormat   none;       // output value format (if valueOutput)
+        source          faceZone;   // type of face source:
+                                    // faceZone, patch, sampledSurface
         sourceName      f0;         // faceZone name, see below
         operation       sum;
         weightField     alpha1;     // optional weight field
@@ -97,6 +98,7 @@ SourceFiles
 #include "fieldValue.H"
 #include "surfaceFieldsFwd.H"
 #include "volFieldsFwd.H"
+#include "surfaceWriter.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -166,11 +168,21 @@ private:
         //- Set faces according to sampledSurface
         void sampledSurfaceFaces(const dictionary&);
 
+        //- Combine faces and points from multiple processors 
+        void combineSurfaceGeometry
+        (
+            faceList& faces,
+            pointField& points
+        ) const;
+
 
 protected:
 
     // Protected data
 
+        //- Surface writer
+        autoPtr<surfaceWriter> surfaceWriterPtr_;
+
         //- Source type
         sourceType source_;
 
@@ -184,7 +196,7 @@ protected:
         label nFaces_;
 
 
-        // If operating on mesh faces (faceZone,patch)
+        // If operating on mesh faces (faceZone, patch)
 
             //- Local list of face IDs
             labelList faceId_;
@@ -196,6 +208,7 @@ protected:
             //  (1 use as is, -1 negate)
             labelList faceSign_;
 
+
         // If operating on sampledSurface
 
             //- underlying sampledSurface
diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
index 450ac18d91db2ec9a305de6b83938b6679318e49..cfdcc08d4fb39214a1d88d8ca8d0a712cb4a07d2 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
@@ -254,31 +254,46 @@ bool Foam::fieldValues::faceSource::writeValues(const word& fieldName)
         combineFields(Sf);
         combineFields(weightField);
 
+        // Write raw values on surface if specified
+        if (surfaceWriterPtr_.valid())
+        {
+            faceList faces;
+            pointField points;
+            combineSurfaceGeometry(faces, points);
+
+            fileName outputDir;
+            if (Pstream::parRun())
+            {
+                // Put in undecomposed case (Note: gives problems for
+                // distributed data running)
+                outputDir = obr_.time().path()/".."/name_;
+            }
+            else
+            {
+                outputDir = obr_.time().path()/name_;
+            }
+
+            outputDir = outputDir/"surface"/obr_.time().timeName();
+
+            surfaceWriterPtr_->write
+            (
+                outputDir,
+                word(sourceTypeNames_[source_]) + "_" + sourceName_,
+                points,
+                faces,
+                fieldName,
+                values,
+                false
+            );
+        }
+
         // apply weight field
         values *= weightField;
 
-
         if (Pstream::master())
         {
             Type result = processValues(values, Sf, weightField);
 
-            if (valueOutput_)
-            {
-                IOField<Type>
-                (
-                    IOobject
-                    (
-                        fieldName + "_" + sourceTypeNames_[source_] + "-"
-                            + sourceName_,
-                        obr_.time().timeName(),
-                        obr_,
-                        IOobject::NO_READ,
-                        IOobject::NO_WRITE
-                    ),
-                    values
-                ).write();
-            }
-
             outputFilePtr_()<< tab << result;
 
             if (log_)