diff --git a/src/functionObjects/field/fieldAverage/fieldAverage.C b/src/functionObjects/field/fieldAverage/fieldAverage.C
index 12ffdcd9a5f7f5d989fba9bd62207b4c59367e05..ad1ce8d969194f48155268859b034274a446b608 100644
--- a/src/functionObjects/field/fieldAverage/fieldAverage.C
+++ b/src/functionObjects/field/fieldAverage/fieldAverage.C
@@ -48,17 +48,17 @@ void Foam::functionObjects::fieldAverage::resetFields()
     {
         if (faItems_[i].mean())
         {
-            if (obr_.found(faItems_[i].meanFieldName()))
+            if (obr().found(faItems_[i].meanFieldName()))
             {
-                obr_.checkOut(*obr_[faItems_[i].meanFieldName()]);
+                obr().checkOut(*obr()[faItems_[i].meanFieldName()]);
             }
         }
 
         if (faItems_[i].prime2Mean())
         {
-            if (obr_.found(faItems_[i].prime2MeanFieldName()))
+            if (obr().found(faItems_[i].prime2MeanFieldName()))
             {
-                obr_.checkOut(*obr_[faItems_[i].prime2MeanFieldName()]);
+                obr().checkOut(*obr()[faItems_[i].prime2MeanFieldName()]);
             }
         }
     }
@@ -108,14 +108,14 @@ void Foam::functionObjects::fieldAverage::initialize()
 
 void Foam::functionObjects::fieldAverage::restart()
 {
-    Log << "    Restarting averaging at time " << obr_.time().timeName()
+    Log << "    Restarting averaging at time " << obr().time().timeName()
         << nl << endl;
 
     totalIter_.clear();
     totalIter_.setSize(faItems_.size(), 1);
 
     totalTime_.clear();
-    totalTime_.setSize(faItems_.size(), obr_.time().deltaTValue());
+    totalTime_.setSize(faItems_.size(), obr().time().deltaTValue());
 
     initialize();
 }
@@ -128,8 +128,8 @@ void Foam::functionObjects::fieldAverage::calcAverages()
         initialize();
     }
 
-    const label currentTimeIndex = obr_.time().timeIndex();
-    const scalar currentTime = obr_.time().value();
+    const label currentTimeIndex = obr().time().timeIndex();
+    const scalar currentTime = obr().time().value();
 
     if (prevTimeIndex_ == currentTimeIndex)
     {
@@ -176,7 +176,7 @@ void Foam::functionObjects::fieldAverage::calcAverages()
     forAll(faItems_, fieldi)
     {
         totalIter_[fieldi]++;
-        totalTime_[fieldi] += obr_.time().deltaTValue();
+        totalTime_[fieldi] += obr().time().deltaTValue();
     }
 
     Log << endl;
@@ -217,11 +217,11 @@ void Foam::functionObjects::fieldAverage::readAveragingProperties()
     totalIter_.setSize(faItems_.size(), 1);
 
     totalTime_.clear();
-    totalTime_.setSize(faItems_.size(), obr_.time().deltaTValue());
+    totalTime_.setSize(faItems_.size(), obr().time().deltaTValue());
 
     if (restartOnRestart_ || restartOnOutput_)
     {
-        Info<< "    Starting averaging at time " << obr_.time().timeName()
+        Info<< "    Starting averaging at time " << obr().time().timeName()
             << nl;
     }
     else
@@ -248,7 +248,7 @@ void Foam::functionObjects::fieldAverage::readAveragingProperties()
             {
                 Info<< "        " << fieldName
                     << ": starting averaging at time "
-                    << obr_.time().timeName() << endl;
+                    << obr().time().timeName() << endl;
             }
         }
     }
diff --git a/src/functionObjects/field/fieldAverage/fieldAverage.H b/src/functionObjects/field/fieldAverage/fieldAverage.H
index 20c8380d7eecb302e1115a2a795bb180309977da..74c4768d3775ce84e1407c92f9a268a83e6a5739 100644
--- a/src/functionObjects/field/fieldAverage/fieldAverage.H
+++ b/src/functionObjects/field/fieldAverage/fieldAverage.H
@@ -30,6 +30,7 @@ Group
 Description
     Calculates average quantities for a user-specified selection of volumetric
     and surface fields.
