diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C
index 76fd9cf0c3583b4613cc06fa5a5db647280ea226..4353b5a7481e9a12bbcadd36c4912fc5634d8121 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C
@@ -24,8 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "kOmegaSST.H"
-#include "bound.H"
-#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -34,91 +32,22 @@ 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));
-}
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F3() const
+void kOmegaSST<BasicTurbulenceModel>::correctNut(const volScalarField& S2)
 {
-    tmp<volScalarField> arg3 = min
+    // Correct the turbulence viscosity
+    kOmegaSSTBase<eddyViscosity<RASModel<BasicTurbulenceModel> > >::correctNut
     (
-        150*(this->mu()/this->rho_)/(omega_*sqr(y_)),
-        scalar(10)
+        S2
     );
 
-    return 1 - tanh(pow4(arg3));
-}
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSST<BasicTurbulenceModel>::kOmegaSST::F23() const
-{
-    tmp<volScalarField> f23(F2());
-
-    if (F3_)
-    {
-        f23() *= 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();
-
+    // Correct the turbulence thermal diffusivity
     BasicTurbulenceModel::correctNut();
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
 template<class BasicTurbulenceModel>
 void kOmegaSST<BasicTurbulenceModel>::correctNut()
 {
@@ -126,53 +55,6 @@ void kOmegaSST<BasicTurbulenceModel>::correctNut()
 }
 
 
-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>
@@ -188,7 +70,7 @@ kOmegaSST<BasicTurbulenceModel>::kOmegaSST
     const word& type
 )
 :
-    eddyViscosity<RASModel<BasicTurbulenceModel> >
+    kOmegaSSTBase<eddyViscosity<RASModel<BasicTurbulenceModel> > >
     (
         type,
         alpha,
@@ -198,286 +80,17 @@ 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_);
+    correctNut();
 
     if (type == typeName)
     {
-        correctNut();
         this->printCoeffs(type);
     }
 }
 
 
-// * * * * * * * * * * * * * * * 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_;
-
-    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_.boundaryField().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()
-        );
-
-        omegaEqn().relax();
-
-        omegaEqn().boundaryManipulate(omega_.boundaryField());
-
-        solve(omegaEqn);
-        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()
-    );
-
-    kEqn().relax();
-    solve(kEqn);
-    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 1bbbedc91df55a2e51dc9f80214520f217564d2d..17327ad9ddeabe4c2cf491f83ea3da1fd39b91f5 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H
@@ -31,6 +31,8 @@ Description
     Implementation of the k-omega-SST turbulence model for
     incompressible and compressible flows.
 
+    Light wrapper around base class.
+
     Turbulence model described in:
     \verbatim
         Menter, F. R. & Esch, T. (2001).
@@ -98,11 +100,15 @@ Description
 SourceFiles
     kOmegaSST.C
 
+SeeAlso
+    kOmegaSSTBase.H
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef kOmegaSST_H
 #define kOmegaSST_H
 
+#include "kOmegaSSTBase.H"
 #include "RASModel.H"
 #include "eddyViscosity.H"
 
@@ -114,110 +120,21 @@ namespace RASModels
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class kOmegaSST Declaration
+                          Class kOmegaSST Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class BasicTurbulenceModel>
 class kOmegaSST
 :
