diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.C
index 1068ea0096..d5f76f5af2 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.C
@@ -32,20 +32,28 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class CloudType>
-Foam::label Foam::FacePostProcessing<CloudType>::applyToFace
+void Foam::FacePostProcessing<CloudType>::applyToFace
 (
-    const label faceI
+    const label faceIn,
+    label& zoneI,
+    label& faceI
 ) const
 {
-    forAll(fZone_, i)
-    {
-        if (fZone_[i] == faceI)
-        {
-            return i;
-        }
-    }
+    const faceZoneMesh& fzm = this->owner().mesh().faceZones();
 
-    return -1;
+    forAll(faceZoneIDs_, i)
+    {
+        const faceZone& fz = fzm[faceZoneIDs_[i]];
+        forAll(fz, j)
+        {
+            if (fz[j] == faceIn)
+            {
+                zoneI = i;
+                faceI = j;
+                return;
+            }
+        }        
+    }
 }
 
 
@@ -53,6 +61,7 @@ template<class CloudType>
 void Foam::FacePostProcessing<CloudType>::write()
 {
     const fvMesh& mesh = this->owner().mesh();
+    const faceZoneMesh& fzm = this->owner().mesh().faceZones();
     const scalar dt = this->owner().time().deltaTValue();
 
     totalTime_ += dt;
@@ -60,133 +69,159 @@ void Foam::FacePostProcessing<CloudType>::write()
     const scalar alpha = (totalTime_ - dt)/totalTime_;
     const scalar beta = dt/totalTime_;
 
-    massTotal_ += mass_;
-
-    massFlux_ = alpha*massFlux_ + beta*mass_/dt;
+    forAll(faceZoneIDs_, zoneI)
+    {
+        massTotal_[zoneI] += mass_[zoneI];
+        massFlux_[zoneI] = alpha*massFlux_[zoneI] + beta*mass_[zoneI]/dt;
+    }
 
     const label procI = Pstream::myProcNo();
 
-    scalarListList allProcMass(Pstream::nProcs());
-    allProcMass[procI] = massTotal_;
-    Pstream::gatherList(allProcMass);
-    scalarField allMass
-    (
-        ListListOps::combine<scalarList>
-        (
-            allProcMass, accessOp<scalarList>()
-        )
-    );
+    Info<< "particleFaceFlux output:" << nl;
 
-    scalarListList allProcMassFlux(Pstream::nProcs());
-    allProcMassFlux[procI] = massFlux_;
-    Pstream::gatherList(allProcMassFlux);
-    scalarField allMassFlux
-    (
-        ListListOps::combine<scalarList>
-        (
-            allProcMassFlux, accessOp<scalarList>()
-        )
-    );
+    List<scalarField> zoneMassTotal(mass_.size());
+    forAll(zoneMassTotal, zoneI)
+    {
+        scalarListList allProcMass(Pstream::nProcs());
+        allProcMass[procI] = massTotal_[zoneI];
+        Pstream::gatherList(allProcMass);
+        zoneMassTotal[zoneI] =
+            ListListOps::combine<scalarList>
+            (
+                allProcMass, accessOp<scalarList>()
+            );
 
-    Info<< "particleFaceFlux output:" << nl
-        << "    total mass      = " << sum(allMass) << nl
-        << "    average mass flux = " << sum(allMassFlux) << nl << endl;
+        const word& zoneName = fzm[faceZoneIDs_[zoneI]].name();
+        Info<< "    " << zoneName << " total mass      = "
+            << sum(zoneMassTotal[zoneI]) << nl;
+    }
+
+    List<scalarField> zoneMassFlux(massFlux_.size());
+    forAll(zoneMassFlux, zoneI)
+    {
+        scalarListList allProcMassFlux(Pstream::nProcs());
+        allProcMassFlux[procI] = massFlux_[zoneI];
+        Pstream::gatherList(allProcMassFlux);
+        zoneMassFlux[zoneI] =
+            ListListOps::combine<scalarList>
+            (
+                allProcMassFlux, accessOp<scalarList>()
+            );
+
+        const word& zoneName = fzm[faceZoneIDs_[zoneI]].name();
+        Info<< "    " << zoneName << " average mass flux = "
+            << sum(zoneMassFlux[zoneI]) << nl;
+    }
+
+    Info<< endl;
 
 
     if (surfaceFormat_ != "none")
     {
-        labelList pointToGlobal;
-        labelList uniqueMeshPointLabels;
-        autoPtr<globalIndex> globalPointsPtr =
-            mesh.globalData().mergePoints
-            (
-                fZone_().meshPoints(),
-                fZone_().meshPointMap(),
-                pointToGlobal,
-                uniqueMeshPointLabels
-            );
+        fileName outputDir = mesh.time().path();
 
-        pointField uniquePoints(mesh.points(), uniqueMeshPointLabels);
-        List<pointField> allProcPoints(Pstream::nProcs());
-        allProcPoints[procI] = uniquePoints;
-        Pstream::gatherList(allProcPoints);
-
-        faceList faces(fZone_().localFaces());
-        forAll(faces, i)
+        if (Pstream::parRun())
         {
-            inplaceRenumber(pointToGlobal, faces[i]);
+            // Put in undecomposed case (Note: gives problems for
+            // distributed data running)
+            outputDir =
+                outputDir/".."/"postProcessing"/cloud::prefix/
+                this->owner().name()/mesh.time().timeName();
         }
-        List<faceList> allProcFaces(Pstream::nProcs());
-        allProcFaces[procI] = faces;
-        Pstream::gatherList(allProcFaces);
-
-        if (Pstream::master())
+        else
         {
-            pointField allPoints
-            (
-                ListListOps::combine<pointField>
+            outputDir =
+                outputDir/"postProcessing"/cloud::prefix/
+                this->owner().name()/mesh.time().timeName();
+        }
+
+        forAll(faceZoneIDs_, zoneI)
+        {
+            const faceZone& fZone = fzm[faceZoneIDs_[zoneI]];
+
+            labelList pointToGlobal;
+            labelList uniqueMeshPointLabels;
+            autoPtr<globalIndex> globalPointsPtr =
+                mesh.globalData().mergePoints
                 (
-                    allProcPoints, accessOp<pointField>()
-                )
-            );
+                    fZone().meshPoints(),
+                    fZone().meshPointMap(),
+                    pointToGlobal,
+                    uniqueMeshPointLabels
+                );
 
-            faceList allFaces
-            (
-                ListListOps::combine<faceList>
+            pointField uniquePoints(mesh.points(), uniqueMeshPointLabels);
+            List<pointField> allProcPoints(Pstream::nProcs());
+            allProcPoints[procI] = uniquePoints;
+            Pstream::gatherList(allProcPoints);
+
+            faceList faces(fZone().localFaces());
+            forAll(faces, i)
+            {
+                inplaceRenumber(pointToGlobal, faces[i]);
+            }
+            List<faceList> allProcFaces(Pstream::nProcs());
+            allProcFaces[procI] = faces;
+            Pstream::gatherList(allProcFaces);
+
+            if (Pstream::master())
+            {
+                pointField allPoints
                 (
-                    allProcFaces, accessOp<faceList>()
-                )
-            );
+                    ListListOps::combine<pointField>
+                    (
+                        allProcPoints, accessOp<pointField>()
+                    )
+                );
 
-            fileName outputDir = mesh.time().path();
+                faceList allFaces
+                (
+                    ListListOps::combine<faceList>
+                    (
+                        allProcFaces, accessOp<faceList>()
+                    )
+                );
 
-            if (Pstream::parRun())
-            {
-                // Put in undecomposed case (Note: gives problems for
-                // distributed data running)
-                outputDir =
-                    outputDir/".."/"postProcessing"/cloud::prefix/
-                    this->owner().name()/mesh.time().timeName();
+                autoPtr<surfaceWriter> writer(surfaceWriter::New(surfaceFormat_));
+                writer->write
+                (
+                    outputDir,
+                    fZone.name(),
+                    allPoints,
+                    allFaces,
+                    "massTotal",
+                    zoneMassTotal[zoneI],
+                    false
+                );
+
+                writer->write
+                (
+                    outputDir,
+                    fZone.name(),
+                    allPoints,
+                    allFaces,
+                    "massFlux",
+                    zoneMassFlux[zoneI],
+                    false
+                );
             }
-            else
-            {
-                outputDir =
-                    outputDir/"postProcessing"/cloud::prefix/
-                    this->owner().name()/mesh.time().timeName();
-            }
-
-            autoPtr<surfaceWriter> writer(surfaceWriter::New(surfaceFormat_));
-            writer->write
-            (
-                outputDir,
-                "massTotal",
-                allPoints,
-                allFaces,
-                "massTotal",
-                allMass,
-                false
-            );
-            writer->write
-            (
-                outputDir,
-                "massFlux",
-                allPoints,
-                allFaces,
-                "massFlux",
-                allMassFlux,
-                false
-            );
         }
     }
 
 
     if (resetOnWrite_)
     {
-        massFlux_ = 0.0;
+        forAll(faceZoneIDs_, zoneI)
+        {
+            massFlux_[zoneI] = 0.0;
+        }
         totalTime_ = 0.0;
     }
 