+    With the %subRegion option, also supports fields on a surfMesh.
 
     Fields are entered as a list of sub-dictionaries, which indicate the type of
     averages to perform, and can be updated during the calculation.  The current
@@ -108,6 +109,7 @@ Usage
         restartPeriod     | Periodic restart period | conditional |
         restartTime       | One-shot reset of the averaging | no | great
         fields            | list of fields and averaging options | yes |
+        subRegion         | name of a sub-region such as a surface name | no |
     \endtable
 
 
@@ -272,10 +274,10 @@ protected:
 
 
         //- Disallow default bitwise copy construct
-        fieldAverage(const fieldAverage&);
+        fieldAverage(const fieldAverage&) = delete;
 
         //- Disallow default bitwise assignment
-        void operator=(const fieldAverage&);
+        void operator=(const fieldAverage&) = delete;
 
 
 public:
diff --git a/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C b/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C
index fd6c50542d03b8c03bd96090bd5e2ce480bd1a53..422a1499391db06f4a70353cdf9d474e6ea843ee 100644
--- a/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C
+++ b/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,6 +26,7 @@ License
 #include "fieldAverageItem.H"
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "surfFields.H"
 #include "OFstream.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -41,11 +42,11 @@ void Foam::functionObjects::fieldAverage::addMeanFieldType(const label fieldi)
 
     Log << "    Reading/initialising field " << meanFieldName << endl;
 
-    if (obr_.foundObject<Type>(meanFieldName))
+    if (foundObject<Type>(meanFieldName))
     {
        // do nothing
     }
-    else if (obr_.found(meanFieldName))
+    else if (obr().found(meanFieldName))
     {
         Log << "    Cannot allocate average field " << meanFieldName
             << " since an object with that name already exists."
@@ -55,18 +56,18 @@ void Foam::functionObjects::fieldAverage::addMeanFieldType(const label fieldi)
     }
     else
     {
-        const Type& baseField = obr_.lookupObject<Type>(fieldName);
+        const Type& baseField = lookupObject<Type>(fieldName);
 
         // Store on registry
-        obr_.store
+        obr().store
         (
             new Type
             (
                 IOobject
                 (
                     meanFieldName,
-                    obr_.time().timeName(obr_.time().startTime().value()),
-                    obr_,
+                    obr().time().timeName(obr().time().startTime().value()),
+                    obr(),
                     restartOnOutput_
                   ? IOobject::NO_READ
                   : IOobject::READ_IF_PRESENT,
@@ -82,24 +83,26 @@ void Foam::functionObjects::fieldAverage::addMeanFieldType(const label fieldi)
 template<class Type>
 void Foam::functionObjects::fieldAverage::addMeanField(const label fieldi)
 {
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
+
     if (faItems_[fieldi].mean())
     {
-        typedef GeometricField<Type, fvPatchField, volMesh>
-            VolFieldType;
-
-        typedef GeometricField<Type, fvsPatchField, surfaceMesh>
-            SurfaceFieldType;
-
         const word& fieldName = faItems_[fieldi].fieldName();
 
-        if (obr_.foundObject<VolFieldType>(fieldName))
+        if (foundObject<VolFieldType>(fieldName))
         {
             addMeanFieldType<VolFieldType>(fieldi);
         }
-        else if (obr_.foundObject<SurfaceFieldType>(fieldName))
+        else if (foundObject<SurfaceFieldType>(fieldName))
         {
             addMeanFieldType<SurfaceFieldType>(fieldi);
         }
+        else if (foundObject<SurfFieldType>(fieldName))
+        {
+            addMeanFieldType<SurfFieldType>(fieldi);
+        }
     }
 }
 
@@ -116,11 +119,11 @@ void Foam::functionObjects::fieldAverage::addPrime2MeanFieldType
 
     Log << "    Reading/initialising field " << prime2MeanFieldName << nl;
 
-    if (obr_.foundObject<Type2>(prime2MeanFieldName))
+    if (foundObject<Type2>(prime2MeanFieldName))
     {
         // do nothing
     }
-    else if (obr_.found(prime2MeanFieldName))
+    else if (obr().found(prime2MeanFieldName))
     {
         Log << "    Cannot allocate average field " << prime2MeanFieldName
             << " since an object with that name already exists."
@@ -130,19 +133,19 @@ void Foam::functionObjects::fieldAverage::addPrime2MeanFieldType
     }
     else
     {
-        const Type1& baseField = obr_.lookupObject<Type1>(fieldName);
-        const Type1& meanField = obr_.lookupObject<Type1>(meanFieldName);
+        const Type1& baseField = lookupObject<Type1>(fieldName);
+        const Type1& meanField = lookupObject<Type1>(meanFieldName);
 
         // Store on registry
-        obr_.store
+        obr().store
         (
             new Type2
             (
                 IOobject
                 (
                     prime2MeanFieldName,
-                    obr_.time().timeName(obr_.time().startTime().value()),
-                    obr_,
+                    obr().time().timeName(obr().time().startTime().value()),
+                    obr(),
                     restartOnOutput_
                   ? IOobject::NO_READ
                   : IOobject::READ_IF_PRESENT,
@@ -160,9 +163,11 @@ void Foam::functionObjects::fieldAverage::addPrime2MeanField(const label fieldi)
 {
     typedef GeometricField<Type1, fvPatchField, volMesh> VolFieldType1;
     typedef GeometricField<Type1, fvsPatchField, surfaceMesh> SurfaceFieldType1;
+    typedef DimensionedField<Type1, surfGeoMesh> SurfFieldType1;
 
     typedef GeometricField<Type2, fvPatchField, volMesh> VolFieldType2;
     typedef GeometricField<Type2, fvsPatchField, surfaceMesh> SurfaceFieldType2;
+    typedef DimensionedField<Type2, surfGeoMesh> SurfFieldType2;
 
     if (faItems_[fieldi].prime2Mean())
     {
@@ -176,17 +181,24 @@ void Foam::functionObjects::fieldAverage::addPrime2MeanField(const label fieldi)
                 << fieldName << nl << exit(FatalError);
         }
 
-        if (obr_.foundObject<VolFieldType1>(fieldName))
+        if (foundObject<VolFieldType1>(fieldName))
         {
             addPrime2MeanFieldType<VolFieldType1, VolFieldType2>(fieldi);
         }
-        else if (obr_.foundObject<SurfaceFieldType1>(fieldName))
+        else if (foundObject<SurfaceFieldType1>(fieldName))
         {
             addPrime2MeanFieldType<SurfaceFieldType1, SurfaceFieldType2>
             (
                 fieldi
             );
         }
+        else if (foundObject<SurfFieldType1>(fieldName))
+        {
+            addPrime2MeanFieldType<SurfFieldType1, SurfFieldType2>
+            (
+                fieldi
+            );
+        }
     }
 }
 
@@ -199,16 +211,16 @@ void Foam::functionObjects::fieldAverage::calculateMeanFieldType
 {
     const word& fieldName = faItems_[fieldi].fieldName();
 
-    if (obr_.foundObject<Type>(fieldName))
+    if (foundObject<Type>(fieldName))
     {
-        const Type& baseField = obr_.lookupObject<Type>(fieldName);
+        const Type& baseField = lookupObject<Type>(fieldName);
 
         Type& meanField = const_cast<Type&>
         (
-            obr_.lookupObject<Type>(faItems_[fieldi].meanFieldName())
+            lookupObject<Type>(faItems_[fieldi].meanFieldName())
         );
 
-        scalar dt = obr_.time().deltaTValue();
+        scalar dt = obr().time().deltaTValue();
         scalar Dt = totalTime_[fieldi];
 
         if (faItems_[fieldi].iterBase())
@@ -239,6 +251,7 @@ void Foam::functionObjects::fieldAverage::calculateMeanFields() const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
 
     forAll(faItems_, i)
     {
@@ -246,6 +259,7 @@ void Foam::functionObjects::fieldAverage::calculateMeanFields() const
         {
             calculateMeanFieldType<VolFieldType>(i);
             calculateMeanFieldType<SurfaceFieldType>(i);
+            calculateMeanFieldType<SurfFieldType>(i);
         }
     }
 }
@@ -259,18 +273,16 @@ void Foam::functionObjects::fieldAverage::calculatePrime2MeanFieldType
 {
     const word& fieldName = faItems_[fieldi].fieldName();
 
-    if (obr_.foundObject<Type1>(fieldName))
+    if (foundObject<Type1>(fieldName))
     {
-        const Type1& baseField = obr_.lookupObject<Type1>(fieldName);
+        const Type1& baseField = lookupObject<Type1>(fieldName);
         const Type1& meanField =
-            obr_.lookupObject<Type1>(faItems_[fieldi].meanFieldName());
+            lookupObject<Type1>(faItems_[fieldi].meanFieldName());
 
-        Type2& prime2MeanField = const_cast<Type2&>
-        (
-            obr_.lookupObject<Type2>(faItems_[fieldi].prime2MeanFieldName())
-        );
+        Type2& prime2MeanField =
+            lookupObjectRef<Type2>(faItems_[fieldi].prime2MeanFieldName());
 
-        scalar dt = obr_.time().deltaTValue();
+        scalar dt = obr().time().deltaTValue();
         scalar Dt = totalTime_[fieldi];
 
         if (faItems_[fieldi].iterBase())
@@ -304,9 +316,11 @@ void Foam::functionObjects::fieldAverage::calculatePrime2MeanFields() const
 {
     typedef GeometricField<Type1, fvPatchField, volMesh> VolFieldType1;
     typedef GeometricField<Type1, fvsPatchField, surfaceMesh> SurfaceFieldType1;
+    typedef DimensionedField<Type1, surfGeoMesh> SurfFieldType1;
 
     typedef GeometricField<Type2, fvPatchField, volMesh> VolFieldType2;
     typedef GeometricField<Type2, fvsPatchField, surfaceMesh> SurfaceFieldType2;
+    typedef DimensionedField<Type2, surfGeoMesh> SurfFieldType2;
 
     forAll(faItems_, i)
     {
@@ -317,6 +331,11 @@ void Foam::functionObjects::fieldAverage::calculatePrime2MeanFields() const
             (
                 i
             );
+
+            calculatePrime2MeanFieldType<SurfFieldType1, SurfFieldType2>
+            (
+                i
+            );
         }
     }
 }
@@ -330,15 +349,13 @@ void Foam::functionObjects::fieldAverage::addMeanSqrToPrime2MeanType
 {
     const word& fieldName = faItems_[fieldi].fieldName();
 
-    if (obr_.foundObject<Type1>(fieldName))
+    if (foundObject<Type1>(fieldName))
     {
         const Type1& meanField =
-            obr_.lookupObject<Type1>(faItems_[fieldi].meanFieldName());
+            lookupObject<Type1>(faItems_[fieldi].meanFieldName());
 
-        Type2& prime2MeanField = const_cast<Type2&>
-        (
-            obr_.lookupObject<Type2>(faItems_[fieldi].prime2MeanFieldName())
-        );
+        Type2& prime2MeanField =
+            lookupObjectRef<Type2>(faItems_[fieldi].prime2MeanFieldName());
 
         prime2MeanField += sqr(meanField);
     }
@@ -350,9 +367,11 @@ void Foam::functionObjects::fieldAverage::addMeanSqrToPrime2Mean() const
 {
     typedef GeometricField<Type1, fvPatchField, volMesh> VolFieldType1;
     typedef GeometricField<Type1, fvsPatchField, surfaceMesh> SurfaceFieldType1;
+    typedef DimensionedField<Type1, surfGeoMesh> SurfFieldType1;
 
     typedef GeometricField<Type2, fvPatchField, volMesh> VolFieldType2;
     typedef GeometricField<Type2, fvsPatchField, surfaceMesh> SurfaceFieldType2;
+    typedef DimensionedField<Type2, surfGeoMesh> SurfFieldType2;
 
     forAll(faItems_, i)
     {
@@ -360,6 +379,7 @@ void Foam::functionObjects::fieldAverage::addMeanSqrToPrime2Mean() const
         {
             addMeanSqrToPrime2MeanType<VolFieldType1, VolFieldType2>(i);
             addMeanSqrToPrime2MeanType<SurfaceFieldType1, SurfaceFieldType2>(i);
+            addMeanSqrToPrime2MeanType<SurfFieldType1, SurfFieldType2>(i);
         }
     }
 }
@@ -371,9 +391,9 @@ void Foam::functionObjects::fieldAverage::writeFieldType
     const word& fieldName
 ) const
 {
-    if (obr_.foundObject<Type>(fieldName))
+    if (foundObject<Type>(fieldName))
     {
-        const Type& f = obr_.lookupObject<Type>(fieldName);
+        const Type& f = lookupObject<Type>(fieldName);
         f.write();
     }
 }
@@ -384,6 +404,7 @@ void Foam::functionObjects::fieldAverage::writeFields() const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
 
     forAll(faItems_, i)
     {
@@ -392,12 +413,14 @@ void Foam::functionObjects::fieldAverage::writeFields() const
             const word& fieldName = faItems_[i].meanFieldName();
             writeFieldType<VolFieldType>(fieldName);
             writeFieldType<SurfaceFieldType>(fieldName);
+            writeFieldType<SurfFieldType>(fieldName);
         }
         if (faItems_[i].prime2Mean())
         {
             const word& fieldName = faItems_[i].prime2MeanFieldName();
             writeFieldType<VolFieldType>(fieldName);
             writeFieldType<SurfaceFieldType>(fieldName);
+            writeFieldType<SurfFieldType>(fieldName);
         }
     }
 }
diff --git a/src/functionObjects/field/fieldsExpression/fieldsExpressionTemplates.C b/src/functionObjects/field/fieldsExpression/fieldsExpressionTemplates.C
index b9b3be11359d67230518659125f9575756f93507..c50b7ed98b72087a105187fe3260cccb393b85e1 100644
--- a/src/functionObjects/field/fieldsExpression/fieldsExpressionTemplates.C
+++ b/src/functionObjects/field/fieldsExpression/fieldsExpressionTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,6 +25,7 @@ License
 
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "surfFields.H"
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
@@ -33,6 +34,7 @@ bool Foam::functionObjects::fieldsExpression::calcFieldTypes(FOType& fo)
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
 
     if (foundObject<VolFieldType>(fieldNames_[0]))
     {
@@ -50,6 +52,14 @@ bool Foam::functionObjects::fieldsExpression::calcFieldTypes(FOType& fo)
             fo.template calcFieldType<SurfaceFieldType>()
         );
     }
+    else if (foundObject<SurfFieldType>(fieldNames_[0]))
+    {
+        return store
+        (
+            resultName_,
+            fo.template calcFieldType<SurfFieldType>()
+        );
+    }
     else
     {
         return false;
diff --git a/src/functionObjects/field/mag/mag.H b/src/functionObjects/field/mag/mag.H
index de2737c31d2294140805b97051f04c83c74e0d5c..2010e25e5a1876bf2374e67b8185b7c0c5fec3e9 100644
--- a/src/functionObjects/field/mag/mag.H
+++ b/src/functionObjects/field/mag/mag.H
@@ -32,6 +32,7 @@ Description
 
     The operation can be applied to any volume or surface fields generating a
     volume or surface scalar field.
+    With the %subRegion option, also supports fields on a surfMesh.
 
 See also
     Foam::functionObjects::fvMeshFunctionObject
diff --git a/src/functionObjects/field/mag/magTemplates.C b/src/functionObjects/field/mag/magTemplates.C
index 0b6035f0a0395352f148416d2748309d1419cda2..dead35a360e6d65e1726e77b020f2853420872e4 100644
--- a/src/functionObjects/field/mag/magTemplates.C
+++ b/src/functionObjects/field/mag/magTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,6 +25,7 @@ License
 
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "surfFields.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -33,6 +34,7 @@ bool Foam::functionObjects::mag::calcMag()
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
 
     if (foundObject<VolFieldType>(fieldName_))
     {
@@ -50,6 +52,14 @@ bool Foam::functionObjects::mag::calcMag()
             Foam::mag(lookupObject<SurfaceFieldType>(fieldName_))
         );
     }
+    else if (foundObject<SurfFieldType>(fieldName_))
+    {
+        return store
+        (
+            resultName_,
+            Foam::mag(lookupObject<SurfFieldType>(fieldName_))
+        );
+    }
     else
     {
         return false;
diff --git a/src/functionObjects/field/magSqr/magSqr.H b/src/functionObjects/field/magSqr/magSqr.H
index f1af984178aad273ae8b95b4e9fd2c36c825e1bd..c112187e70e0744113b7d70b461c49d733339800 100644
--- a/src/functionObjects/field/magSqr/magSqr.H
+++ b/src/functionObjects/field/magSqr/magSqr.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +32,7 @@ Description
 
     The operation can be applied to any volume or surface field generating a
     volume or surface scalar field.
+    With the %subRegion option, also supports fields on a surfMesh.
 
 See also
     Foam::functionObjects::fvMeshFunctionObject
diff --git a/src/functionObjects/field/magSqr/magSqrTemplates.C b/src/functionObjects/field/magSqr/magSqrTemplates.C
index a5b1011c395da901838a38e99a7f8f743d6ea2b9..a3fee019bb30614e410592ddd31a52919ec520ed 100644
--- a/src/functionObjects/field/magSqr/magSqrTemplates.C
+++ b/src/functionObjects/field/magSqr/magSqrTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,6 +25,7 @@ License
 
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "surfFields.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -33,6 +34,7 @@ bool Foam::functionObjects::magSqr::calcMagSqr()
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
 
     if (foundObject<VolFieldType>(fieldName_))
     {
@@ -50,6 +52,14 @@ bool Foam::functionObjects::magSqr::calcMagSqr()
             Foam::magSqr(lookupObject<SurfaceFieldType>(fieldName_))
         );
     }
+    else if (foundObject<SurfFieldType>(fieldName_))
+    {
+        return store
+        (
+            resultName_,
+            Foam::magSqr(lookupObject<SurfFieldType>(fieldName_))
+        );
+    }
     else
     {
         return false;
diff --git a/src/functionObjects/field/readFields/readFieldsTemplates.C b/src/functionObjects/field/readFields/readFieldsTemplates.C
index 07d73044c796d8e54dc70c286dc832a1a003cfcd..326b8f2156215c0c790c254e98c2f9de836e5919 100644
--- a/src/functionObjects/field/readFields/readFieldsTemplates.C
+++ b/src/functionObjects/field/readFields/readFieldsTemplates.C
@@ -26,6 +26,7 @@ License
 #include "readFields.H"
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "surfFields.H"
 #include "Time.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -35,6 +36,7 @@ bool Foam::functionObjects::readFields::loadField(const word& fieldName)
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+    typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
 
     if (foundObject<VolFieldType>(fieldName))
     {
@@ -49,6 +51,12 @@ bool Foam::functionObjects::readFields::loadField(const word& fieldName)
             << " " << fieldName << " already exists in database"
             << " already in database" << endl;
     }
+    else if (foundObject<SurfFieldType>(fieldName))
+    {
+        DebugInfo<< "readFields: " << SurfFieldType::typeName
+            << " " << fieldName << " already exists in database"
+            << " already in database" << endl;
+    }
     else
     {
         IOobject fieldHeader
@@ -76,6 +84,23 @@ bool Foam::functionObjects::readFields::loadField(const word& fieldName)
             mesh_.objectRegistry::store(sfPtr);
             return true;
         }
+        else if (fieldHeader.typeHeaderOk<SurfFieldType>(true))
+        {
+            if (isA<surfMesh>(obr()))
+            {
+                const surfMesh& s = dynamicCast<const surfMesh>(obr());
+
+                // Store field on surfMesh database
+                Log << "    Reading " << fieldName << endl;
+                SurfFieldType* sfPtr(new SurfFieldType(fieldHeader, s));
+                s.store(sfPtr);
+                return true;
+            }
+            else
+            {
+                return false;
+            }
+        }
     }
 
     return false;