-    public eddyViscosity<RASModel<BasicTurbulenceModel> >
+    public kOmegaSSTBase<eddyViscosity<RASModel<BasicTurbulenceModel> > >
 {
-    // Private Member Functions
-
-        // Disallow default bitwise copy construct and assignment
-        kOmegaSST(const kOmegaSST&);
-        kOmegaSST& 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
+    // Protecetd Member Functions
 
+        virtual void correctNut(const volScalarField& S2);
         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 +167,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();
 };
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.C b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.C
new file mode 100644
index 0000000000000000000000000000000000000000..19943c0c1bb48e2c3fefd878df4aadcb405504c1
--- /dev/null
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.C
@@ -0,0 +1,480 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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 "bound.H"
+#include "wallDist.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::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 BasicEddyViscosityModel>
+tmp<volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::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 BasicEddyViscosityModel>
+tmp<volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::F3() const
+{
+    tmp<volScalarField> arg3 = min
+    (
+        150*(this->mu()/this->rho_)/(omega_*sqr(y_)),
+        scalar(10)
+    );
+
+    return 1 - tanh(pow4(arg3));
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::F23() const
+{
+    tmp<volScalarField> f23(F2());
+
+    if (F3_)
+    {
+        f23() *= F3();
+    }
+
+    return f23;
+}
+
+
+template<class BasicEddyViscosityModel>
+void kOmegaSSTBase<BasicEddyViscosityModel>::correctNut
+(
+    const volScalarField& S2
+)
+{
+    // Correct the turbulence viscosity
+    this->nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
+    this->nut_.correctBoundaryConditions();
+}
+
+
+template<class BasicEddyViscosityModel>
+void kOmegaSSTBase<BasicEddyViscosityModel>::correctNut()
+{
+    correctNut(2*magSqr(symm(fvc::grad(this->U_))));
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::kSource() const
+{
+    return tmp<fvScalarMatrix>
+    (
+        new fvScalarMatrix
+        (
+            k_,
+            dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
+        )
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::omegaSource() const
+{
+    return tmp<fvScalarMatrix>
+    (
+        new fvScalarMatrix
+        (
+            omega_,
+            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
+        )
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::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 BasicEddyViscosityModel>
+kOmegaSSTBase<BasicEddyViscosityModel>::kOmegaSSTBase
+(
+    const word& type,
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+:
+    BasicEddyViscosityModel
+    (
+        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 BasicEddyViscosityModel>
+bool kOmegaSSTBase<BasicEddyViscosityModel>::read()
+{
+    if (BasicEddyViscosityModel::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 BasicEddyViscosityModel>
+void kOmegaSSTBase<BasicEddyViscosityModel>::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_;
+
+    BasicEddyViscosityModel::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_.boundaryField().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()
+        );
+
+        omegaEqn().relax();
+
+        omegaEqn().boundaryManipulate(omega_.boundaryField());
+
+        solve(omegaEqn);
+        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()
+    );
+
+    kEqn().relax();
+    solve(kEqn);
+    bound(k_, this->kMin_);
+
+    correctNut(S2);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.H b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.H
new file mode 100644
index 0000000000000000000000000000000000000000..eaa4bd76526274ba15f7808782e79748e5523270
--- /dev/null
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.H
@@ -0,0 +1,320 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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::kOmegaSSTBase
+
+Description
+    Base class 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
+        kOmegaSSTBaseCoeffs
+        {
+            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
+    kOmegaSSTBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kOmegaSSTBase_H
+#define kOmegaSSTBase_H
+
+#include "RASModel.H"
+#include "eddyViscosity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class kOmegaSSTBase Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicEddyViscosityModel>
+class kOmegaSSTBase
+:
+    public BasicEddyViscosityModel
+{
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        kOmegaSSTBase(const kOmegaSSTBase&);
+        kOmegaSSTBase& operator=(const kOmegaSSTBase&);
+
+
+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_;
+
+
+    // Protected 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);
+
+        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:
+
+    typedef typename BasicEddyViscosityModel::alphaField alphaField;
+    typedef typename BasicEddyViscosityModel::rhoField rhoField;
+    typedef typename BasicEddyViscosityModel::transportModel transportModel;
+
+
+    // Constructors
+
+        //- Construct from components
+        kOmegaSSTBase
+        (
+            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 ~kOmegaSSTBase()
+    {}
+
+
+    // 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/turbulenceModelDoc.H b/src/TurbulenceModels/turbulenceModels/doc/turbulenceModelDoc.H
similarity index 100%
rename from src/TurbulenceModels/turbulenceModels/turbulenceModelDoc.H
rename to src/TurbulenceModels/turbulenceModels/doc/turbulenceModelDoc.H