diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
index db2d4140700320b205c91b2c77a538a17314b15d..29fc1f86762e5ba6a75ef70ad2db2cd119d183c1 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -108,15 +108,18 @@ makeLESModel(Smagorinsky);
 #include "WALE.H"
 makeLESModel(WALE);
 
-#include "dynamicLagrangian.H"
-makeLESModel(dynamicLagrangian);
-
 #include "kEqn.H"
 makeLESModel(kEqn);
 
 #include "dynamicKEqn.H"
 makeLESModel(dynamicKEqn);
 
+#include "dynamicLagrangian.H"
+makeLESModel(dynamicLagrangian);
+
+#include "kOmegaSSTDES.H"
+makeLESModel(kOmegaSSTDES);
+
 #include "SpalartAllmarasDES.H"
 makeLESModel(SpalartAllmarasDES);
 
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C
index 1a3dc38620ba382c81b1f38a4541d59321077e75..44c2f5950043e6438be3a397aaa674f8ba333fbf 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -100,15 +100,18 @@ makeLESModel(Smagorinsky);
 #include "WALE.H"
 makeLESModel(WALE);
 
-#include "dynamicLagrangian.H"
-makeLESModel(dynamicLagrangian);
-
 #include "kEqn.H"
 makeLESModel(kEqn);
 
 #include "dynamicKEqn.H"
 makeLESModel(dynamicKEqn);
 
+#include "dynamicLagrangian.H"
+makeLESModel(dynamicLagrangian);
+
+#include "kOmegaSSTDES.H"
+makeLESModel(kOmegaSSTDES);
+
 #include "SpalartAllmarasDES.H"
 makeLESModel(SpalartAllmarasDES);
 