-    mass_ = 0.0;
+    forAll(mass_, zoneI)
+    {
+        mass_[zoneI] = 0.0;
+    }
 
     // writeProperties();
 }
@@ -202,7 +237,7 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
 )
 :
     CloudFunctionObject<CloudType>(dict, owner, typeName),
-    fZone_(owner.mesh().faceZones()[this->coeffDict().lookup("faceZone")]),
+    faceZoneIDs_(),
     surfaceFormat_(this->coeffDict().lookup("surfaceFormat")),
     resetOnWrite_(this->coeffDict().lookup("resetOnWrite")),
     totalTime_(0.0),
@@ -210,15 +245,32 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
     massTotal_(),
     massFlux_()
 {
-    label allFaces = returnReduce(fZone_().size(), sumOp<scalar>());
-    Info<< "        Number of faces = " << allFaces << endl;
+    wordList faceZoneNames(this->coeffDict().lookup("faceZones"));
+    mass_.setSize(faceZoneNames.size());
+    massTotal_.setSize(faceZoneNames.size());
+    massFlux_.setSize(faceZoneNames.size());
 
-    mass_.setSize(fZone_.size(), 0.0);
+    DynamicList<label> zoneIDs;
+    const faceZoneMesh& fzm = owner.mesh().faceZones();
+    forAll(faceZoneNames, i)
+    {
+        const word& zoneName = faceZoneNames[i];
+        label zoneI = fzm.findZoneID(zoneName);
+        if (zoneI != -1)
+        {
+            zoneIDs.append(zoneI);
+            const faceZone& fz = fzm[zoneI];
+            label nFaces = returnReduce(fz.size(), sumOp<label>());
+            mass_[i].setSize(nFaces, 0.0);
+            massTotal_[i].setSize(nFaces, 0.0);
+            massFlux_[i].setSize(nFaces, 0.0);
+            Info<< "        " << zoneName << " faces: " << nFaces;
+        }
+    }
+
+    faceZoneIDs_.transfer(zoneIDs);
 
     // readProperties(); AND initialise mass... fields
-
-    massTotal_.setSize(fZone_.size(), 0.0);
-    massFlux_.setSize(fZone_.size(), 0.0);
 }
 
 
