diff --git a/applications/utilities/postProcessing/noise/Make/options b/applications/utilities/postProcessing/noise/Make/options
index b29b939580..2fc476e741 100644
--- a/applications/utilities/postProcessing/noise/Make/options
+++ b/applications/utilities/postProcessing/noise/Make/options
@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/randomProcesses/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude
 
 EXE_LIBS = \
     -lrandomProcesses \
diff --git a/applications/utilities/postProcessing/noise/noise.C b/applications/utilities/postProcessing/noise/noise.C
index 1b0ca07fec..5c3ad2a139 100644
--- a/applications/utilities/postProcessing/noise/noise.C
+++ b/applications/utilities/postProcessing/noise/noise.C
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,157 +28,104 @@ Group
     grpPostProcessingUtilities
 
 Description
-    Utility to perform noise analysis of pressure data using the noiseFFT
-    library.
+    Utility to perform noise analysis of pressure data.
 
-    Control settings are read from the $FOAM_CASE/system/noiseDict dictionary,
-    or user-specified dictionary using the -dict option.  Pressure data is
-    read using a CSV reader:
+    The utility provides a light wrapper around the run-time selectable
+    noise model.  Current options include:
+    - point, and
+    - surface noise.
 
-    \heading Usage
+    \heading Example usage
 
     \verbatim
-    pRef        101325;
-    N           65536;
-    nw          100;
-    f1          25;
-    fU          10000;
-    graphFormat raw;
+    noiseModel      surfaceNoise; // pointNoise
 
-    pressureData
+    surfaceNoiseCoeffs
     {
-        fileName        "pressureData"
-        nHeaderLine         1;          // number of header lines
-        refColumn           0;          // reference column index
-        componentColumns    (1);        // component column indices
-        separator           " ";        // optional (defaults to ",")
-        mergeSeparators     no;         // merge multiple separators
-        outOfBounds         clamp;      // optional out-of-bounds handling
-        interpolationScheme linear;     // optional interpolation scheme
+        windowModel     Hanning;
+
+        HanningCoeffs
+        {
+            // Window overlap percentage
+            overlapPercent  50;
+            symmetric       yes;
+            extended        yes;
+
+            // Optional number of windows, default = all available
+            nWindow         5;
+        }
+
+
+        // Input file
+        inputFile   "postProcessing/faceSource1/surface/patch/patch.case";
+
+        // Surface reader
+        reader      ensight;
+
+        // Surface writer
+        writer      ensight;
+
+        // Collate times for ensight output - ensures geometry is only written once
+        writeOptions
+        {
+            ensight
+            {
+                collateTimes 1;
+            }
+        }
+
+        // Number of samples in sampling window
+        // Must be a power of 2, default = 2^16 (=65536)
+        N               4096;
+
+        // Write interval for FFT data, default = 1
+        fftWriteInterval 100;
     }
     \endverbatim
 
-    where
-    \table
-        Property    | Description                   | Required  | Default value
-        pRef        | Reference pressure            | no        | 0
-        N           | Number of samples in sampling window | no | 65536
-        nw          | Number of sampling windows    | no        | 100
-        fl          | Lower frequency band          | no        | 25
-        fU          | Upper frequency band          | no        | 10000
-        graphFormat | Output graph format          | no        | raw
-    \endtable
-
-    Current graph outputs include:
-    - FFT of the pressure data
-    - narrow-band PFL (pressure-fluctuation level) spectrum
-    - one-third-octave-band PFL spectrum
-    - one-third-octave-band pressure spectrum
 
 SeeAlso
-    CSV.H
     noiseFFT.H
+    noiseModel.H
+    windowModel.H
 
 \*---------------------------------------------------------------------------*/
 
-
-#include "noiseFFT.H"
 #include "argList.H"
 #include "Time.H"
-#include "functionObjectFile.H"
-#include "CSV.H"
+#include "noiseModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 using namespace Foam;
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-Foam::scalar checkUniformTimeStep(const scalarField& t)
-{
-    // check that a uniform time step has been applied
-    scalar deltaT = -1.0;
-    if (t.size() > 1)
-    {
-        for (label i = 1; i < t.size(); i++)
-        {
-            scalar dT = t[i] - t[i-1];
-            if (deltaT < 0)
-            {
-                deltaT = dT;
-            }
-
-            if (mag(deltaT - dT) > SMALL)
-            {
-                FatalErrorInFunction
-                    << "Unable to process data with a variable time step"
-                    << exit(FatalError);
-            }
-        }
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "Unable to create FFT with a single value"
-            << exit(FatalError);
-    }
-
-    return deltaT;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
-    argList::noParallel();
     #include "addDictOption.H"
     #include "setRootCase.H"
     #include "createTime.H"
-    #include "createFields.H"
 
-    Info<< "Reading data file" << endl;
-    Function1Types::CSV<scalar> pData("pressure", dict, "Data");
-
-    // time history data
-    const scalarField t(pData.x());
-
-    // pressure data
-    const scalarField p(pData.y());
-
-    if (t.size() < N)
+    word dictName("noiseDict");
+    if (args.optionFound("dict"))
     {
-        FatalErrorInFunction
-            << "Block size N = " << N
-            << " is larger than number of data = " << t.size()
-            << exit(FatalError);
+        dictName = args["dict"];
     }
 
-    Info<< "    read " << t.size() << " values" << nl << endl;
+    IOdictionary dict
+    (
+        IOobject
+        (
+            dictName,
+            runTime.system(),
+            runTime,
+            IOobject::MUST_READ
+        )
+    );
 
-
-    Info<< "Creating noise FFT" << endl;
-    noiseFFT nfft(checkUniformTimeStep(t), p);
-
-    nfft -= pRef;
-
-    fileName baseFileName(pData.fName().lessExt());
-
-    graph Pf(nfft.RMSmeanPf(N, min(nfft.size()/N, nw)));
-    Info<< "    Creating graph for " << Pf.title() << endl;
-    Pf.write(baseFileName + graph::wordify(Pf.title()), graphFormat);
-
-    graph Lf(nfft.Lf(Pf));
-    Info<< "    Creating graph for " << Lf.title() << endl;
-    Lf.write(baseFileName + graph::wordify(Lf.title()), graphFormat);
-
-    graph Ldelta(nfft.Ldelta(Lf, f1, fU));
-    Info<< "    Creating graph for " << Ldelta.title() << endl;
-    Ldelta.write(baseFileName + graph::wordify(Ldelta.title()), graphFormat);
-
-    graph Pdelta(nfft.Pdelta(Pf, f1, fU));
-    Info<< "    Creating graph for " << Pdelta.title() << endl;
-    Pdelta.write(baseFileName + graph::wordify(Pdelta.title()), graphFormat);
+    autoPtr<noiseModel> model(noiseModel::New(dict));
+    model->calculate();
 
     Info<< nl << "End\n" << endl;
 
diff --git a/etc/cshrc b/etc/cshrc
index 30138fdb36..f45c2174a0 100644
--- a/etc/cshrc
+++ b/etc/cshrc
@@ -31,7 +31,7 @@
 #------------------------------------------------------------------------------
 
 setenv WM_PROJECT OpenFOAM
-setenv WM_PROJECT_VERSION plus
+setenv WM_PROJECT_VERSION stage
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
diff --git a/src/randomProcesses/Make/files b/src/randomProcesses/Make/files
index 90cb360357..a16bde251c 100644
--- a/src/randomProcesses/Make/files
+++ b/src/randomProcesses/Make/files
@@ -1,25 +1,31 @@
-Kmesh                       = Kmesh
-
-fft                         = fft
-
-processes                   = processes
-UOprocess                   = $(processes)/UOprocess
-
-turbulence                  = turbulence
-
-noise                       = noise
-
+Kmesh = Kmesh
 $(Kmesh)/Kmesh.C
 
+fft = fft
 $(fft)/fft.C
 $(fft)/fftRenumber.C
 $(fft)/calcEk.C
 $(fft)/kShellIntegration.C
 
+processes = processes
+UOprocess = $(processes)/UOprocess
 $(UOprocess)/UOprocess.C
 
+turbulence = turbulence
 $(turbulence)/turbGen.C
 
-$(noise)/noiseFFT.C
+noise = noise
+$(noise)/noiseFFT/noiseFFT.C
+$(noise)/noiseModels/noiseModel/noiseModel.C
+$(noise)/noiseModels/noiseModel/noiseModelNew.C
+$(noise)/noiseModels/pointNoise/pointNoise.C
+$(noise)/noiseModels/surfaceNoise/surfaceNoise.C
+
+
+windowModels = windowModels
+$(windowModels)/windowModel/windowModel.C
+$(windowModels)/windowModel/windowModelNew.C
+$(windowModels)/Hanning/Hanning.C
+
 
 LIB = $(FOAM_LIBBIN)/librandomProcesses
diff --git a/src/randomProcesses/Make/options b/src/randomProcesses/Make/options
index 71b7873964..86d0d7baee 100644
--- a/src/randomProcesses/Make/options
+++ b/src/randomProcesses/Make/options
@@ -1,5 +1,9 @@
 EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude
 
 LIB_LIBS = \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lsampling \
+    -lsurfMesh
diff --git a/src/randomProcesses/noise/noiseFFT.C b/src/randomProcesses/noise/noiseFFT.C
deleted file mode 100644
index 0f09cdf817..0000000000
--- a/src/randomProcesses/noise/noiseFFT.C
+++ /dev/null
@@ -1,454 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "noiseFFT.H"
-#include "IFstream.H"
-#include "DynamicList.H"
-#include "fft.H"
-#include "SubField.H"
-#include "mathematicalConstants.H"
-
-// * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * //
-
-Foam::scalar Foam::noiseFFT::p0 = 2e-5;
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::noiseFFT::noiseFFT
-(
-    const scalar deltat,
-    const scalarField& pressure
-)
-:
-    scalarField(pressure),
-    deltat_(deltat)
-{}
-
-
-Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
-:
-    scalarField(),
-    deltat_(0.0)
-{
-    // Construct pressure data file
-    IFstream pFile(pFileName);
-
-    // Check pFile stream is OK
-    if (!pFile.good())
-    {
-        FatalErrorInFunction
-            << "Cannot read file " << pFileName
-            << exit(FatalError);
-    }
-
-    if (skip)
-    {
-        scalar dummyt, dummyp;
-
-        for (label i=0; i<skip; i++)
-        {
-            pFile >> dummyt;
-
-            if (!pFile.good() || pFile.eof())
-            {
-                FatalErrorInFunction
-                    << "Number of points in file " << pFileName
-                    << " is less than the number to be skipped = " << skip
-                    << exit(FatalError);
-            }
-
-            pFile >> dummyp;
-        }
-    }
-
-    scalar t = 0, T = 0;
-    DynamicList<scalar> pData(100000);
-    label i = 0;
-
-    while (!(pFile >> t).eof())
-    {
-        T = t;
-        pFile >> pData(i++);
-    }
-
-    deltat_ = T/pData.size();
-
-    this->transfer(pData);
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::graph Foam::noiseFFT::pt() const
-{
-    scalarField t(size());
-    forAll(t, i)
-    {
-        t[i] = i*deltat_;
-    }
-
-    return graph
-    (
-        "p(t)",
-        "t [s]",
-        "p(t) [Pa]",
-        t,
-        *this
-    );
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::noiseFFT::window
-(
-    const label N,
-    const label ni
-) const
-{
-    label windowOffset = N;
-
-    if ((N + ni*windowOffset) > size())
-    {
-        FatalErrorInFunction
-            << "Requested window is outside set of data" << endl
-            << "number of data = " << size() << endl
-            << "size of window = " << N << endl
-            << "window  = " << ni
-            << exit(FatalError);
-    }
-
-    tmp<scalarField> tpw(new scalarField(N));
-    scalarField& pw = tpw.ref();
-
-    label offset = ni*windowOffset;
-
-    forAll(pw, i)
-    {
-        pw[i] = operator[](i + offset);
-    }
-
-    return tpw;
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::noiseFFT::Hanning(const label N) const
-{
-    scalarField t(N);
-    forAll(t, i)
-    {
-        t[i] = i*deltat_;
-    }
-
-    scalar T = N*deltat_;
-
-    return 2*(0.5 - 0.5*cos(constant::mathematical::twoPi*t/T));
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
-(
-    const tmp<scalarField>& tpn
-) const
-{
-    tmp<scalarField> tPn2
-    (
-        mag
-        (
-            fft::reverseTransform
-            (
-                ReComplexField(tpn),
-                labelList(1, tpn().size())
-            )
-        )
-    );
-
-    tpn.clear();
-
-    tmp<scalarField> tPn
-    (
-        new scalarField
-        (
-            scalarField::subField(tPn2(), tPn2().size()/2)
-        )
-    );
-    scalarField& Pn = tPn.ref();
-
-    Pn *= 2.0/sqrt(scalar(tPn2().size()));
-    Pn[0] /= 2.0;
-
-    return tPn;
-}
-
-
-Foam::graph Foam::noiseFFT::meanPf
-(
-    const label N,
-    const label nw
-) const
-{
-    if (N > size())
-    {
-        FatalErrorInFunction
-            << "Requested window is outside set of data" << nl
-            << "number of data = " << size() << nl
-            << "size of window = " << N << nl
-            << "window  = " << nw
-            << exit(FatalError);
-    }
-
-    scalarField MeanPf(N/2, 0.0);
-
-    scalarField Hwf(Hanning(N));
-
-    for (label wi=0; wi<nw; ++wi)
-    {
-        MeanPf += Pf(Hwf*window(N, wi));
-    }
-
-    MeanPf /= nw;
-
-    scalarField f(MeanPf.size());
-    scalar deltaf = 1.0/(N*deltat_);
-
-    forAll(f, i)
-    {
-        f[i] = i*deltaf;
-    }
-
-    return graph
-    (
-        "P(f)",
-        "f [Hz]",
-        "P(f) [Pa]",
-        f,
-        MeanPf
-    );
-}
-
-
-Foam::graph Foam::noiseFFT::RMSmeanPf
-(
-    const label N,
-    const label nw
-) const
-{
-    if (N > size())
-    {
-        FatalErrorInFunction
-            << "Requested window is outside set of data" << endl
-            << "number of data = " << size() << endl
-            << "size of window = " << N << endl
-            << "window  = " << nw
-            << exit(FatalError);
-    }
-
-    scalarField RMSMeanPf(N/2, 0.0);
-
-    scalarField Hwf(Hanning(N));
-
-    for (label wi=0; wi<nw; ++wi)
-    {
-        RMSMeanPf += sqr(Pf(Hwf*window(N, wi)));
-    }
-
-    RMSMeanPf = sqrt(RMSMeanPf/nw);
-
-    scalarField f(RMSMeanPf.size());
-    scalar deltaf = 1.0/(N*deltat_);
-
-    forAll(f, i)
-    {
-        f[i] = i*deltaf;
-    }
-
-    return graph
-    (
-        "P(f)",
-        "f [Hz]",
-        "P(f) [Pa]",
-        f,
-        RMSMeanPf
-    );
-}
-
-
-Foam::graph Foam::noiseFFT::Lf(const graph& gPf) const
-{
-    return graph
-    (
-        "L(f)",
-        "f [Hz]",
-        "L(f) [dB]",
-        gPf.x(),
-        20*log10(gPf.y()/p0)
-    );
-}
-
-
-Foam::graph Foam::noiseFFT::Ldelta
-(
-    const graph& gLf,
-    const scalar f1,
-    const scalar fU
-) const
-{
-    const scalarField& f = gLf.x();
-    const scalarField& Lf = gLf.y();
-
-    scalarField ldelta(Lf.size(), 0.0);
-    scalarField fm(ldelta.size());
-
-    scalar fratio = cbrt(2.0);
-    scalar deltaf = 1.0/(2*Lf.size()*deltat_);
-
-    scalar fl = f1/sqrt(fratio);
-    scalar fu = fratio*fl;
-
-    label istart = label(fl/deltaf);
-    label j = 0;
-
-    for (label i = istart; i<Lf.size(); i++)
-    {
-        scalar fmi = sqrt(fu*fl);
-
-        if (fmi > fU + 1) break;
-
-        if (f[i] >= fu)
-        {
-            fm[j] = fmi;
-            ldelta[j] = 10*log10(ldelta[j]);
-
-            j++;
-
-            fl = fu;
-            fu *= fratio;
-        }
-
-        ldelta[j] += pow(10, Lf[i]/10.0);
-    }
-
-    fm.setSize(j);
-    ldelta.setSize(j);
-
-    return graph
-    (
-        "Ldelta",
-        "fm [Hz]",
-        "Ldelta [dB]",
-        fm,
-        ldelta
-    );
-}
-
-
-Foam::graph Foam::noiseFFT::Pdelta
-(
-    const graph& gPf,
-    const scalar f1,
-    const scalar fU
-) const
-{
-    const scalarField& f = gPf.x();
-    const scalarField& Pf = gPf.y();
-
-    scalarField pdelta(Pf.size(), 0.0);
-    scalarField fm(pdelta.size());
-
-    scalar fratio = cbrt(2.0);
-    scalar deltaf = 1.0/(2*Pf.size()*deltat_);
-
-    scalar fl = f1/sqrt(fratio);
-    scalar fu = fratio*fl;
-
-    label istart = label(fl/deltaf + 1.0 - SMALL);
-    label j = 0;
-
-    for (label i = istart; i<Pf.size(); i++)
-    {
-        scalar fmi = sqrt(fu*fl);
-
-        if (fmi > fU + 1) break;
-
-        if (f[i] >= fu)
-        {
-            fm[j] = fmi;
-            pdelta[j] = sqrt((2.0/3.0)*pdelta[j]);
-
-            j++;
-
-            fl = fu;
-            fu *= fratio;
-        }
-
-        pdelta[j] += sqr(Pf[i]);
-    }
-
-    fm.setSize(j);
-    pdelta.setSize(j);
-
-    return graph
-    (
-        "Pdelta",
-        "fm [Hz]",
-        "Pdelta [dB]",
-        fm,
-        pdelta
-    );
-}
-
-
-Foam::scalar Foam::noiseFFT::Lsum(const graph& gLf) const
-{
-    const scalarField& Lf = gLf.y();
-
-    scalar lsum = 0.0;
-
-    forAll(Lf, i)
-    {
-        lsum += pow(10, Lf[i]/10.0);
-    }
-
-    lsum = 10*log10(lsum);
-
-    return lsum;
-}
-
-
-Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
-{
-    return p0*pow(10.0, db/20.0);
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::noiseFFT::dbToPa
-(
-    const tmp<scalarField>& db
-) const
-{
-    return p0*pow(10.0, db/20.0);
-}
-
-
-// ************************************************************************* //
diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.C b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
new file mode 100644
index 0000000000..0d36c47d16
--- /dev/null
+++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
@@ -0,0 +1,519 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "noiseFFT.H"
+#include "IFstream.H"
+#include "DynamicList.H"
+#include "fft.H"
+#include "SubField.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::noiseFFT::p0 = 2e-5;
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseFFT::frequencies
+(
+    const label N,
+    const scalar deltaT
+)
+{
+    tmp<scalarField> tf(new scalarField(N/2, 0));
+    scalarField& f = tf();
+
+    scalar deltaf = 1.0/(N*deltaT);
+    forAll(f, i)
+    {
+        f[i] = i*deltaf;
+    }
+
+    return tf;
+}
+
+
+void Foam::noiseFFT::octaveFrequenciesIDs
+(
+    const scalarField& f,
+    const scalar fLower,
+    const scalar fUpper,
+    const scalar octave,
+    labelList& freqBandIDs
+)
+{
+    // Set the indices of to the lower frequency bands for the input frequency
+    // range. Ensure that the centre frequency passes though 1000 Hz
+
+    // Low frequency bound given by:
+    //     fLow = f0*(2^(bandI/octave/2))
+    // Centre frequency given by:
+    //     fCentre = f0*(2^(bandI/octave))
+
+    scalar f0 = 1000;
+    scalar minFrequency = max(fLower, min(f));
+
+    // Lower frequency band limit
+    label band0Low = ceil(2*octave*log(minFrequency/f0)/log(2.0));
+
+    // Centre frequency band limit
+    //label band0Centre = ceil(octave*log(fLower/f0)/log(2.0));
+
+    scalar fLowerBand = f0*pow(2, band0Low/octave/2);
+    scalar fRatio = pow(2, 1.0/octave);
+
+    bool complete = false;
+    DynamicList<label> bandIDs(f.size());
+    forAll(f, i)
+    {
+        while (f[i] >= fLowerBand)
+        {
+            bandIDs.append(i);
+            fLowerBand *= fRatio;
+
+            if (fLowerBand > fUpper)
+            {
+                complete = true;
+                break;
+            }
+        }
+
+        if (complete) break;
+    }
+
+    freqBandIDs.transfer(bandIDs);
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseFFT::octaveFrequencies
+(
+    const scalarField& f,
+    const scalar fLower,
+    const scalar fUpper,
+    const scalar octave
+)
+{
+    labelList freqBandIDs;
+    octaveFrequenciesIDs(f, fLower, fUpper, octave, freqBandIDs);
+    tmp<scalarField> tf(new scalarField(f, freqBandIDs));
+    return tf;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::noiseFFT::noiseFFT
+(
+    const scalar deltaT,
+    const scalarField& pressure
+)
+:
+    scalarField(pressure),
+    deltaT_(deltaT)
+{}
+
+
+Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
+:
+    scalarField(),
+    deltaT_(0.0)
+{
+    // Construct pressure data file
+    IFstream pFile(pFileName);
+
+    // Check pFile stream is OK
+    if (!pFile.good())
+    {
+        FatalErrorInFunction
+            << "Cannot read file " << pFileName
+            << exit(FatalError);
+    }
+
+    if (skip)
+    {
+        scalar dummyt, dummyp;
+
+        for (label i = 0; i < skip; i++)
+        {
+            pFile >> dummyt;
+
+            if (!pFile.good() || pFile.eof())
+            {
+                FatalErrorInFunction
+                    << "Number of points in file " << pFileName
+                    << " is less than the number to be skipped = " << skip
+                    << exit(FatalError);
+            }
+
+            pFile >> dummyp;
+        }
+    }
+
+    scalar t = 0, T0 = 0, T1 = 0;
+    DynamicList<scalar> pData(100000);
+    label i = 0;
+
+    while (!(pFile >> t).eof())
+    {
+        if (i == 0)
+        {
+            T0 = t;
+        }
+
+        T1 = t;
+        pFile >> pData(i++);
+    }
+
+    // Note: assumes fixed time step
+    deltaT_ = (T1 - T0)/pData.size();
+
+    this->transfer(pData);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::graph Foam::noiseFFT::pt() const
+{
+    scalarField t(size());
+    forAll(t, i)
+    {
+        t[i] = i*deltaT_;
+    }
+
+    return graph
+    (
+        "p(t)",
+        "t [s]",
+        "p(t) [Pa]",
+        t,
+        *this
+    );
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
+(
+    const tmp<scalarField>& tpn
+) const
+{
+    tmp<scalarField> tPn2
+    (
+        mag
+        (
+            fft::reverseTransform
+            (
+                ReComplexField(tpn),
+                labelList(1, tpn().size())
+            )
+        )
+    );
+
+    tpn.clear();
+
+    tmp<scalarField> tPn
+    (
+        new scalarField
+        (
+            scalarField::subField(tPn2(), tPn2().size()/2)
+        )
+    );
+    scalarField& Pn = tPn();
+
+    Pn *= 2.0/sqrt(scalar(tPn2().size()));
+    Pn[0] /= 2.0;
+
+    return tPn;
+}
+
+
+Foam::graph Foam::noiseFFT::meanPf(const windowModel& window) const
+{
+    const label N = window.nSamples();
+    const label nWindow = window.nWindow();
+
+    scalarField meanPf(N/2, 0.0);
+
+    for (label windowI = 0; windowI < nWindow; ++windowI)
+    {
+        meanPf += Pf(window.apply<scalar>(*this, windowI));
+    }
+
+    meanPf /= scalar(nWindow);
+
+    scalar deltaf = 1.0/(N*deltaT_);
+    scalarField f(meanPf.size());
+    forAll(f, i)
+    {
+        f[i] = i*deltaf;
+    }
+
+    return graph
+    (
+        "P(f)",
+        "f [Hz]",
+        "P(f) [Pa]",
+        f,
+        meanPf
+    );
+}
+
+
+Foam::graph Foam::noiseFFT::RMSmeanPf(const windowModel& window) const
+{
+    const label N = window.nSamples();
+    const label nWindow = window.nWindow();
+
+    scalarField RMSMeanPf(N/2, 0.0);
+
+    for (label windowI = 0; windowI < nWindow; ++windowI)
+    {
+        RMSMeanPf += sqr(Pf(window.apply<scalar>(*this, windowI)));
+    }
+
+    RMSMeanPf = sqrt(RMSMeanPf/scalar(nWindow));
+
+    scalar deltaf = 1.0/(N*deltaT_);
+    scalarField f(RMSMeanPf.size());
+    forAll(f, i)
+    {
+        f[i] = i*deltaf;
+    }
+
+    return graph
+    (
+        "P(f)",
+        "f [Hz]",
+        "P(f) [Pa]",
+        f,
+        RMSMeanPf
+    );
+}
+
+
+Foam::graph Foam::noiseFFT::PSDf(const windowModel& window) const
+{
+    const label N = window.nSamples();
+    const label nWindow = window.nWindow();
+
+    scalarField psd(N/2, 0.0);
+
+    for (label windowI = 0; windowI < nWindow; ++windowI)
+    {
+        psd += 0.5*sqr(Pf(window.apply<scalar>(*this, windowI)));
+    }
+
+    scalar deltaf = 1.0/(N*deltaT_);
+
+    psd /= nWindow*deltaf;
+
+    scalarField f(psd.size());
+
+    forAll(f, i)
+    {
+        f[i] = i*deltaf;
+    }
+
+    return graph
+    (
+        "PSD(f)",
+        "f [Hz]",
+        "PSD(f) [PaPa_Hz]",
+        f,
+        psd
+    );
+}
+
+
+Foam::graph Foam::noiseFFT::PSD(const graph& gPSD) const
+{
+    return graph
+    (
+        "PSD(dB)",
+        "f [Hz]",
+        "PSD_dB(f) [dB]",
+        gPSD.x(),
+        10*log10(gPSD.y()/sqr(p0))
+    );
+}
+
+
+Foam::graph Foam::noiseFFT::Lf(const graph& gPf) const
+{
+    return graph
+    (
+        "L(f)",
+        "f [Hz]",
+        "L(f) [dB]",
+        gPf.x(),
+        20*log10(gPf.y()/p0)
+    );
+}
+
+
+Foam::graph Foam::noiseFFT::Ldelta
+(
+    const graph& gLf,
+    const labelList& freqBandIDs
+) const
+{
+    if (freqBandIDs.size() < 2)
+    {
+        WarningInFunction
+            << "Octave frequency bands are not defined "
+            << "- skipping Ldelta calculation"
+            << endl;
+
+        return graph
+        (
+            "Ldelta",
+            "fm [Hz]",
+            "Ldelta [dB]",
+            scalarField(),
+            scalarField()
+        );
+    }
+
+    const scalarField& f = gLf.x();
+    const scalarField& Lf = gLf.y();
+
+    scalarField ldelta(freqBandIDs.size() - 1, 0.0);
+    scalarField fm(freqBandIDs.size() -1, 0.0);
+
+    for (label bandI = 0; bandI < freqBandIDs.size() - 1; bandI++)
+    {
+        label f0 = freqBandIDs[bandI];
+        label f1 = freqBandIDs[bandI+1];
+        fm[bandI] = f[f0];
+
+        if (f0 == f1) continue;
+
+        for (label freqI = f0; freqI < f1; freqI++)
+        {
+            ldelta[bandI] += pow(10, Lf[freqI]/10.0);
+        }
+        ldelta[bandI] = 10*log10(ldelta[bandI]);
+    }
+
+    return graph
+    (
+        "Ldelta",
+        "fm [Hz]",
+        "Ldelta [dB]",
+        fm,
+        ldelta
+    );
+}
+
+
+Foam::graph Foam::noiseFFT::Pdelta
+(
+    const graph& gPf,
+    const labelList& freqBandIDs
+) const
+{
+    if (freqBandIDs.size() < 2)
+    {
+        WarningInFunction
+            << "Octave frequency bands are not defined "
+            << "- skipping Pdelta calculation"
+            << endl;
+
+        return graph
+        (
+            "Pdelta",
+            "fm [Hz]",
+            "Pdelta [dB]",
+            scalarField(),
+            scalarField()
+        );
+    }
+
+    const scalarField& f = gPf.x();
+    const scalarField& Pf = gPf.y();
+
+    scalarField pdelta(freqBandIDs.size() - 1, 0.0);
+    scalarField fm(pdelta.size());
+
+    for (label bandI = 0; bandI < freqBandIDs.size() - 1; bandI++)
+    {
+        label f0 = freqBandIDs[bandI];
+        label f1 = freqBandIDs[bandI+1];
+        fm[bandI] = f[f0];
+
+        if (f0 == f1) continue;
+
+        for (label freqI = f0; freqI < f1; freqI++)
+        {
+            pdelta[bandI] += sqr(Pf[freqI]);
+        }
+        pdelta[bandI] = sqrt((2.0/3.0)*pdelta[bandI]);
+    }
+
+    return graph
+    (
+        "Pdelta",
+        "fm [Hz]",
+        "Pdelta [dB]",
+        fm,
+        pdelta
+    );
+}
+
+
+Foam::scalar Foam::noiseFFT::Lsum(const graph& gLf) const
+{
+    const scalarField& Lf = gLf.y();
+
+    scalar lsum = 0.0;
+
+    forAll(Lf, i)
+    {
+        lsum += pow(10, Lf[i]/10.0);
+    }
+
+    lsum = 10*log10(lsum);
+
+    return lsum;
+}
+
+
+Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
+{
+    return p0*pow(10.0, db/20.0);
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseFFT::dbToPa
+(
+    const tmp<scalarField>& db
+) const
+{
+    return p0*pow(10.0, db/20.0);
+}
+
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/noise/noiseFFT.H b/src/randomProcesses/noise/noiseFFT/noiseFFT.H
similarity index 54%
rename from src/randomProcesses/noise/noiseFFT.H
rename to src/randomProcesses/noise/noiseFFT/noiseFFT.H
index 085702964c..ab580cbcdf 100644
--- a/src/randomProcesses/noise/noiseFFT.H
+++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,11 +25,26 @@ Class
     Foam::noiseFFT
 
 Description
-    FFT of the pressure field
+    Performs FFT of pressure field to generate noise data.
+
+    General functionality:
+    - Pf: fft of the pressure data
+    - meanPf: multi-window mean fft
+    - RMSmeanPf: multi-window RMS mean fft
+    - PSDf: multi-window power spectral density (PSD) in frequency domain
+    - PSD: multi-window power spectral density (PSD) in dB
+    - Lf: narrow-band pressure-fluctuation level (PFL) in frequency domain
+
+    Octave-based data:
+    - Ldelta: PFL spectrum
+    - Pdelta: pressure spectrum
 
 SourceFiles
     noiseFFT.C
 
+SeeAlso
+    windowModel.H
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef noiseFFT_H
@@ -37,6 +52,7 @@ SourceFiles
 
 #include "scalarField.H"
 #include "graph.H"
+#include "windowModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,7 +60,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class noiseFFT Declaration
+                          Class noiseFFT Declaration
 \*---------------------------------------------------------------------------*/
 
 class noiseFFT
@@ -54,7 +70,7 @@ class noiseFFT
     // Private data
 
         //- Time spacing of the raw data
-        scalar deltat_;
+        scalar deltaT_;
 
 
 public:
@@ -81,34 +97,72 @@ public:
 
     // Member Functions
 
+        //- Return the FFT frequencies
+        static tmp<scalarField> frequencies
+        (
+            const label N,
+            const scalar deltaT
+        );
+
+        //- Return a list of the frequency indices wrt f field that
+        //  correspond to the bands limits for a given octave
+        static void octaveFrequenciesIDs
+        (
+            const scalarField& f,
+            const scalar fLower,
+            const scalar fUpper,
+            const scalar octave,
+            labelList& freqBandIDs
+        );
+
+        //- Return the 1/octave octave frequency bounds
+        static tmp<scalarField> octaveFrequencies
+        (
+            const scalarField& f,
+            const scalar fLower,
+            const scalar fUpper,
+            const scalar octave
+        );
+
         //- Return the graph of p(t)
         graph pt() const;
 
-        //- Return the nth window
-        tmp<scalarField> window(const label N, const label n) const;
-
-        //- Return the Hanning window function
-        tmp<scalarField> Hanning(const label N) const;
-
         //- Return the fft of the given pressure data
         tmp<scalarField> Pf(const tmp<scalarField>& pn) const;
 
-        //- Return the multi-window mean fft of the complete pressure data
-        graph meanPf(const label N, const label nw) const;
+        //- Return the multi-window mean fft of the complete pressure data [Pa]
+        graph meanPf(const windowModel& window) const;
 
-        //- Return the multi-window RMS mean fft of the complete pressure data
-        graph RMSmeanPf(const label N, const label nw) const;
+        //- Return the multi-window RMS mean fft of the complete pressure
+        //  data [Pa]
+        graph RMSmeanPf(const windowModel& window) const;
 
-        //- Return the narrow-band PFL (pressure-fluctuation level) spectrum
+        //- Return the multi-window PSD (power spectral density) of the complete
+        //  pressure data [Pa^2/Hz]
+        graph PSDf(const windowModel& window) const;
+
+        //- Return the PSD [dB]
+        graph PSD(const graph& gPSD) const;
+
+        //- Return the narrow-band PFL (pressure-fluctuation level)
+        //  spectrum [dB]
         graph Lf(const graph& gPf) const;
 
-        //- Return the one-third-octave-band PFL spectrum
-        //  starting at octave with mean frequency f1
-        graph Ldelta(const graph& gLf, const scalar f1, const scalar fU) const;
+        //- Return the octave-band PFL spectrum starting at octave
+        //  frequencies given by the supplied frequency bands [dB]
+        graph Ldelta
+        (
+            const graph& gLf,
+            const labelList& freqBandIDs
+        ) const;
 
-        //- Return the one-third-octave-band pressure spectrum
-        //  starting at octave with mean frequency f1
-        graph Pdelta(const graph& gLf, const scalar f1, const scalar fU) const;
+        //- Return the octave-band pressure spectrum at octave
+        //  frequencies given by the supplied frequency bands [dB]
+        graph Pdelta
+        (
+            const graph& gLf, 
+            const labelList& freqBandIDs
+        ) const;
 
         //- Return the total PFL as the sum of Lf over all frequencies
         scalar Lsum(const graph& gLf) const;
diff --git a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C
new file mode 100644
index 0000000000..c32c4ec21d
--- /dev/null
+++ b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C
@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "noiseModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(noiseModel, 0);
+    defineRunTimeSelectionTable(noiseModel, dictionary);
+}
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+Foam::scalar Foam::noiseModel::checkUniformTimeStep
+(
+    const scalarList& times
+) const
+{
+    scalar deltaT = -1.0;
+
+    if (times.size() > 1)
+    {
+        for (label i = 1; i < times.size(); i++)
+        {
+            scalar dT = times[i] - times[i-1];
+
+            if (deltaT < 0)
+            {
+                deltaT = dT;
+            }
+
+            if (mag(deltaT - dT) > SMALL)
+            {
+                FatalErrorInFunction
+                    << "Unable to process data with a variable time step"
+                    << exit(FatalError);
+            }
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unable to create FFT with a single value"
+            << exit(FatalError);
+    }
+
+    return deltaT;
+}
+
+
+Foam::label Foam::noiseModel::findStartTimeIndex
+(
+    const instantList& allTimes,
+    const scalar startTime
+) const
+{
+    forAll(allTimes, timeI)
+    {
+        const instant& i = allTimes[timeI];
+
+        if (i.value() >= startTime)
+        {
+            return timeI;
+        }
+    }
+
+    return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::noiseModel::noiseModel(const dictionary& dict)
+:
+    dict_(dict),
+    pRef_(dict.lookupOrDefault("pRef", 0)),
+    nSamples_(dict.lookupOrDefault("N", 65536)),
+    fLower_(dict.lookupOrDefault("fl", 25)),
+    fUpper_(dict.lookupOrDefault("fu", 10000)),
+    startTime_(dict.lookupOrDefault("startTime", 0)),
+    windowModelPtr_(windowModel::New(dict, nSamples_))
+{
+    // Check number of samples  - must be a power of 2 for our FFT
+    bool powerOf2 = ((nSamples_ != 0) && !(nSamples_ & (nSamples_ - 1)));
+    if (!powerOf2)
+    {
+        FatalIOErrorInFunction(dict)
+            << "N: Number of samples in sampling windows must be a "
+            << "power of 2"
+            << exit(FatalIOError);
+    }
+
+    if (fLower_ < 0)
+    {
+        FatalIOErrorInFunction(dict)
+            << "fl: lower frequency bound must be greater than zero"
+            << exit(FatalIOError);
+
+    }
+
+    if (fUpper_ < 0)
+    {
+        FatalIOErrorInFunction(dict)
+            << "fu: upper frequency bound must be greater than zero"
+            << exit(FatalIOError);
+
+    }
+
+    if (fUpper_ < fLower_)
+    {
+        FatalIOErrorInFunction(dict)
+            << "fu: upper frequency bound must be greater than lower "
+            << "frequency bound (fl)"
+            << exit(FatalIOError);
+
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::noiseModel::~noiseModel()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H
new file mode 100644
index 0000000000..f7e739b7d0
--- /dev/null
+++ b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H
@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::noiseModel
+
+Description
+    Base class for noise models.
+
+    Data is read from a dictionary, e.g.
+
+    \verbatim
+    pRef            0;
+    N               4096;
+    fl              25;
+    fu              25;
+    startTime       0;
+    \endverbatim
+
+    where
+    \table
+        Property    | Description                   | Required  | Default value
+        pRef        | Reference pressure            | no        | 0
+        N           | Number of samples in sampling window | no | 65536 (2^16)
+        fl          | Lower frequency bounds        | no        | 25
+        fu          | Upper frequency bounds        | no        | 10000
+    \endtable
+
+Note
+    The number of samples in the sampling window must be a power of 2
+
+
+SourceFiles
+    noiseModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noiseModel_H
+#define noiseModel_H
+
+#include "dictionary.H"
+#include "scalarList.H"
+#include "instantList.H"
+#include "windowModel.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class noiseModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class noiseModel
+{
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        noiseModel(const noiseModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const noiseModel&);
+
+
+protected:
+
+    // Protected Data
+
+        //- Copy of dictionary used for construction
+        const dictionary dict_;
+
+        //- Reference pressure
+        scalar pRef_;
+
+        //- Number of samples in sampling window, default = 2^16
+        label nSamples_;
+
+        //- Lower frequency limit, default = 25Hz
+        scalar fLower_;
+
+        //- Upper frequency limit, default = 10kHz
+        scalar fUpper_;
+
+        //- Start time, default = 0s
+        scalar startTime_;
+
+        //- Window model
+        autoPtr<windowModel> windowModelPtr_;
+        
+
+
+    // Protected Member Functions
+
+        //- Check and return uniform time step
+        scalar checkUniformTimeStep
+        (
+            const scalarList& times
+        ) const; 
+
+        //- Find and return start time index
+        label findStartTimeIndex
+        (
+            const instantList& allTimes,
+            const scalar startTime
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("noiseModel");
+
+    //- Run time selection table  
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        noiseModel,
+        dictionary,
+        (
+            const dictionary& dict
+        ),
+        (dict)
+    );
+
+    //- Selector
+    static autoPtr<noiseModel> New(const dictionary& dict);
+
+    //- Constructor
+    noiseModel(const dictionary& dict);
+
+    //- Destructor
+    virtual ~noiseModel();
+
+
+    // Public Member Functions
+
+        //- Abstract call to calculate
+        virtual void calculate() = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModelNew.C b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModelNew.C
new file mode 100644
index 0000000000..3f5fa4f25f
--- /dev/null
+++ b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModelNew.C
@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "noiseModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::noiseModel> Foam::noiseModel::New(const dictionary& dict)
+{
+    const word modelType(dict.lookup("noiseModel"));
+
+    Info<< "Selecting noiseModel " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown noiseModel type "
+            << modelType << nl << nl
+            << "Valid noiseModel types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<noiseModel>(cstrIter()(dict.subDict(modelType + "Coeffs")));
+}
+
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C
new file mode 100644
index 0000000000..2014d721e3
--- /dev/null
+++ b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C
@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointNoise.H"
+#include "noiseFFT.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace noiseModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(pointNoise, 0);
+addToRunTimeSelectionTable(noiseModel, pointNoise, dictionary);
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void pointNoise::filterTimeData
+(
+    const CSV<scalar>& pData,
+    scalarField& t,
+    scalarField& p
+)
+{
+    const scalarField t0(pData.x());
+    const scalarField p0(pData.y());
+
+    DynamicList<scalar> tf(t0.size());
+    DynamicList<scalar> pf(t0.size());
+
+    forAll(t0, timeI)
+    {
+        if (t0[timeI] >= startTime_)
+        {
+            tf.append(t0[timeI]);
+            pf.append(p0[timeI]);
+        }
+    }
+
+    t.transfer(tf);
+    p.transfer(pf);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void pointNoise::calculate()
+{
+    // Point data only handled by master
+    if (!Pstream::master())
+    {
+        return;
+    }
+
+    Info<< "Reading data file" << endl;
+
+    CSV<scalar> pData("pressure", dict_, "Data");
+    fileName baseFileName(pData.fName().lessExt());
+
+    // Time and pressure history data
+    scalarField t, p;
+    filterTimeData(pData, t, p);
+    Info<< "    read " << t.size() << " values" << nl << endl;
+
+    Info<< "Creating noise FFT" << endl;
+    const scalar deltaT = checkUniformTimeStep(t);
+
+    // Determine the windowing
+    windowModelPtr_->validate(t.size());
+
+    // Create the fft
+    noiseFFT nfft(deltaT, p);
+
+    nfft -= pRef_;
+
+    graph Pf(nfft.RMSmeanPf(windowModelPtr_()));
+    Info<< "    Creating graph for " << Pf.title() << endl;
+    Pf.write(baseFileName + graph::wordify(Pf.title()), graphFormat_);
+
+    graph Lf(nfft.Lf(Pf));
+    Info<< "    Creating graph for " << Lf.title() << endl;
+    Lf.write(baseFileName + graph::wordify(Lf.title()), graphFormat_);
+
+    graph PSDf(nfft.PSDf(windowModelPtr_()));
+    Info<< "    Creating graph for " << PSDf.title() << endl;
+    PSDf.write(baseFileName + graph::wordify(PSDf.title()), graphFormat_);
+
+    graph PSD(nfft.PSD(PSDf));
+    Info<< "    Creating graph for " << PSD.title() << endl;
+    PSD.write(baseFileName + graph::wordify(PSD.title()), graphFormat_);
+
+    labelList octave13BandIDs;
+    noiseFFT::octaveFrequenciesIDs(Pf.x(), fLower_, fUpper_, 3, octave13BandIDs);
+
+    graph Ldelta(nfft.Ldelta(Lf, octave13BandIDs));
+    Info<< "    Creating graph for " << Ldelta.title() << endl;
+    Ldelta.write(baseFileName + graph::wordify(Ldelta.title()), graphFormat_);
+
+    graph Pdelta(nfft.Pdelta(Pf, octave13BandIDs));
+    Info<< "    Creating graph for " << Pdelta.title() << endl;
+    Pdelta.write(baseFileName + graph::wordify(Pdelta.title()), graphFormat_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+pointNoise::pointNoise(const dictionary& dict)
+:
+    noiseModel(dict),
+    graphFormat_(dict.lookupOrDefault<word>("graphFormat", "raw"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+pointNoise::~pointNoise()
+{}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace noiseModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H
new file mode 100644
index 0000000000..020e462238
--- /dev/null
+++ b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H
@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::noiseModels::pointNoise
+
+Description
+    Perform noise analysis on point-based pressure data.
+
+    Input data is read from a dictionary, e.g.
+
+    \verbatim
+    // Pressure reference
+    pRef            0;
+
+    // Number of samples in sampling window
+    // Must be a power of 2, default = 2^16 (=65536)
+    N               4096;
+
+    // Lower frequency bounds 
+    fl              25;
+
+    // Upper frequency bounds 
+    fu              25;
+
+    // Start time
+    startTime       0;
+
+    windowModel     <modelType>
+    <modelType>Coeffs
+    {
+        ...
+    }
+
+    // Pressure data supplied in CSV file format
+    csvFileData
+    {
+        fileName        "pressureData";
+        nHeaderLine     1;
+        refColumn       0;
+        componentColumns (1);
+        separator       " ";
+        mergeSeparators yes;
+    }
+
+    graphFormat     raw;
+
+    \endverbatim
+
+SourceFiles
+    pointNoise.C
+
+SeeAlso
+    noiseModel.H
+    windowModel.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointNoise_H
+#define pointNoise_H
+
+#include "noiseModel.H"
+#include "CSV.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace noiseModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class pointNoise Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointNoise
+:
+    public noiseModel
+{
+
+protected:
+
+    // Protected data
+
+        //- Graph format
+        word graphFormat_;
+
+
+    // Protected Member Functions
+
+        void filterTimeData
+        (
+            const CSV<scalar>& pData,
+            scalarField& t,
+            scalarField& p
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("pointNoise");
+
+    //- Constructor
+    pointNoise(const dictionary& dict);
+
+    //- Destructor
+    virtual ~pointNoise();
+
+
+    // Public Member Functions
+
+        //- Calculate
+        virtual void calculate();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace noiseModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C
new file mode 100644
index 0000000000..1dcc319766
--- /dev/null
+++ b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C
@@ -0,0 +1,490 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceNoise.H"
+#include "surfaceReader.H"
+#include "surfaceWriter.H"
+#include "noiseFFT.H"
+#include "graph.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace noiseModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(surfaceNoise, 0);
+addToRunTimeSelectionTable(noiseModel, surfaceNoise, dictionary);
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void surfaceNoise::initialise(const dictionary& dict)
+{
+    // All reading performed on the master processor only
+    if (Pstream::master())
+    {
+        // Create the surface reader
+        const word readerType(dict.lookup("reader"));
+        readerPtr_.reset(surfaceReader::New(readerType, inputFileName_).ptr());
+
+        // Find the index of the pressure data
+        const word pName(dict.lookupOrDefault<word>("pName", "p"));
+        const List<word> fieldNames(readerPtr_->fieldNames(0));
+        pIndex_ = findIndex(fieldNames, pName);
+        if (pIndex_ == -1)
+        {
+            FatalErrorInFunction
+                << "Unable to find pressure field name " << pName
+                << " in list of available fields: " << fieldNames
+                << exit(FatalError);
+        }
+
+        // Create the surface writer
+        // - Could be done later, but since this utility can process a lot of
+        //   data we can ensure that the user-input is correct prior to doing
+        //   the heavy lifting
+        const word writerType(dict.lookup("writer"));
+        dictionary optDict
+        (
+            dict.subOrEmptyDict("writeOptions").subOrEmptyDict(writerType)
+        );
+        writerPtr_.reset(surfaceWriter::New(writerType, optDict).ptr());
+
+        // Set the time range
+        const instantList allTimes = readerPtr_->times();
+        startTimeIndex_ = findStartTimeIndex(allTimes, startTime_);
+
+        // Determine the windowing
+        label nAvailableTimes = allTimes.size() - startTimeIndex_;
+        label nRequiredTimes = windowModelPtr_->validate(nAvailableTimes);
+
+        // Restrict times
+        times_.setSize(nRequiredTimes);
+        forAll(times_, timeI)
+        {
+            times_[timeI] = allTimes[timeI + startTimeIndex_].value();
+        }
+        deltaT_ = checkUniformTimeStep(times_);
+
+        const meshedSurface& surf = readerPtr_->geometry();
+        nFace_ = surf.size();
+    }
+
+    Pstream::scatter(pIndex_);
+    Pstream::scatter(times_);
+    Pstream::scatter(startTimeIndex_);
+    Pstream::scatter(deltaT_);
+    Pstream::scatter(nFace_);
+}
+
+
+void surfaceNoise::readSurfaceData
+(
+    const labelList& procFaceOffset,
+    List<scalarField>& pData
+)
+{
+    // Data is stored as pressure on surface at a given time.  Now we need to
+    // pivot the data so that we have the complete pressure time history per
+    // surface face.  In serial mode, this results in all pressure data being
+    // loaded into memory (!)
+
+    Info << "Reading pressure data" << endl;
+
+    if (Pstream::parRun())
+    {
+        PstreamBuffers pBufs(Pstream::nonBlocking);
+
+        // Procedure:
+        // 1. Master processor reads pressure data for all faces for all times
+        // 2. Master sends each processor a subset of faces
+        // 3. Local processor reconstructs the full time history for the subset
+        //    of faces
+        // Note: reading all data on master to avoid potential NFS problems...
+
+        const label myProcI = Pstream::myProcNo();
+        const label nLocalFace =
+            procFaceOffset[myProcI + 1] - procFaceOffset[myProcI];
+
+        // Complete pressure time history data for subset of faces
+        pData.setSize(nLocalFace);
+        const label nTimes = times_.size();
+        forAll(pData, faceI)
+        {
+            pData[faceI].setSize(nTimes);
+        }
+
+        // Read and send data
+        forAll(times_, i)
+        {
+            pBufs.clear();
+
+            if (Pstream::master())
+            {
+                label timeI = i + startTimeIndex_;
+
+                Info<< "    time: " << times_[i] << endl;
+
+                // Read pressure at all faces for time timeI
+                scalarField p(readerPtr_->field(timeI, pIndex_, scalar(0)));
+
+                // Send subset of faces to each processor
+                for (label procI = 0; procI < Pstream::nProcs(); procI++)
+                {
+                    label face0 = procFaceOffset[procI];
+                    label nLocalFace =
+                        procFaceOffset[procI + 1] - procFaceOffset[procI];
+
+                    UOPstream toProc(procI, pBufs);
+                    toProc << SubList<scalar>(p, nLocalFace, face0);
+                }
+            }
+
+            pBufs.finishedSends();
+
+            // Receive data from the master
+            UIPstream fromProc(0, pBufs);
+
+            scalarList pSlice(fromProc);
+
+            forAll(pSlice, faceI)
+            {
+                pData[faceI][i] = pSlice[faceI];
+            }
+        }
+    }
+    else
+    {
+        const label nLocalFace = procFaceOffset[0];
+
+        pData.setSize(nLocalFace);
+        forAll(times_, timeI)
+        {
+            forAll(pData, faceI)
+            {
+                pData[faceI].setSize(times_.size());
+            }
+        }
+
+        forAll(times_, i)
+        {
+            label timeI = i + startTimeIndex_;
+
+            Info<< "    time: " << times_[i] << endl;
+            const scalarField p(readerPtr_->field(timeI, pIndex_, scalar(0)));
+
+            forAll(p, faceI)
+            {
+                pData[faceI][i] = p[faceI];
+            }
+        }
+    }
+
+    Info<< "Read "
+        << returnReduce(pData.size(), sumOp<label>())
+        << " pressure traces with "
+        << times_.size()
+        << " time values" << nl << endl;
+}
+
+
+void surfaceNoise::writeSurfaceData
+(
+    const word& fName,
+    const word& groupName,
+    const word& title,
+    const scalar freq,
+    const scalarField& data,
+    const labelList& procFaceOffset
+) const
+{
+    Info<< "    processing " << title << " for frequency " << freq << endl;
+
+    fileName outDir
+    (
+        fileName("postProcessing")/"noise"/groupName/Foam::name(freq)
+    );
+
+    if (Pstream::parRun())
+    {
+        // Collect the surface data so that we can output the surfaces
+
+        PstreamBuffers pBufs(Pstream::nonBlocking);
+
+        if (!Pstream::master())
+        {
+            UOPstream toProc(0, pBufs);
+            toProc << data;
+        }
+
+        pBufs.finishedSends();
+
+        if (Pstream::master())
+        {
+            const meshedSurface& surf = readerPtr_->geometry();
+
+            scalarField allData(surf.size());
+
+            forAll(data, faceI)
+            {
+                // Master procFaceOffset is zero...
+                allData[faceI] = data[faceI];
+            }
+
+            for (label procI = 1; procI < Pstream::nProcs(); procI++)
+            {
+                UIPstream fromProc(procI, pBufs);
+                scalarList dataSlice(fromProc);
+                forAll(dataSlice, i)
+                {
+                    label faceI = procFaceOffset[procI] + i;
+                    allData[faceI] = dataSlice[i];
+                }
+            }
+
+            fileName outFileName = writerPtr_->write
+            (
+                outDir,
+                fName,
+                surf.points(),
+                surf.faces(),
+                title,
+                allData,
+                false
+            );
+        }
+    }
+    else
+    {
+        const meshedSurface& surf = readerPtr_->geometry();
+
+        writerPtr_->write
+        (
+            outDir,
+            fName,
+            surf.points(),
+            surf.faces(),
+            title,
+            data,
+            false
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+surfaceNoise::surfaceNoise(const dictionary& dict)
+:
+    noiseModel(dict),
+    inputFileName_(dict.lookup("inputFile")),
+    pIndex_(0),
+    times_(),
+    deltaT_(0),
+    startTimeIndex_(0),
+    nFace_(0),
+    fftWriteInterval_(dict.lookupOrDefault("fftWriteInterval", 1))
+{
+    initialise(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+surfaceNoise::~surfaceNoise()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void surfaceNoise::calculate()
+{
+    // Container for pressure time history data per face
+    List<scalarField> pData;
+
+    // Processor procFaceOffsets
+    labelList procFaceOffset;
+    if (Pstream::parRun())
+    {
+        const label nProcs = Pstream::nProcs();
+        const label nFacePerProc = floor(nFace_/nProcs) + 1;
+
+        procFaceOffset.setSize(nProcs + 1, 0);
+        for (label i = 1; i < procFaceOffset.size(); i++)
+        {
+            procFaceOffset[i] = min(i*nFacePerProc, nFace_);
+        }
+    }
+    else
+    {
+        procFaceOffset.setSize(1, nFace_);
+    }
+
+    // Read pressure data from file
+    readSurfaceData(procFaceOffset, pData);
+
+    // Process the pressure data, and store results as surface values per
+    // frequency so that it can be output using the surface writer
+
+    Info<< "Creating noise FFTs" << endl;
+
+    // Storage for FFT data
+    const label nLocalFace = pData.size();
+    const scalarField freq1(noiseFFT::frequencies(nSamples_, deltaT_));
+    const label nFFT = freq1.size()/fftWriteInterval_;
+    List<scalarField> surfPf(nFFT);
+    List<scalarField> surfLf(nFFT);
+    List<scalarField> surfPSD(nFFT);
+    forAll(surfPf, freqI)
+    {
+        surfPf[freqI].setSize(nLocalFace);
+        surfLf[freqI].setSize(nLocalFace);
+        surfPSD[freqI].setSize(nLocalFace);
+    }
+
+    // Storage for 1/3 octave data
+    labelList octave13BandIDs;
+    noiseFFT::octaveFrequenciesIDs(freq1, fLower_, fUpper_, 3, octave13BandIDs);
+
+    List<scalarField> surfPdelta(octave13BandIDs.size() - 1);
+    List<scalarField> surfLdelta(octave13BandIDs.size() - 1);
+    forAll(surfPdelta, freqI)
+    {
+        surfPdelta[freqI].setSize(nLocalFace);
+        surfLdelta[freqI].setSize(nLocalFace);
+    }
+
+    forAll(pData, faceI)
+    {
+        const scalarField& p = pData[faceI];
+
+        noiseFFT nfft(deltaT_, p);
+
+        nfft -= pRef_;
+        graph Pf(nfft.RMSmeanPf(windowModelPtr_()));
+        graph Lf(nfft.Lf(Pf));
+        graph PSDf(nfft.PSDf(windowModelPtr_()));
+        graph PSD(nfft.PSD(PSDf));
+
+        // Store the frequency results in slot for face of surface
+        forAll(surfPf, i)
+        {
+            label freqI = (i + 1)*fftWriteInterval_ - 1;
+            surfPf[i][faceI] = Pf.y()[freqI];
+            surfLf[i][faceI] = Lf.y()[freqI];
+            surfPSD[i][faceI] = PSD.y()[freqI];
+        }
+
+        graph Pdelta(nfft.Pdelta(Pf, octave13BandIDs));
+        graph Ldelta(nfft.Ldelta(Lf, octave13BandIDs));
+
+        // Store the 1/3 octave results in slot for face of surface
+        forAll(surfPdelta, freqI)
+        {
+            surfPdelta[freqI][faceI] = Pdelta.y()[freqI];
+            surfLdelta[freqI][faceI] = Ldelta.y()[freqI];
+        }
+
+        // Free the storage for p
+//        p.clear();
+    }
+
+    Info<< "Writing fft surface data" << endl;
+
+    forAll(surfPf, i)
+    {
+        label freqI = i*fftWriteInterval_;
+        const word& fName = inputFileName_.name(true);
+        const word gName = "fft";
+        writeSurfaceData
+        (
+            fName,
+            gName,
+            "Pf",
+            freq1[freqI],
+            surfPf[i],
+            procFaceOffset
+        );
+        writeSurfaceData
+        (
+            fName,
+            gName,
+            "Lf",
+            freq1[freqI],
+            surfLf[i],
+            procFaceOffset
+        );
+        writeSurfaceData
+        (
+            fName,
+            gName,
+            "PSD",
+            freq1[freqI],
+            surfPSD[i],
+            procFaceOffset
+        );
+    }
+
+    Info<< "Writing one-third octave surface data" << endl;
+
+    forAll(surfPdelta, i)
+    {
+        const word& fName = inputFileName_.name(true);
+        const word gName = "oneThirdOctave";
+        writeSurfaceData
+        (
+            fName,
+            gName,
+            "Pdelta",
+            octave13BandIDs[i],
+            surfPdelta[i],
+            procFaceOffset
+        );
+        writeSurfaceData
+        (
+            fName,
+            gName,
+            "Ldelta",
+            octave13BandIDs[i],
+            surfLdelta[i],
+            procFaceOffset
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace noiseModels
+} // End namespace Foam
+
+// ************************************************************************* //
+
+
+
diff --git a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H
new file mode 100644
index 0000000000..e2416e838a
--- /dev/null
+++ b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H
@@ -0,0 +1,201 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::noiseModels::surfaceNoise
+
+Description
+    Perform noise analysis on surface-based pressure data.
+
+    Input data is read from a dictionary, e.g.
+
+    \verbatim
+    // Pressure reference
+    pRef            0;
+
+    // Number of samples in sampling window
+    // Must be a power of 2, default = 2^16 (=65536)
+    N               4096;
+
+    // Lower frequency bounds 
+    fl              25;
+
+    // Upper frequency bounds 
+    fu              25;
+
+    // Start time
+    startTime       0;
+
+    windowModel     <modelType>
+    <modelType>Coeffs
+    {
+        ...
+    }
+
+    // Input file
+    inputFile   "postProcessing/faceSource1/surface/patch/patch.case";
+
+    // Write interval for FFT data, default = 1
+    fftWriteInterval 100;
+
+    // Surface reader
+    reader      ensight;
+
+    // Surface writer
+    writer      ensight;
+
+    // Collate times for ensight output - ensures geometry is only written once
+    writeOptions
+    {
+        ensight
+        {
+            collateTimes 1;
+        }
+    }
+    \endverbatim
+
+SourceFiles
+    surfaceNoise.C
+
+SeeAlso
+    noiseModel.H
+
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfaceNoise_H
+#define surfaceNoise_H
+
+#include "noiseModel.H"
+#include "labelList.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class surfaceReader;
+class surfaceWriter;
+
+namespace noiseModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class surfaceNoise Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfaceNoise
+:
+    public noiseModel
+{
+
+protected:
+
+    // Protected Data
+
+        //- Input file name
+        const fileName inputFileName_;
+
+        //- Index of pressure field in reader field list
+        label pIndex_;
+
+        //- Sample times
+        scalarList times_;
+
+        //- Time step (constant)
+        scalar deltaT_;
+
+        //- Start time index
+        label startTimeIndex_;
+
+        //- Number of surface faces
+        label nFace_;
+
+        //- Frequency data output interval, default = 1
+        //  nSamples/2 data points are returned from the FFT, which can 
+        //  result in a very large number of output files (1 per frequency)
+        label fftWriteInterval_;
+     
+        //- Pointer to the surface reader
+        mutable autoPtr<surfaceReader> readerPtr_;
+
+        //- Pointer to the surface writer
+        autoPtr<surfaceWriter> writerPtr_;
+
+
+    // Protected Member Functions
+
+        //- Initialise
+        void initialise(const dictionary& dict);
+
+        //- Read surface data
+        void readSurfaceData
+        (
+            const labelList& procFaceOffset,
+            List<scalarField>& pData
+        );
+
+        //- Write surface data to file
+        void writeSurfaceData
+        (
+            const word& fName,
+            const word& groupName,
+            const word& title,
+            const scalar freq,
+            const scalarField& data,
+            const labelList& procFaceOffset
+        ) const;
+
+
+
+public:
+
+    //- Runtime type information
+    TypeName("surfaceNoise");
+
+    //- Constructor
+    surfaceNoise(const dictionary& dict);
+
+    //- Destructor
+    virtual ~surfaceNoise();
+
+
+    // Public Member Functions
+
+        //- Calculate
+        virtual void calculate();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace noiseModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/windowModels/Hanning/Hanning.C b/src/randomProcesses/windowModels/Hanning/Hanning.C
new file mode 100644
index 0000000000..19a48258ac
--- /dev/null
+++ b/src/randomProcesses/windowModels/Hanning/Hanning.C
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Hanning.H"
+#include "addToRunTimeSelectionTable.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace windowModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Hanning, 0);
+addToRunTimeSelectionTable(windowModel, Hanning, dictionary);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Hanning::Hanning(const dictionary& dict, const label nSamples)
+:
+    windowModel(dict, nSamples),
+    symmetric_(readBool(dict.lookup("symmetric"))),
+    extended_(readBool(dict.lookup("extended"))),
+    alpha_(dict.lookupOrDefault("alpha", 0.5))
+{
+    // Extend range if required
+    label offset = extended_ ? 1 : 0;
+    scalar m = nSamples - 1 + 2*offset;
+    
+    scalarField t(nSamples);
+    forAll(t, i)
+    {
+        t[i] = i + offset;
+    }
+
+    scalarField& wf = *this;
+    wf = (1 - alpha_)*cos(constant::mathematical::twoPi*t/m);
+
+    // Reset second half of window if symmetric
+    if (symmetric_)
+    {
+        label midPointI = 0;
+        if (nSamples % 2 == 0)
+        {
+            midPointI = nSamples/2;
+        }
+        else
+        {
+            midPointI = (nSamples + 1)/2;
+        }
+
+        for (label i = 0; i < midPointI; i++)
+        {
+            wf[nSamples - i - 1] = wf[i];
+        }
+    }
+
+    scalar sumSqr = sum(sqr(wf));
+
+    // Normalisation (if required)
+    wf *= sqrt(nSamples/sumSqr);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Hanning::~Hanning()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Hanning::symmetric() const
+{
+    return symmetric_;
+}
+
+
+bool Hanning::extended() const
+{
+    return extended_;
+}
+
+
+Foam::scalar Hanning::alpha() const
+{
+    return alpha_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace windowModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/windowModels/Hanning/Hanning.H b/src/randomProcesses/windowModels/Hanning/Hanning.H
new file mode 100644
index 0000000000..8a0ca618ed
--- /dev/null
+++ b/src/randomProcesses/windowModels/Hanning/Hanning.H
@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::windowModels::Hanning
+
+Description
+    Hanning window
+
+    The window is described by the function
+    \f[
+        wf = (1 - \alpha) cos(2 \pi t/m);
+    \f]
+
+    Where:
+    \vartable
+        \alpha  | Coefficient with a default value of 0.5
+        t       | time
+        m       | window width
+    \endvartable
+
+    The window can be further manipulated by the controls:
+    - \c symmetric: force the window to be symmetric
+    - \c extended: extend the window by 1 element at start and end to produce
+      non-zero values at the start and end positions.  Note: window is
+      normalised to preserve energy content
+
+SourceFiles
+    Hanning.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Hanning_H
+#define Hanning_H
+
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "windowModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace windowModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Hanning Declaration
+\*---------------------------------------------------------------------------*/
+
+class Hanning
+:
+    public windowModel
+{
+
+protected:
+
+    // Protected Member Data
+
+        //- Symmetric switch
+        bool symmetric_;
+
+        //- Extended switch
+        bool extended_;
+
+        //- Window coefficient, default = 0.5
+        scalar alpha_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Hanning");
+
+
+    //- Construct from dictionary
+    Hanning(const dictionary& dict, const label nSamples);
+
+    //- Destuctor
+    virtual ~Hanning();
+
+
+    // Public Member Functions
+
+        //- Return the symmetric flag
+        bool symmetric() const;
+
+        //- Return the extended flag
+        bool extended() const;
+
+        //- Return the window coefficient
+        scalar alpha() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace windowModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/windowModels/windowModel/windowModel.C b/src/randomProcesses/windowModels/windowModel/windowModel.C
new file mode 100644
index 0000000000..9b0e0658ff
--- /dev/null
+++ b/src/randomProcesses/windowModels/windowModel/windowModel.C
@@ -0,0 +1,100 @@
+#include "windowModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(windowModel, 0);
+    defineRunTimeSelectionTable(windowModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::windowModel::windowModel(const dictionary& dict, const label nSamples)
+:
+    scalarField(nSamples),
+    nOverlapSamples_(0),
+    nWindow_(dict.lookupOrDefault("nWindow", -1))
+{
+    scalar prc = readScalar(dict.lookup("overlapPercent"));
+    nOverlapSamples_ = floor(prc/scalar(100)*nSamples);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::windowModel::~windowModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::label Foam::windowModel::nSamples() const
+{
+    return size();
+}
+
+
+Foam::label Foam::windowModel::nWindow() const
+{
+    return nWindow_;
+}
+
+
+Foam::label Foam::windowModel::nWindowsTotal(label nSamplesTotal) const
+{
+    const label nSamples = this->nSamples();
+
+    return floor((nSamplesTotal - nSamples)/(nSamples - nOverlapSamples_)) + 1;
+}
+
+
+Foam::label Foam::windowModel::validate(const label nSamplesTotal)
+{
+    label nSamples = this->nSamples();
+
+    if (nSamplesTotal < nSamples)
+    {
+        FatalErrorInFunction
+            << "Block size N = " << nSamples
+            << " is larger than total number of data points = " << nSamplesTotal
+            << exit(FatalError);
+    }
+
+    label nWindowAvailable = nWindowsTotal(nSamplesTotal);
+
+    if (nWindow_ == -1)
+    {
+        nWindow_ = nWindowAvailable;
+    }
+
+    if (nWindow_ > nWindowAvailable)
+    {
+        FatalErrorInFunction
+            << "Number of data points calculated with " << nWindow_
+            << " windows greater than the total number of data points"
+            << nl
+            << "    Block size, N = " << nSamples << nl
+            << "    Total number of data points = " << nSamplesTotal << nl
+            << "    Maximum number of windows = " << nWindowAvailable << nl
+            << "    Requested number of windows = " << nWindow_
+            << exit(FatalError);
+    }
+
+    label nRequiredSamples =
+        nWindow_*nSamples - (nWindow_ - 1)*nOverlapSamples_;
+
+    Info<< "Windowing:" << nl
+        << "    Total samples              : " << nSamplesTotal << nl
+        << "    Samples per window         : " << nSamples << nl
+        << "    Number of windows          : " << nWindow_ << nl
+        << "    Overlap size               : " << nOverlapSamples_ << nl
+        << "    Required number of samples : " << nRequiredSamples
+        << endl;
+
+    return nRequiredSamples;
+}
+
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/windowModels/windowModel/windowModel.H b/src/randomProcesses/windowModels/windowModel/windowModel.H
new file mode 100644
index 0000000000..b9e1845fe0
--- /dev/null
+++ b/src/randomProcesses/windowModels/windowModel/windowModel.H
@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::windowModel
+
+Description
+    Base class for windowing models
+
+SourceFiles
+    noiseFFT.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef windowModel_H
+#define windowModel_H
+
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class windowModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class windowModel
+:
+    public scalarField
+{
+
+protected:
+
+    // Protected Data
+
+        //- Number of overlap samples per window
+        label nOverlapSamples_;
+
+        //- Number of windows
+        label nWindow_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("windowModel");
+
+    // Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        windowModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            const label nSamples
+        ),
+        (dict, nSamples)
+    );
+
+
+    //- Construct from dictionary
+    windowModel(const dictionary& dict, const label nSamples);
+
+
+    // Selectors
+
+        //- Return a reference to the selected window model
+        static autoPtr<windowModel> New
+        (
+            const dictionary& dict,
+            const label nSamples
+        );
+
+
+    //- Destuctor
+    virtual ~windowModel();
+
+
+    // Public Member Functions
+
+        //- Return the number of samples in the window
+        label nSamples() const;
+
+        //- Return the number of windows
+        label nWindow() const;
+
+        //- Return the total number of windows for a given number of samples
+        label nWindowsTotal(label nSamplesTotal) const;
+
+        //- Validate that the window is applicable to the data set size, and
+        //  return the number of required data points
+        label validate(label n);
+
+        //- Return the windowed data
+        template<class Type>
+        tmp<Field<Type> > apply
+        (
+            const Field<Type>& fld,
+            const label windowI
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "windowModelTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/windowModels/windowModel/windowModelNew.C b/src/randomProcesses/windowModels/windowModel/windowModelNew.C
new file mode 100644
index 0000000000..68aed5d766
--- /dev/null
+++ b/src/randomProcesses/windowModels/windowModel/windowModelNew.C
@@ -0,0 +1,63 @@
+#/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "windowModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::windowModel> Foam::windowModel::New
+(
+    const dictionary& dict,
+    const label nSamples
+)
+{
+    const word modelType(dict.lookup("windowModel"));
+
+    Info<< "Selecting windowModel " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "windowModel::New(const dictionary&, const label)"
+        )   << "Unknown windowModel type "
+            << modelType << nl << nl
+            << "Valid windowModel types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<windowModel>
+    (
+        cstrIter()(dict.subDict(modelType + "Coeffs"), nSamples)
+    );
+
+}
+
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C b/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C
new file mode 100644
index 0000000000..1343b55bb9
--- /dev/null
+++ b/src/randomProcesses/windowModels/windowModel/windowModelTemplates.C
@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::windowModel::apply
+(
+    const Field<Type>& fld,
+    const label windowI
+) const
+{
+    label nSamples = this->nSamples();
+
+    if (nSamples > fld.size())
+    {
+        FatalErrorIn
+        (
+            "template<class Type> "
+            "Foam::tmp<Foam::Field<Type> > Foam::windowModel::apply"
+            "("
+                "const Field<Type>&, "
+                "const label"
+            ") const"
+        )
+            << "Number of samples in sampling window is greater than the "
+            << "size of the input field" << nl
+            << "    input field size       = " << fld.size() << nl
+            << "    window size            = " << nSamples << nl
+            << "    requested window index = " << windowI
+            << exit(FatalError);
+    }
+
+
+    tmp<Field<Type> > tresult(new Field<Type>(nSamples, pTraits<Type>::zero));
+    Field<Type>& result = tresult();
+
+    label nWindow = nWindowsTotal(fld.size());
+    if (windowI >= nWindow)
+    {
+        FatalErrorIn
+        (
+            "template<class Type> "
+            "Foam::tmp<Foam::Field<Type> > Foam::windowModel::apply"
+            "("
+                "const Field<Type>&, "
+                "const label"
+            ") const"
+        )
+            << "Requested window " << windowI << " outside of range. "
+            << "Number of available windows is " << nWindow
+            << abort(FatalError);
+    }
+
+    label windowOffset = windowI*(nSamples - nOverlapSamples_);
+
+    const scalarField& wf = *this;
+    result = wf*SubField<Type>(fld, nSamples, windowOffset);
+
+    return tresult;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 03fab92a74..3416dd38d1 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -51,6 +51,12 @@ $(surfWriters)/starcd/starcdSurfaceWriter.C
 $(surfWriters)/vtk/vtkSurfaceWriter.C
 $(surfWriters)/boundaryData/boundaryDataSurfaceWriter.C
 
+surfReaders = sampledSurface/readers
+
+$(surfReaders)/surfaceReader.C
+$(surfReaders)/surfaceReaderNew.C
+$(surfReaders)/ensight/ensightSurfaceReader.C
+
 graphField/writePatchGraph.C
 graphField/writeCellGraph.C
 graphField/makeGraph.C
diff --git a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C
new file mode 100644
index 0000000000..6d17bb5aa0
--- /dev/null
+++ b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C
@@ -0,0 +1,416 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "ensightSurfaceReader.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(ensightSurfaceReader, 0);
+    addToRunTimeSelectionTable(surfaceReader, ensightSurfaceReader, fileName);
+}
+
+
+void Foam::ensightSurfaceReader::skip(const label n, IFstream& is) const
+{
+    label i = 0;
+    token t;
+    while (is.good() && (i < n))
+    {
+        is >> t;
+        i++;
+
+        if (debug)
+        {
+            Info<< "Skipping token " << t << endl;
+        }
+    }
+
+    if (i != n)
+    {
+        WarningInFunction
+            << "Requested to skip " << n << "tokens, but stream exited after "
+            << i << " tokens. Last token read: " << t
+            << endl;
+    }
+}
+
+
+void Foam::ensightSurfaceReader::debugSection
+(
+    const word& expected,
+    IFstream& is
+) const
+{
+    word actual(is);
+
+    if (expected != actual)
+    {
+        FatalIOErrorInFunction(is)
+            << "Expected section header '" << expected
+            << "' but read the word '" << actual << "'"
+            << exit(FatalIOError);
+    }
+}
+
+
+void Foam::ensightSurfaceReader::readCase(IFstream& is)
+{
+    if (debug)
+    {
+        InfoInFunction<< endl;
+    }
+
+    if (!is.good())
+    {
+        FatalErrorInFunction
+            << "Cannot read file " << is.name()
+            << exit(FatalError);
+    }
+
+    // Read the file
+    debugSection("FORMAT", is);
+    skip(3, is); // type: ensight gold
+
+    debugSection("GEOMETRY", is);
+    readSkip(is, 2, meshFileName_);
+
+    debugSection("VARIABLE", is);
+
+    DynamicList<word> fieldNames(10);
+    DynamicList<word> fieldFileNames(10);
+    word fieldName;
+    word fieldFileName;
+    while (is.good())
+    {
+        word primitiveType(is); // scalar, vector
+
+        if (primitiveType == "TIME")
+        {
+            break;
+        }
+
+        readSkip(is, 3, fieldName); // p, U etc
+        fieldNames.append(fieldName);
+
+        is >> fieldFileName; // surfaceName.****.fieldName
+        fieldFileNames.append(fieldFileName);
+    }
+    fieldNames_.transfer(fieldNames);
+    fieldFileNames_.transfer(fieldFileNames);
+
+    if (debug)
+    {
+        Info<< "fieldNames: " << fieldNames << nl
+            << "fieldFileNames: " << fieldFileNames << endl;
+    }
+
+    // Start reading time information
+    skip(3, is); // time set: 1
+    readSkip(is, 3, nTimeSteps_);
+    readSkip(is, 3, timeStartIndex_);
+    readSkip(is, 2, timeIncrement_);
+
+    if (debug)
+    {
+        Info<< "nTimeSteps: " << nTimeSteps_ << nl
+            << "timeStartIndex: " << timeStartIndex_ << nl
+            << "timeIncrement: " << timeIncrement_ << endl;
+    }
+
+    // Read the time values
+    skip(2, is);  // time values:
+    timeValues_.setSize(nTimeSteps_);
+    for (label i = 0; i < nTimeSteps_; i++)
+    {
+        scalar t(readScalar(is));
+
+        timeValues_[i].value() = t;
+        // TODO: use character representation of t directly instead of
+        // regenerating from scalar value
+        timeValues_[i].name() = Foam::name(t);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::ensightSurfaceReader::ensightSurfaceReader(const fileName& fName)
+:
+    surfaceReader(fName),
+    baseDir_(fName.path()),
+    meshFileName_(),
+    fieldNames_(),
+    fieldFileNames_(),
+    nTimeSteps_(0),
+    timeStartIndex_(0),
+    timeIncrement_(1),
+    timeValues_(),
+    surfPtr_(NULL)
+{
+    IFstream is(fName);
+    readCase(is);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::ensightSurfaceReader::~ensightSurfaceReader()
+{}
+
+
+// * * * * * * * * * * * * * Public Member Functions   * * * * * * * * * * * //
+
+const Foam::meshedSurface& Foam::ensightSurfaceReader::geometry()
+{
+    if (debug)
+    {
+        InfoInFunction<< endl;
+    }
+
+    if (!surfPtr_.valid())
+    {
+        IFstream is(baseDir_/meshFileName_);
+
+        if (!is.good())
+        {
+            FatalErrorInFunction
+                << "Cannot read file " << is.name()
+                << exit(FatalError);
+        }
+
+        token t;
+        while (is.good())
+        {
+            is >> t;
+
+            if (t.isWord())
+            {
+                word wordToken = t.wordToken();
+                if (wordToken == "coordinates")
+                {
+                    break;
+                }
+            }
+        }
+
+        label nPoints(readLabel(is));
+
+        if (debug)
+        {
+            Info<< "nPoints: " << nPoints << endl;
+        }
+
+        pointField points(nPoints);
+        {
+            scalarField x(nPoints);
+            for (label dir = 0; dir < 3; dir++)
+            {
+                forAll(points, pointI)
+                {
+                    x[pointI] = readScalar(is);
+                }
+
+                points.replace(dir, x);
+            }
+        }
+
+
+        // Read faces - may be a mix of tris, quads and polys
+        DynamicList<face> faces(ceil(nPoints/3));
+
+        while (is.good())
+        {
+            token t(is);
+
+            if (is.eof())
+            {
+                break;
+            }
+
+            word faceType(t.wordToken());
+
+            if (debug)
+            {
+                Info<< "faceType: " << faceType << endl;
+            }
+
+            label nFace(readLabel(is));
+
+            if (debug)
+            {
+                Info<< "nFace: " << nFace << endl;
+            }
+
+            if (faceType == "tria3")
+            {
+                label np = 3;
+                for (label faceI = 0; faceI < nFace; faceI++)
+                {
+                    face f(np);
+                    for (label fpI = 0; fpI < np; fpI++)
+                    {
+                        f[fpI] = readLabel(is);
+                    }
+
+                    faces.append(f);
+                }
+            }
+            else if (faceType == "quad4")
+            {
+                label np = 4;
+                for (label faceI = 0; faceI < nFace; faceI++)
+                {
+                    face f(np);
+                    for (label fpI = 0; fpI < np; fpI++)
+                    {
+                        f[fpI] = readLabel(is);
+                    }
+
+                    faces.append(f);
+                }
+            }
+            else if (faceType == "nsided")
+            {
+                labelList np(nFace);
+                for (label faceI = 0; faceI < nFace; faceI++)
+                {
+                    np[faceI] = readLabel(is);
+                }
+                for (label faceI = 0; faceI < nFace; faceI++)
+                {
+                    face f(np[faceI]);
+                    for (label fpI = 0; fpI < f.size(); fpI++)
+                    {
+                        f[fpI] = readLabel(is);
+                    }
+
+                    faces.append(f);
+                }
+            }
+            else if (faceType != "")
+            {
+                WarningInFunction
+                    << "Unknown face type: " << faceType
+                    << ".  Aborting read and continuing with current elements "
+                    << "only" << endl;
+            }
+        }
+
+        if (debug)
+        {
+            Info<< "read nFaces: " << faces.size() << endl;
+        }
+
+        forAll(faces, faceI)
+        {
+            face& f = faces[faceI];
+
+            forAll(f, fpI)
+            {
+                f[fpI]--;
+            }
+        }
+
+        surfPtr_.reset(new meshedSurface(xferMove(points), faces.xfer()));
+    }
+
+    return surfPtr_();
+}
+
+
+Foam::instantList Foam::ensightSurfaceReader::times() const
+{
+    return timeValues_;
+}
+
+
+Foam::wordList Foam::ensightSurfaceReader::fieldNames
+(
+    const label timeIndex
+) const
+{
+    return fieldNames_;
+}
+
+
+Foam::tmp<Foam::Field<Foam::scalar> > Foam::ensightSurfaceReader::field
+(
+    const label timeIndex,
+    const label fieldIndex,
+    const scalar& refValue
+) const
+{
+    return readField<scalar>(timeIndex, fieldIndex);
+}
+
+
+Foam::tmp<Foam::Field<Foam::vector> > Foam::ensightSurfaceReader::field
+(
+    const label timeIndex,
+    const label fieldIndex,
+    const vector& refValue
+) const
+{
+    return readField<vector>(timeIndex, fieldIndex);
+}
+
+
+Foam::tmp<Foam::Field<Foam::sphericalTensor> >
+Foam::ensightSurfaceReader::field
+(
+    const label timeIndex,
+    const label fieldIndex,
+    const sphericalTensor& refValue
+) const
+{
+    return readField<sphericalTensor>(timeIndex, fieldIndex);
+}
+
+
+Foam::tmp<Foam::Field<Foam::symmTensor> > Foam::ensightSurfaceReader::field
+(
+    const label timeIndex,
+    const label fieldIndex,
+    const symmTensor& refValue
+) const
+{
+    return readField<symmTensor>(timeIndex, fieldIndex);
+}
+
+
+Foam::tmp<Foam::Field<Foam::tensor> > Foam::ensightSurfaceReader::field
+(
+    const label timeIndex,
+    const label fieldIndex,
+    const tensor& refValue
+) const
+{
+    return readField<tensor>(timeIndex, fieldIndex);
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H
new file mode 100644
index 0000000000..e0b9afd453
--- /dev/null
+++ b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H
@@ -0,0 +1,191 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::ensightensightSurfaceReader
+
+Description
+    Ensight format surface reader
+
+SourceFiles
+    ensightSurfaceReader.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ensightSurfaceReader_H
+#define ensightSurfaceReader_H
+
+#include "surfaceReader.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class ensightSurfaceReader Declaration
+\*---------------------------------------------------------------------------*/
+
+class ensightSurfaceReader
+:
+    public surfaceReader
+{
+protected:
+
+    // Protected Data
+
+        //- Base directory
+        fileName baseDir_;
+
+        //- Name of mesh file
+        word meshFileName_;
+
+        //- Field names
+        List<word> fieldNames_;
+
+        //- Field file names
+        List<word> fieldFileNames_;
+
+        //- Number of time steps
+        label nTimeSteps_;
+
+        //- Start time index
+        label timeStartIndex_;
+
+        //- Time increment
+        label timeIncrement_;
+
+        //- Times
+        instantList timeValues_;
+
+        //- Pointer to the surface
+        autoPtr<meshedSurface> surfPtr_;
+
+
+    // Protected Member Functions
+
+        //- Read and check a section header
+        void debugSection(const word& expected, IFstream& is) const;
+
+        //- Read the case file
+        void readCase(IFstream& is);
+
+        //- Helper function to skip forward n steps in stream
+        void skip(const label n, IFstream& is) const;
+
+        //- Helper function to return Type after skipping n tokens
+        template<class Type>
+        void readSkip(IFstream& is, const label nSkip, Type& value) const;
+
+        //- Helper function to return a field
+        template<class Type>
+        tmp<Field<Type> > readField
+        (
+            const label timeIndex,
+            const label fieldIndex
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ensight");
+
+    // Constructors
+
+        //- Construct from fileName
+        ensightSurfaceReader(const fileName& fName);
+
+
+    //- Destructor
+    virtual ~ensightSurfaceReader();
+
+
+    // Member Functions
+
+        //- Return a reference to the surface geometry
+        virtual const meshedSurface& geometry();
+
+        //- Return a list of the available times
+        virtual instantList times() const;
+
+        //- Return a list of the available fields at a given time
+        virtual wordList fieldNames(const label timeIndex) const;
+
+        //- Return a scalar field at a given time
+        virtual tmp<Field<scalar> > field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const scalar& refValue = pTraits<scalar>::zero
+        ) const;
+
+        //- Return a scalar field at a given time
+        virtual tmp<Field<vector> > field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const vector& refValue = pTraits<vector>::zero
+        ) const;
+
+        //- Return a sphericalTensor field at a given time
+        virtual tmp<Field<sphericalTensor> > field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const sphericalTensor& reValue = pTraits<sphericalTensor>::zero
+        ) const;
+
+        //- Return a symmTensor field at a given time
+        virtual tmp<Field<symmTensor> > field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const symmTensor& reValue = pTraits<symmTensor>::zero
+        ) const;
+
+        //- Return a tensor field at a given time
+        virtual tmp<Field<tensor> > field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const tensor& reValue = pTraits<tensor>::zero
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ensightSurfaceReaderTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C
new file mode 100644
index 0000000000..edb89955a5
--- /dev/null
+++ b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C
@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include <iomanip>
+#include <sstream>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::ensightSurfaceReader::readSkip
+(
+    IFstream& is,
+    const label nSkip,
+    Type& value
+) const
+{
+    skip(nSkip, is);
+
+    is  >> value;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::ensightSurfaceReader::readField
+(
+    const label timeIndex,
+    const label fieldIndex
+) const
+{
+    if (debug)
+    {
+        InfoInFunction<< endl;
+    }
+
+    const word& fieldName(fieldNames_[fieldIndex]);
+    const label fileIndex = timeStartIndex_ + timeIndex*timeIncrement_;
+
+    fileName fieldFileName(fieldFileNames_[fieldIndex]);
+
+    std::ostringstream oss;
+    oss << std::setfill('0') << std::setw(4) << fileIndex;
+    const word indexStr = oss.str();
+    fieldFileName.replace("****", indexStr);
+
+    IFstream is(baseDir_/fieldFileName);
+
+    if (!is.good())
+    {
+        FatalErrorInFunction
+            << "Cannot read file " << is.name()
+            << " for field " << fieldName
+            << exit(FatalError);
+    }
+
+    // Check that data type is as expected
+    word primitiveType(is);
+
+    if (debug)
+    {
+        Info<< "primitiveType: " << primitiveType << endl;
+    }
+
+    if (primitiveType != pTraits<Type>::typeName)
+    {
+        FatalIOErrorInFunction(is)
+            << "Expected " << pTraits<Type>::typeName << "values "
+            << "but found type " << primitiveType
+            << exit(FatalIOError);
+    }
+
+    tmp<Field<Type> > tField(new Field<Type>());
+
+    label n;
+    if (surfPtr_.valid())
+    {
+        n = surfPtr_->size();
+    }
+    else
+    {
+        n = 1000;
+    }
+
+    Type value;
+    word wValue;
+    label iValue;
+
+    // Read header info: part index, e.g. part 1
+    is  >> wValue >> iValue;
+
+    // Read data file
+    // - Assume that file contains a mix of words and numbers, and that all
+    //   numbers relate to face values, e.g. header comprises of words and
+    //   element types are also words, e.g. tria3, quad4, nsided
+    DynamicList<Type> values(n);
+    while (is.good())
+    {
+        token t(is);
+
+        if (is.eof())
+        {
+            break;
+        }
+
+        if (t.isWord())
+        {
+            wValue = t.wordToken();
+        }
+        else
+        {
+            is.putBack(t);
+            is  >> value;
+            values.append(value);
+        }
+    }
+
+    tField().transfer(values);
+
+    return tField;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/readers/surfaceReader.C b/src/sampling/sampledSurface/readers/surfaceReader.C
new file mode 100644
index 0000000000..fb90c1d80f
--- /dev/null
+++ b/src/sampling/sampledSurface/readers/surfaceReader.C
@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceReader.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(surfaceReader, 0);
+    defineRunTimeSelectionTable(surfaceReader, fileName);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceReader::surfaceReader(const fileName& fName)
+:
+    fileName_(fName)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfaceReader::~surfaceReader()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/readers/surfaceReader.H b/src/sampling/sampledSurface/readers/surfaceReader.H
new file mode 100644
index 0000000000..2ccdc470de
--- /dev/null
+++ b/src/sampling/sampledSurface/readers/surfaceReader.H
@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::surfaceReader
+
+Description
+    Base class for surface readers
+
+SourceFiles
+    surfaceReader.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfaceReader_H
+#define surfaceReader_H
+
+#include "typeInfo.H"
+#include "autoPtr.H"
+#include "MeshedSurfaces.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class surfaceReader Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfaceReader
+{
+protected:
+
+    //- File name
+    fileName fileName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("surfaceReader");
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            surfaceReader,
+            fileName,
+            (
+                const fileName& fName
+            ),
+            (fName)
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected surfaceReader
+        static autoPtr<surfaceReader> New
+        (
+            const word& readType,
+            const fileName& fName
+        );
+
+
+    // Constructors
+
+        //- Construct from fileName
+        surfaceReader(const fileName& fName);
+
+
+    //- Destructor
+    virtual ~surfaceReader();
+
+
+    // Member Functions
+
+        //- Return a reference to the surface geometry
+        virtual const meshedSurface& geometry() = 0;
+
+        //- Return a list of the available times
+        virtual instantList times() const = 0;
+
+        //- Return a list of the available fields at a given time
+        virtual wordList fieldNames(const label timeIndex) const = 0;
+
+        //- Return a scalar field at a given time
+        virtual tmp<Field<scalar> > field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const scalar& refValue = pTraits<scalar>::zero
+        ) const = 0;
+
+        //- Return a vector field at a given time
+        virtual tmp<Field<vector> > field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const vector& refValue = pTraits<vector>::zero
+        ) const = 0;
+
+        //- Return a sphericalTensor field at a given time
+        virtual tmp<Field<sphericalTensor> > field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const sphericalTensor& reValue = pTraits<sphericalTensor>::zero
+        ) const = 0;
+
+        //- Return a symmTensor field at a given time
+        virtual tmp<Field<symmTensor> > field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const symmTensor& reValue = pTraits<symmTensor>::zero
+        ) const = 0;
+
+        //- Return a tensor field at a given time
+        virtual tmp<Field<tensor> > field
+        (
+            const label timeIndex,
+            const label fieldIndex,
+            const tensor& reValue = pTraits<tensor>::zero
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/readers/surfaceReaderNew.C b/src/sampling/sampledSurface/readers/surfaceReaderNew.C
new file mode 100644
index 0000000000..415179e273
--- /dev/null
+++ b/src/sampling/sampledSurface/readers/surfaceReaderNew.C
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceReader.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::surfaceReader> Foam::surfaceReader::New
+(
+    const word& readerType,
+    const fileName& fName
+)
+{
+    fileNameConstructorTable::iterator cstrIter =
+        fileNameConstructorTablePtr_->find(readerType);
+
+    if (cstrIter == fileNameConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown reader type \"" << readerType << "\"\n\n"
+            << "Valid reader types: "
+            << fileNameConstructorTablePtr_->sortedToc() << nl
+            << exit(FatalError);
+    }
+
+    return autoPtr<surfaceReader>(cstrIter()(fName));
+}
+
+
+// ************************************************************************* //