diff --git a/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C b/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C
new file mode 100644
index 0000000000000000000000000000000000000000..dbe579d0cf4fb302780d0ab27d9c5f33be468a1f
--- /dev/null
+++ b/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C
@@ -0,0 +1,506 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "kOmegaSSTBase.H"
+#include "fvOptions.H"
+#include "bound.H"
+#include "wallDist.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+tmp<volScalarField>
+kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::kOmegaSST::F1
+(
+    const volScalarField& CDkOmega
+) const
+{
+    tmp<volScalarField> CDkOmegaPlus = max
+    (
+        CDkOmega,
+        dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
+    );
+
+    tmp<volScalarField> arg1 = min
+    (
+        min
+        (
+            max
+            (
+                (scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
+                scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
+            ),
+            (4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
+        ),
+        scalar(10)
+    );
+
+    return tanh(pow4(arg1));
+}
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+tmp<volScalarField>
+kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::kOmegaSST::F2() const
+{
+    tmp<volScalarField> arg2 = min
+    (
+        max
+        (
+            (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
+            scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
+        ),
+        scalar(100)
+    );
+
+    return tanh(sqr(arg2));
+}
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+tmp<volScalarField>
+kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::kOmegaSST::F3() const
+{
+    tmp<volScalarField> arg3 = min
+    (
+        150*(this->mu()/this->rho_)/(omega_*sqr(y_)),
+        scalar(10)
+    );
+
+    return 1 - tanh(pow4(arg3));
+}
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+tmp<volScalarField>
+kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::kOmegaSST::F23() const
+{
+    tmp<volScalarField> f23(F2());
+
+    if (F3_)
+    {
+        f23.ref() *= F3();
+    }
+
+    return f23;
+}
+
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+void kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::correctNut
+(
+    const volScalarField& S2,
+    const volScalarField& F2
+)
+{
+    this->nut_ = a1_*k_/max(a1_*omega_, b1_*F2*sqrt(S2));
+    this->nut_.correctBoundaryConditions();
+    fv::options::New(this->mesh_).correct(this->nut_);
+
+    BasicTurbulenceModel::correctNut();
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+void kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::correctNut()
+{
+    correctNut(2*magSqr(symm(fvc::grad(this->U_))), F23());
+}
+
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+tmp<volScalarField> kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::epsilonByk
+(
+    const volScalarField& F1,
+    const volScalarField& F2
+) const
+{
+    return betaStar_*omega_;
+}
+
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+tmp<fvScalarMatrix>
+kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::kSource() const
+{
+    return tmp<fvScalarMatrix>
+    (
+        new fvScalarMatrix
+        (
+            k_,
+            dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
+        )
+    );
+}
+
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+tmp<fvScalarMatrix>
+kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::omegaSource() const
+{
+    return tmp<fvScalarMatrix>
+    (
+        new fvScalarMatrix
+        (
+            omega_,
+            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
+        )
+    );
+}
+
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+tmp<fvScalarMatrix> kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::Qsas
+(
+    const volScalarField& S2,
+    const volScalarField& gamma,
+    const volScalarField& beta
+) const
+{
+    return tmp<fvScalarMatrix>
+    (
+        new fvScalarMatrix
+        (
+            omega_,
+            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::kOmegaSST
+(
+    const word& type,
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+:
+    TurbulenceModel
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    alphaK1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "alphaK1",
+            this->coeffDict_,
+            0.85
+        )
+    ),
+    alphaK2_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "alphaK2",
+            this->coeffDict_,
+            1.0
+        )
+    ),
+    alphaOmega1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "alphaOmega1",
+            this->coeffDict_,
+            0.5
+        )
+    ),
+    alphaOmega2_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "alphaOmega2",
+            this->coeffDict_,
+            0.856
+        )
+    ),
+    gamma1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "gamma1",
+            this->coeffDict_,
+            5.0/9.0
+        )
+    ),
+    gamma2_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "gamma2",
+            this->coeffDict_,
+            0.44
+        )
+    ),
+    beta1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "beta1",
+            this->coeffDict_,
+            0.075
+        )
+    ),
+    beta2_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "beta2",
+            this->coeffDict_,
+            0.0828
+        )
+    ),
+    betaStar_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "betaStar",
+            this->coeffDict_,
+            0.09
+        )
+    ),
+    a1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "a1",
+            this->coeffDict_,
+            0.31
+        )
+    ),
+    b1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "b1",
+            this->coeffDict_,
+            1.0
+        )
+    ),
+    c1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "c1",
+            this->coeffDict_,
+            10.0
+        )
+    ),
+    F3_
+    (
+        Switch::lookupOrAddToDict
+        (
+            "F3",
+            this->coeffDict_,
+            false
+        )
+    ),
+
+    y_(wallDist::New(this->mesh_).y()),
+
+    k_
+    (
+        IOobject
+        (
+            IOobject::groupName("k", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+    omega_
+    (
+        IOobject
+        (
+            IOobject::groupName("omega", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    )
+{
+    bound(k_, this->kMin_);
+    bound(omega_, this->omegaMin_);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+bool kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::read()
+{
+    if (BasicTurbulenceModel::read())
+    {
+        alphaK1_.readIfPresent(this->coeffDict());
+        alphaK2_.readIfPresent(this->coeffDict());
+        alphaOmega1_.readIfPresent(this->coeffDict());
+        alphaOmega2_.readIfPresent(this->coeffDict());
+        gamma1_.readIfPresent(this->coeffDict());
+        gamma2_.readIfPresent(this->coeffDict());
+        beta1_.readIfPresent(this->coeffDict());
+        beta2_.readIfPresent(this->coeffDict());
+        betaStar_.readIfPresent(this->coeffDict());
+        a1_.readIfPresent(this->coeffDict());
+        b1_.readIfPresent(this->coeffDict());
+        c1_.readIfPresent(this->coeffDict());
+        F3_.readIfPresent("F3", this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+void kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::correct()
+{
+    if (!this->turbulence_)
+    {
+        return;
+    }
+
+    // Local references
+    const alphaField& alpha = this->alpha_;
+    const rhoField& rho = this->rho_;
+    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
+    const volVectorField& U = this->U_;
+    volScalarField& nut = this->nut_;
+    fv::options& fvOptions(fv::options::New(this->mesh_));
+
+    BasicTurbulenceModel::correct();
+
+    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
+
+    tmp<volTensorField> tgradU = fvc::grad(U);
+    volScalarField S2(2*magSqr(symm(tgradU())));
+    volScalarField GbyNu((tgradU() && dev(twoSymm(tgradU()))));
+    volScalarField G(this->GName(), nut*GbyNu);
+    tgradU.clear();
+
+    // Update omega and G at the wall
+    omega_.boundaryFieldRef().updateCoeffs();
+
+    volScalarField CDkOmega
+    (
+        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
+    );
+
+    volScalarField F1(this->F1(CDkOmega));
+    volScalarField F23(this->F23());
+
+    {
+        volScalarField gamma(this->gamma(F1));
+        volScalarField beta(this->beta(F1));
+
+        // Turbulent frequency equation
+        tmp<fvScalarMatrix> omegaEqn
+        (
+            fvm::ddt(alpha, rho, omega_)
+          + fvm::div(alphaRhoPhi, omega_)
+          - fvm::laplacian(alpha*rho*DomegaEff(F1), omega_)
+         ==
+            alpha*rho*gamma
+           *min
+            (
+                GbyNu,
+                (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23*sqrt(S2))
+            )
+          - fvm::SuSp((2.0/3.0)*alpha*rho*gamma*divU, omega_)
+          - fvm::Sp(alpha*rho*beta*omega_, omega_)
+          - fvm::SuSp
+            (
+                alpha*rho*(F1 - scalar(1))*CDkOmega/omega_,
+                omega_
+            )
+          + Qsas(S2, gamma, beta)
+          + omegaSource()
+          + fvOptions(alpha, rho, omega_)
+        );
+
+        omegaEqn.ref().relax();
+        fvOptions.constrain(omegaEqn.ref());
+        omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
+        solve(omegaEqn);
+        fvOptions.correct(omega_);
+        bound(omega_, this->omegaMin_);
+    }
+
+    // Turbulent kinetic energy equation
+    tmp<fvScalarMatrix> kEqn
+    (
+        fvm::ddt(alpha, rho, k_)
+      + fvm::div(alphaRhoPhi, k_)
+      - fvm::laplacian(alpha*rho*DkEff(F1), k_)
+     ==
+        min(alpha*rho*G, (c1_*betaStar_)*alpha*rho*k_*omega_)
+      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
+      - fvm::Sp(alpha*rho*epsilonByk(F1, F23), k_)
+      + kSource()
+      + fvOptions(alpha, rho, k_)
+    );
+
+    kEqn.ref().relax();
+    fvOptions.constrain(kEqn.ref());
+    solve(kEqn);
+    fvOptions.correct(k_);
+    bound(k_, this->kMin_);
+
+    correctNut(S2, F23);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.H b/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.H
new file mode 100644
index 0000000000000000000000000000000000000000..5adf9961ab1c793287603c211fe563e73302847e
--- /dev/null
+++ b/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.H
@@ -0,0 +1,331 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::kOmegaSST
+
+Group
+    grpTurbulence
+
+Description
+    Implementation of the k-omega-SST turbulence model for
+    incompressible and compressible flows.
+
+    Turbulence model described in:
+    \verbatim
+        Menter, F. R. & Esch, T. (2001).
+        Elements of Industrial Heat Transfer Prediction.
+        16th Brazilian Congress of Mechanical Engineering (COBEM).
+    \endverbatim
+
+    with updated coefficients from
+    \verbatim
+        Menter, F. R., Kuntz, M., and Langtry, R. (2003).
+        Ten Years of Industrial Experience with the SST Turbulence Model.
+        Turbulence, Heat and Mass Transfer 4, ed: K. Hanjalic, Y. Nagano,
+        & M. Tummers, Begell House, Inc., 625 - 632.
+    \endverbatim
+
+    but with the consistent production terms from the 2001 paper as form in the
+    2003 paper is a typo, see
+    \verbatim
+        http://turbmodels.larc.nasa.gov/sst.html
+    \endverbatim
+
+    and the addition of the optional F3 term for rough walls from
+    \verbatim
+        Hellsten, A. (1998).
+        "Some Improvements in Menter’s k-omega-SST turbulence model"
+        29th AIAA Fluid Dynamics Conference, AIAA-98-2554.
+    \endverbatim
+
+    Note that this implementation is written in terms of alpha diffusion
+    coefficients rather than the more traditional sigma (alpha = 1/sigma) so
+    that the blending can be applied to all coefficuients in a consistent
+    manner.  The paper suggests that sigma is blended but this would not be
+    consistent with the blending of the k-epsilon and k-omega models.
+
+    Also note that the error in the last term of equation (2) relating to
+    sigma has been corrected.
+
+    Wall-functions are applied in this implementation by using equations (14)
+    to specify the near-wall omega as appropriate.
+
+    The blending functions (15) and (16) are not currently used because of the
+    uncertainty in their origin, range of applicability and that if y+ becomes
+    sufficiently small blending u_tau in this manner clearly becomes nonsense.
+
+    The default model coefficients are
+    \verbatim
+        kOmegaSSTCoeffs
+        {
+            alphaK1     0.85;
+            alphaK2     1.0;
+            alphaOmega1 0.5;
+            alphaOmega2 0.856;
+            beta1       0.075;
+            beta2       0.0828;
+            betaStar    0.09;
+            gamma1      5/9;
+            gamma2      0.44;
+            a1          0.31;
+            b1          1.0;
+            c1          10.0;
+            F3          no;
+        }
+    \endverbatim
+
+SourceFiles
+    kOmegaSST.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kOmegaSSTBase_H
+#define kOmegaSSTBase_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class kOmegaSST Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+class kOmegaSST
+:
+    public TurbulenceModel
+{
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        kOmegaSST(const kOmegaSST&);
+        void operator=(const kOmegaSST&);
+
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar alphaK1_;
+            dimensionedScalar alphaK2_;
+
+            dimensionedScalar alphaOmega1_;
+            dimensionedScalar alphaOmega2_;
+
+            dimensionedScalar gamma1_;
+            dimensionedScalar gamma2_;
+
+            dimensionedScalar beta1_;
+            dimensionedScalar beta2_;
+
+            dimensionedScalar betaStar_;
+
+            dimensionedScalar a1_;
+            dimensionedScalar b1_;
+            dimensionedScalar c1_;
+
+            Switch F3_;
+
+
+        // Fields
+
+            //- Wall distance
+            //  Note: different to wall distance in parent RASModel
+            //  which is for near-wall cells only
+            const volScalarField& y_;
+
+            volScalarField k_;
+            volScalarField omega_;
+
+
+    // Private Member Functions
+
+        tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
+        tmp<volScalarField> F2() const;
+        tmp<volScalarField> F3() const;
+        tmp<volScalarField> F23() const;
+
+        tmp<volScalarField> blend
+        (
+            const volScalarField& F1,
+            const dimensionedScalar& psi1,
+            const dimensionedScalar& psi2
+        ) const
+        {
+            return F1*(psi1 - psi2) + psi2;
+        }
+
+        tmp<volScalarField> alphaK(const volScalarField& F1) const
+        {
+            return blend(F1, alphaK1_, alphaK2_);
+        }
+
+        tmp<volScalarField> alphaOmega(const volScalarField& F1) const
+        {
+            return blend(F1, alphaOmega1_, alphaOmega2_);
+        }
+
+        tmp<volScalarField> beta(const volScalarField& F1) const
+        {
+            return blend(F1, beta1_, beta2_);
+        }
+
+        tmp<volScalarField> gamma(const volScalarField& F1) const
+        {
+            return blend(F1, gamma1_, gamma2_);
+        }
+
+        void correctNut(const volScalarField& S2, const volScalarField& F2);
+
+
+    // Protected Member Functions
+
+        virtual void correctNut();
+
+        //- Return epsilon/k which for standard RAS is betaStar*omega
+        virtual tmp<volScalarField> epsilonByk
+        (
+            const volScalarField& F1,
+            const volScalarField& F2
+        ) const;
+
+        virtual tmp<fvScalarMatrix> kSource() const;
+
+        virtual tmp<fvScalarMatrix> omegaSource() const;
+
+        virtual tmp<fvScalarMatrix> Qsas
+        (
+            const volScalarField& S2,
+            const volScalarField& gamma,
+            const volScalarField& beta
+        ) const;
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    // Constructors
+
+        //- Construct from components
+        kOmegaSST
+        (
+            const word& type,
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName
+        );
+
+
+    //- Destructor
+    virtual ~kOmegaSST()
+    {}
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return the effective diffusivity for k
+        tmp<volScalarField> DkEff(const volScalarField& F1) const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu())
+            );
+        }
+
+        //- Return the effective diffusivity for omega
+        tmp<volScalarField> DomegaEff(const volScalarField& F1) const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    "DomegaEff",
+                    alphaOmega(F1)*this->nut_ + this->nu()
+                )
+            );
+        }
+
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k_;
+        }
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const
+        {
+              return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        "epsilon",
+                        this->mesh_.time().timeName(),
+                        this->mesh_
+                    ),
+                    betaStar_*k_*omega_,
+                    omega_.boundaryField().types()
+                )
+            );
+        }
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> omega() const
+        {
+            return omega_;
+        }
+
+        //- Solve the turbulence equations and correct the turbulence viscosity
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#ifdef NoRepository
+    #include "kOmegaSSTBase.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#endif
+
+// ************************************************************************* //
diff --git a/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.C b/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.C
index 319febe26f2ee3b7240d7728997ab679d88959cf..c5cc7ebdb70fc625cd9e5604035ae4dd83030b6a 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.C
+++ b/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.C
@@ -80,6 +80,28 @@ Foam::LESModel<BasicTurbulenceModel>::LESModel
         )
     ),
 
+    epsilonMin_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "epsilonMin",
+            LESDict_,
+            kMin_.dimensions()/dimTime,
+            SMALL
+        )
+    ),
+
+    omegaMin_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "omegaMin",
+            LESDict_,
+            dimless/dimTime,
+            SMALL
+        )
+    ),
+
     delta_
     (
         LESdelta::New
diff --git a/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.H b/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.H
index 1449f4b8f64cc0b978652e8a4d9e02ae2fef4010..297e86b3d3002ec4e06d4a47a5f0583b93480d65 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.H
+++ b/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.H
@@ -78,6 +78,12 @@ protected:
         //- Lower limit of k
         dimensionedScalar kMin_;
 
+        //- Lower limit of epsilon
+        dimensionedScalar epsilonMin_;
+
+        //- Lower limit for omega
+        dimensionedScalar omegaMin_;
+
         //- Run-time selectable delta model
         autoPtr<Foam::LESdelta> delta_;
 
diff --git a/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.C b/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.C
new file mode 100644
index 0000000000000000000000000000000000000000..0ad21dd8cb49221e7f47d4b9602c51876bbe234b
--- /dev/null
+++ b/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.C
@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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 "kOmegaSSTDES.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::Lt() const
+{
+    return sqrt(this->k_)/(this->betaStar_*this->omega_);
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::FDES
+(
+    const volScalarField& F1,
+    const volScalarField& F2
+) const
+{
+    switch (FSST_)
+    {
+        case 0:
+            return max(Lt()/(CDES_*this->delta()), scalar(1));
+        case 1:
+            return max(Lt()*(1 - F1)/(CDES_*this->delta()), scalar(1));
+        case 2:
+            return max(Lt()*(1 - F2)/(CDES_*this->delta()), scalar(1));
+        default:
+            FatalErrorInFunction
+                << "Incorrect FSST = " << FSST_ << ", should be 0, 1 or 2"
+                << exit(FatalError);
+            return F1;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::epsilonByk
+(
+    const volScalarField& F1,
+    const volScalarField& F2
+) const
+{
+    return this->betaStar_*this->omega_*FDES(F1, F2);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+kOmegaSSTDES<BasicTurbulenceModel>::kOmegaSSTDES
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    kOmegaSST
+    <
+        LESeddyViscosity<BasicTurbulenceModel>,
+        BasicTurbulenceModel
+    >
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    CDES_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "CDES",
+            this->coeffDict_,
+            0.61
+        )
+    ),
+    FSST_(this->coeffDict_.lookupOrDefault("FSST", 2))
+{
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool kOmegaSSTDES<BasicTurbulenceModel>::read()
+{
+    if
+    (
+        kOmegaSST<LESeddyViscosity<BasicTurbulenceModel>, BasicTurbulenceModel>
+        ::read()
+    )
+    {
+        CDES_.readIfPresent(this->coeffDict());
+        this->coeffDict().readIfPresent("FSST", FSST_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.H b/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.H
new file mode 100644
index 0000000000000000000000000000000000000000..014f30e9bfb9539d79ee1e485aea2214ed936a74
--- /dev/null
+++ b/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.H
@@ -0,0 +1,174 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+Class
+    Foam::LESModels::kOmegaSST
+
+Description
+    Implementation of the k-omega-SST-DES turbulence model for
+    incompressible and compressible flows.
+
+    DES model described in:
+    \verbatim
+        Menter, F. R., Kuntz, M., and Langtry, R. (2003).
+        Ten Years of Industrial Experience with the SST Turbulence Model.
+        Turbulence, Heat and Mass Transfer 4, ed: K. Hanjalic, Y. Nagano,
+        & M. Tummers, Begell House, Inc., 625 - 632.
+    \endverbatim
+
+    Optional support for zonal filtering based on F1 or F2 is provided as
+    described in the paper.
+
+    For further details of the implementation of the base k-omega-SST model
+    see Foam::kOmegaSST.
+
+Group
+    grpLESTurbulence
+
+SeeAlso
+    Foam::kOmegaSST
+
+SourceFiles
+    kOmegaSST.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kOmegaSSTDES_H
+#define kOmegaSSTDES_H
+
+#include "kOmegaSSTBase.H"
+#include "LESModel.H"
+#include "LESeddyViscosity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class kOmegaSST Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class kOmegaSSTDES
+:
+    public Foam::kOmegaSST
+    <
+        LESeddyViscosity<BasicTurbulenceModel>,
+        BasicTurbulenceModel
+    >
+{
+
+protected:
+
+    // Protected data
+
+        // Model constants
+
+            //- DES coefficient
+            dimensionedScalar CDES_;
+
+            //- Zonal filter choice
+            //
+            //  - 0: no filtering
+            //  - 1: (1 - F1)
+            //  - 2: (1 - F2)
+            direction FSST_;
+
+
+    // Protected Member Functions
+
+        //- Return the turbulent length-scale
+        tmp<volScalarField> Lt() const;
+
+        //- The DES dissipation-rate multiplier with options zonal filtering
+        //  based on either F1 or F2
+        virtual tmp<volScalarField> FDES
+        (
+            const volScalarField& F1,
+            const volScalarField& F2
+        ) const;
+
+        //- Return epsilon/k which for standard RAS is betaStar*omega
+        virtual tmp<volScalarField> epsilonByk
+        (
+            const volScalarField& F1,
+            const volScalarField& F2
+        ) const;
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("kOmegaSSTDES");
+
+
+    // Constructors
+
+        //- Construct from components
+        kOmegaSSTDES
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~kOmegaSSTDES()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#ifdef NoRepository
+    #include "kOmegaSSTDES.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#endif
+
+// ************************************************************************* //
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C
index 85f73846551d237bb5ec5a7c3dcf2c1b74174aba..f26547b509de8adc13a9e1422ffe36d08f79aaf1 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,9 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "kOmegaSST.H"
-#include "fvOptions.H"
-#include "bound.H"
-#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -35,146 +32,6 @@ namespace Foam
 namespace RASModels
 {
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F1
-(
-    const volScalarField& CDkOmega
-) const
-{
-    tmp<volScalarField> CDkOmegaPlus = max
-    (
-        CDkOmega,
-        dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
-    );
-
-    tmp<volScalarField> arg1 = min
-    (
-        min
-        (
-            max
-            (
-                (scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
-                scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
-            ),
-            (4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
-        ),
-        scalar(10)
-    );
-
-    return tanh(pow4(arg1));
-}
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F2() const
-{
-    tmp<volScalarField> arg2 = min
-    (
-        max
-        (
-            (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
-            scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_)
-        ),
-        scalar(100)
-    );
-
-    return tanh(sqr(arg2));
-}
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F3() const
-{
-    tmp<volScalarField> arg3 = min
-    (
-        150*(this->mu()/this->rho_)/(omega_*sqr(y_)),
-        scalar(10)
-    );
-
-    return 1 - tanh(pow4(arg3));
-}
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F23() const
-{
-    tmp<volScalarField> f23(F2());
-
-    if (F3_)
-    {
-        f23.ref() *= F3();
-    }
-
-    return f23;
-}
-
-
-template<class BasicTurbulenceModel>
-void kOmegaSST<BasicTurbulenceModel>::correctNut(const volScalarField& S2)
-{
-    this->nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
-    this->nut_.correctBoundaryConditions();
-    fv::options::New(this->mesh_).correct(this->nut_);
-
-    BasicTurbulenceModel::correctNut();
-}
-
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class BasicTurbulenceModel>
-void kOmegaSST<BasicTurbulenceModel>::correctNut()
-{
-    correctNut(2*magSqr(symm(fvc::grad(this->U_))));
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<fvScalarMatrix> kOmegaSST<BasicTurbulenceModel>::kSource() const
-{
-    return tmp<fvScalarMatrix>
-    (
-        new fvScalarMatrix
-        (
-            k_,
-            dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
-        )
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<fvScalarMatrix> kOmegaSST<BasicTurbulenceModel>::omegaSource() const
-{
-    return tmp<fvScalarMatrix>
-    (
-        new fvScalarMatrix
-        (
-            omega_,
-            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
-        )
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<fvScalarMatrix> kOmegaSST<BasicTurbulenceModel>::Qsas
-(
-    const volScalarField& S2,
-    const volScalarField& gamma,
-    const volScalarField& beta
-) const
-{
-    return tmp<fvScalarMatrix>
-    (
-        new fvScalarMatrix
-        (
-            omega_,
-            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
-        )
-    );
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
@@ -190,7 +47,11 @@ kOmegaSST<BasicTurbulenceModel>::kOmegaSST
     const word& type
 )
 :
-    eddyViscosity<RASModel<BasicTurbulenceModel>>
+    Foam::kOmegaSST
+    <
+        eddyViscosity<RASModel<BasicTurbulenceModel>>,
+        BasicTurbulenceModel
+    >
     (
         type,
         alpha,
@@ -200,156 +61,8 @@ kOmegaSST<BasicTurbulenceModel>::kOmegaSST
         phi,
         transport,
         propertiesName
-    ),
-
-    alphaK1_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "alphaK1",
-            this->coeffDict_,
-            0.85
-        )
-    ),
-    alphaK2_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "alphaK2",
-            this->coeffDict_,
-            1.0
-        )
-    ),
-    alphaOmega1_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "alphaOmega1",
-            this->coeffDict_,
-            0.5
-        )
-    ),
-    alphaOmega2_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "alphaOmega2",
-            this->coeffDict_,
-            0.856
-        )
-    ),
-    gamma1_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "gamma1",
-            this->coeffDict_,
-            5.0/9.0
-        )
-    ),
-    gamma2_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "gamma2",
-            this->coeffDict_,
-            0.44
-        )
-    ),
-    beta1_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "beta1",
-            this->coeffDict_,
-            0.075
-        )
-    ),
-    beta2_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "beta2",
-            this->coeffDict_,
-            0.0828
-        )
-    ),
-    betaStar_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "betaStar",
-            this->coeffDict_,
-            0.09
-        )
-    ),
-    a1_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "a1",
-            this->coeffDict_,
-            0.31
-        )
-    ),
-    b1_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "b1",
-            this->coeffDict_,
-            1.0
-        )
-    ),
-    c1_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "c1",
-            this->coeffDict_,
-            10.0
-        )
-    ),
-    F3_
-    (
-        Switch::lookupOrAddToDict
-        (
-            "F3",
-            this->coeffDict_,
-            false
-        )
-    ),
-
-    y_(wallDist::New(this->mesh_).y()),
-
-    k_
-    (
-        IOobject
-        (
-            IOobject::groupName("k", U.group()),
-            this->runTime_.timeName(),
-            this->mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        this->mesh_
-    ),
-    omega_
-    (
-        IOobject
-        (
-            IOobject::groupName("omega", U.group()),
-            this->runTime_.timeName(),
-            this->mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        this->mesh_
     )
 {
-    bound(k_, this->kMin_);
-    bound(omega_, this->omegaMin_);
-
     if (type == typeName)
     {
         this->printCoeffs(type);
@@ -357,133 +70,6 @@ kOmegaSST<BasicTurbulenceModel>::kOmegaSST
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class BasicTurbulenceModel>
-bool kOmegaSST<BasicTurbulenceModel>::read()
-{
-    if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
-    {
-        alphaK1_.readIfPresent(this->coeffDict());
-        alphaK2_.readIfPresent(this->coeffDict());
-        alphaOmega1_.readIfPresent(this->coeffDict());
-        alphaOmega2_.readIfPresent(this->coeffDict());
-        gamma1_.readIfPresent(this->coeffDict());
-        gamma2_.readIfPresent(this->coeffDict());
-        beta1_.readIfPresent(this->coeffDict());
-        beta2_.readIfPresent(this->coeffDict());
-        betaStar_.readIfPresent(this->coeffDict());
-        a1_.readIfPresent(this->coeffDict());
-        b1_.readIfPresent(this->coeffDict());
-        c1_.readIfPresent(this->coeffDict());
-        F3_.readIfPresent("F3", this->coeffDict());
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-template<class BasicTurbulenceModel>
-void kOmegaSST<BasicTurbulenceModel>::correct()
-{
-    if (!this->turbulence_)
-    {
-        return;
-    }
-
-    // Local references
-    const alphaField& alpha = this->alpha_;
-    const rhoField& rho = this->rho_;
-    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
-    const volVectorField& U = this->U_;
-    volScalarField& nut = this->nut_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
-
-    eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
-
-    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
-
-    tmp<volTensorField> tgradU = fvc::grad(U);
-    volScalarField S2(2*magSqr(symm(tgradU())));
-    volScalarField GbyNu((tgradU() && dev(twoSymm(tgradU()))));
-    volScalarField G(this->GName(), nut*GbyNu);
-    tgradU.clear();
-
-    // Update omega and G at the wall
-    omega_.boundaryFieldRef().updateCoeffs();
-
-    volScalarField CDkOmega
-    (
-        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
-    );
-
-    volScalarField F1(this->F1(CDkOmega));
-
-    {
-        volScalarField gamma(this->gamma(F1));
-        volScalarField beta(this->beta(F1));
-
-        // Turbulent frequency equation
-        tmp<fvScalarMatrix> omegaEqn
-        (
-            fvm::ddt(alpha, rho, omega_)
-          + fvm::div(alphaRhoPhi, omega_)
-          - fvm::laplacian(alpha*rho*DomegaEff(F1), omega_)
-         ==
-            alpha*rho*gamma
-           *min
-            (
-                GbyNu,
-                (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2))
-            )
-          - fvm::SuSp((2.0/3.0)*alpha*rho*gamma*divU, omega_)
-          - fvm::Sp(alpha*rho*beta*omega_, omega_)
-          - fvm::SuSp
-            (
-                alpha*rho*(F1 - scalar(1))*CDkOmega/omega_,
-                omega_
-            )
-          + Qsas(S2, gamma, beta)
-          + omegaSource()
-          + fvOptions(alpha, rho, omega_)
-        );
-
-        omegaEqn.ref().relax();
-        fvOptions.constrain(omegaEqn.ref());
-        omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
-        solve(omegaEqn);
-        fvOptions.correct(omega_);
-        bound(omega_, this->omegaMin_);
-    }
-
-    // Turbulent kinetic energy equation
-    tmp<fvScalarMatrix> kEqn
-    (
-        fvm::ddt(alpha, rho, k_)
-      + fvm::div(alphaRhoPhi, k_)
-      - fvm::laplacian(alpha*rho*DkEff(F1), k_)
-     ==
-        min(alpha*rho*G, (c1_*betaStar_)*alpha*rho*k_*omega_)
-      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
-      - fvm::Sp(alpha*rho*betaStar_*omega_, k_)
-      + kSource()
-      + fvOptions(alpha, rho, k_)
-    );
-
-    kEqn.ref().relax();
-    fvOptions.constrain(kEqn.ref());
-    solve(kEqn);
-    fvOptions.correct(k_);
-    bound(k_, this->kMin_);
-
-    correctNut(S2);
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace RASModels
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H
index 79f29427199e3ac0e1429c734a83f0dd7a2ac25c..351b0cee5958419b06fb29f035833e88f0c4e5ea 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,73 +27,8 @@ Class
 Group
     grpRASTurbulence
 
-Description
-    Implementation of the k-omega-SST turbulence model for
-    incompressible and compressible flows.
-
-    Turbulence model described in:
-    \verbatim
-        Menter, F. R. & Esch, T. (2001).
-        Elements of Industrial Heat Transfer Prediction.
-        16th Brazilian Congress of Mechanical Engineering (COBEM).
-    \endverbatim
-
-    with updated coefficients from
-    \verbatim
-        Menter, F. R., Kuntz, M., and Langtry, R. (2003).
-        Ten Years of Industrial Experience with the SST Turbulence Model.
-        Turbulence, Heat and Mass Transfer 4, ed: K. Hanjalic, Y. Nagano,
-        & M. Tummers, Begell House, Inc., 625 - 632.
-    \endverbatim
-
-    but with the consistent production terms from the 2001 paper as form in the
-    2003 paper is a typo, see
-    \verbatim
-        http://turbmodels.larc.nasa.gov/sst.html
-    \endverbatim
-
-    and the addition of the optional F3 term for rough walls from
-    \verbatim
-        Hellsten, A. (1998).
-        "Some Improvements in Menter’s k-omega-SST turbulence model"
-        29th AIAA Fluid Dynamics Conference, AIAA-98-2554.
-    \endverbatim
-
-    Note that this implementation is written in terms of alpha diffusion
-    coefficients rather than the more traditional sigma (alpha = 1/sigma) so
-    that the blending can be applied to all coefficuients in a consistent
-    manner.  The paper suggests that sigma is blended but this would not be
-    consistent with the blending of the k-epsilon and k-omega models.
-
-    Also note that the error in the last term of equation (2) relating to
-    sigma has been corrected.
-
-    Wall-functions are applied in this implementation by using equations (14)
-    to specify the near-wall omega as appropriate.
-
-    The blending functions (15) and (16) are not currently used because of the
-    uncertainty in their origin, range of applicability and that if y+ becomes
-    sufficiently small blending u_tau in this manner clearly becomes nonsense.
-
-    The default model coefficients are
-    \verbatim
-        kOmegaSSTCoeffs
-        {
-            alphaK1     0.85;
-            alphaK2     1.0;
-            alphaOmega1 0.5;
-            alphaOmega2 0.856;
-            beta1       0.075;
-            beta2       0.0828;
-            betaStar    0.09;
-            gamma1      5/9;
-            gamma2      0.44;
-            a1          0.31;
-            b1          1.0;
-            c1          10.0;
-            F3          no;
-        }
-    \endverbatim
+SeeAlso
+    Foam::kOmegaSST
 
 SourceFiles
     kOmegaSST.C
@@ -103,6 +38,7 @@ SourceFiles
 #ifndef kOmegaSST_H
 #define kOmegaSST_H
 
+#include "kOmegaSSTBase.H"
 #include "RASModel.H"
 #include "eddyViscosity.H"
 
@@ -114,111 +50,18 @@ namespace RASModels
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class kOmegaSST Declaration
+                          Class kOmegaSST Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class BasicTurbulenceModel>
 class kOmegaSST
 :
-    public eddyViscosity<RASModel<BasicTurbulenceModel>>
+    public Foam::kOmegaSST
+    <
+        eddyViscosity<RASModel<BasicTurbulenceModel>>,
+        BasicTurbulenceModel
+    >
 {
-    // Private Member Functions
-
-        // Disallow default bitwise copy construct and assignment
-        kOmegaSST(const kOmegaSST&);
-        void operator=(const kOmegaSST&);
-
-
-protected:
-
-    // Protected data
-
-        // Model coefficients
-
-            dimensionedScalar alphaK1_;
-            dimensionedScalar alphaK2_;
-
-            dimensionedScalar alphaOmega1_;
-            dimensionedScalar alphaOmega2_;
-
-            dimensionedScalar gamma1_;
-            dimensionedScalar gamma2_;
-
-            dimensionedScalar beta1_;
-            dimensionedScalar beta2_;
-
-            dimensionedScalar betaStar_;
-
-            dimensionedScalar a1_;
-            dimensionedScalar b1_;
-            dimensionedScalar c1_;
-
-            Switch F3_;
-
-
-        // Fields
-
-            //- Wall distance
-            //  Note: different to wall distance in parent RASModel
-            //  which is for near-wall cells only
-            const volScalarField& y_;
-
-            volScalarField k_;
-            volScalarField omega_;
-
-
-    // Private Member Functions
-
-        tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
-        tmp<volScalarField> F2() const;
-        tmp<volScalarField> F3() const;
-        tmp<volScalarField> F23() const;
-
-        tmp<volScalarField> blend
-        (
-            const volScalarField& F1,
-            const dimensionedScalar& psi1,
-            const dimensionedScalar& psi2
-        ) const
-        {
-            return F1*(psi1 - psi2) + psi2;
-        }
-
-        tmp<volScalarField> alphaK(const volScalarField& F1) const
-        {
-            return blend(F1, alphaK1_, alphaK2_);
-        }
-
-        tmp<volScalarField> alphaOmega(const volScalarField& F1) const
-        {
-            return blend(F1, alphaOmega1_, alphaOmega2_);
-        }
-
-        tmp<volScalarField> beta(const volScalarField& F1) const
-        {
-            return blend(F1, beta1_, beta2_);
-        }
-
-        tmp<volScalarField> gamma(const volScalarField& F1) const
-        {
-            return blend(F1, gamma1_, gamma2_);
-        }
-
-        void correctNut(const volScalarField& S2);
-
-
-    // Protected Member Functions
-
-        virtual void correctNut();
-        virtual tmp<fvScalarMatrix> kSource() const;
-        virtual tmp<fvScalarMatrix> omegaSource() const;
-        virtual tmp<fvScalarMatrix> Qsas
-        (
-            const volScalarField& S2,
-            const volScalarField& gamma,
-            const volScalarField& beta
-        ) const;
-
 
 public:
 
@@ -250,68 +93,6 @@ public:
     //- Destructor
     virtual ~kOmegaSST()
     {}
-
-
-    // Member Functions
-
-        //- Re-read model coefficients if they have changed
-        virtual bool read();
-
-        //- Return the effective diffusivity for k
-        tmp<volScalarField> DkEff(const volScalarField& F1) const
-        {
-            return tmp<volScalarField>
-            (
-                new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu())
-            );
-        }
-
-        //- Return the effective diffusivity for omega
-        tmp<volScalarField> DomegaEff(const volScalarField& F1) const
-        {
-            return tmp<volScalarField>
-            (
-                new volScalarField
-                (
-                    "DomegaEff",
-                    alphaOmega(F1)*this->nut_ + this->nu()
-                )
-            );
-        }
-
-        //- Return the turbulence kinetic energy
-        virtual tmp<volScalarField> k() const
-        {
-            return k_;
-        }
-
-        //- Return the turbulence kinetic energy dissipation rate
-        virtual tmp<volScalarField> epsilon() const
-        {
-              return tmp<volScalarField>
-            (
-                new volScalarField
-                (
-                    IOobject
-                    (
-                        "epsilon",
-                        this->mesh_.time().timeName(),
-                        this->mesh_
-                    ),
-                    betaStar_*k_*omega_,
-                    omega_.boundaryField().types()
-                )
-            );
-        }
-
-        //- Return the turbulence kinetic energy dissipation rate
-        virtual tmp<volScalarField> omega() const
-        {
-            return omega_;
-        }
-
-        //- Solve the turbulence equations and correct the turbulence viscosity
-        virtual void correct();
 };