@@ -229,7 +281,7 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
 )
 :
     CloudFunctionObject<CloudType>(pff),
-    fZone_(pff.fZone_.clone(pff.fZone_.zoneMesh())),
+    faceZoneIDs_(pff.faceZoneIDs_),
     surfaceFormat_(pff.surfaceFormat_),
     resetOnWrite_(pff.resetOnWrite_),
     totalTime_(pff.totalTime_),
@@ -266,11 +318,13 @@ void Foam::FacePostProcessing<CloudType>::postFace(const parcelType& p)
      || this->owner().solution().transient()
     )
     {
-        const label faceI = applyToFace(p.face());
+        label zoneI = -1;
+        label faceI = -1;
+        applyToFace(p.face(), zoneI, faceI);
 
-        if (faceI != -1)
+        if ((zoneI != -1) && (faceI != -1))
         {
-            mass_[faceI] += p.mass()*p.nParticle();
+            mass_[zoneI][faceI] += p.mass()*p.nParticle();
         }
     }
 }
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.H
index de239097e6..85a80b1a4a 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.H
@@ -64,8 +64,8 @@ class FacePostProcessing
             typedef typename CloudType::parcelType parcelType;
 
 
-        //- Face zone
-        const faceZone& fZone_;
+        //- Face zone IDs
+        labelList faceZoneIDs_;
 
         //- Surface output format
         const word surfaceFormat_;
@@ -77,19 +77,24 @@ class FacePostProcessing
         scalar totalTime_;
 
         //- Mass storage
-        scalarField mass_;
+        List<scalarField> mass_;
 
         //- Mass total storage
-        scalarField massTotal_;
+        List<scalarField> massTotal_;
 
         //- Mass flux storage
-        scalarField massFlux_;
+        List<scalarField> massFlux_;
 
 
     // Private Member Functions
 
-        //- Return index into massFlux_ list if valid face, else -1
-        label applyToFace(const label faceI) const;
+        //- Return index into storage lists if valid zone and face
+        void applyToFace
+        (
+            const label faceIn,
+            label& zoneI, label&
+            faceI
+        ) const;
 
 
 protected: