diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
index d26226feda..9b4c373420 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
@@ -129,7 +129,12 @@ Foam::kineticTheoryModel::kineticTheoryModel
         U.mesh(),
         dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0)
     )
-{}
+{
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -140,6 +145,199 @@ Foam::kineticTheoryModel::~kineticTheoryModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::kineticTheoryModel::read()
+{
+    if
+    (
+        eddyViscosity
+        <
+            RASModel<PhaseIncompressibleTurbulenceModel<phaseModel> >
+        >::read()
+    )
+    {
+        this->coeffDict().lookup("equilibrium") >> equilibrium_;
+        e_.readIfPresent(this->coeffDict());
+        alphaMax_.readIfPresent(this->coeffDict());
+        alphaMinFriction_.readIfPresent(this->coeffDict());
+
+        viscosityModel_->read();
+        conductivityModel_->read();
+        radialModel_->read();
+        granularPressureModel_->read();
+        frictionalStressModel_->read();
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::k() const
+{
+    notImplemented("kineticTheoryModel::k()");
+    return nut_;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::epsilon() const
+{
+    notImplemented("kineticTheoryModel::epsilon()");
+    return nut_;
+}
+
+
+Foam::tmp<Foam::volSymmTensorField> Foam::kineticTheoryModel::R() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                IOobject::groupName("R", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+          - (this->nut_)*dev(twoSymm(fvc::grad(this->U_)))
+          - (lambda_*fvc::div(this->phi_))*symmTensor::I
+        )
+    );
+}
+
+
+/*
+Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::pp() const
+{
+
+    // Particle pressure coefficient
+    // Coefficient in front of Theta (Eq. 3.22, p. 45)
+    volScalarField PsCoeff
+    (
+        granularPressureModel_->granularPressureCoeff
+        (
+            alpha,
+            gs0,
+            rho,
+            e_
+        )
+    );
+
+    // Frictional pressure
+    volScalarField pf
+    (
+        frictionalStressModel_->frictionalPressure
+        (
+            alpha,
+            alphaMinFriction_,
+            alphaMax_
+        )
+    );
+
+    // Return total particle pressure
+    return PsCoeff*Theta_ + pf;
+}
+*/
+
+
+Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::pPrime() const
+{
+    // Local references
+    const volScalarField& alpha = this->alpha_;
+    const volScalarField& rho = phase_.rho();
+
+    return
+    (
+        Theta_
+       *granularPressureModel_->granularPressureCoeffPrime
+        (
+            alpha,
+            radialModel_->g0(alpha, alphaMinFriction_, alphaMax_),
+            radialModel_->g0prime(alpha, alphaMinFriction_, alphaMax_),
+            rho,
+            e_
+        )
+     +  frictionalStressModel_->frictionalPressurePrime
+        (
+            alpha,
+            alphaMinFriction_,
+            alphaMax_
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::kineticTheoryModel::pPrimef() const
+{
+    // Local references
+    const volScalarField& alpha = this->alpha_;
+    const volScalarField& rho = phase_.rho();
+
+    return fvc::interpolate
+    (
+        Theta_
+       *granularPressureModel_->granularPressureCoeffPrime
+        (
+            alpha,
+            radialModel_->g0(alpha, alphaMinFriction_, alphaMax_),
+            radialModel_->g0prime(alpha, alphaMinFriction_, alphaMax_),
+            rho,
+            e_
+        )
+     +  frictionalStressModel_->frictionalPressurePrime
+        (
+            alpha,
+            alphaMinFriction_,
+            alphaMax_
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volSymmTensorField> Foam::kineticTheoryModel::devRhoReff() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                IOobject::groupName("devRhoReff", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+          - (this->rho_*this->nut_)
+           *dev(twoSymm(fvc::grad(this->U_)))
+          - ((this->rho_*lambda_)*fvc::div(this->phi_))*symmTensor::I
+        )
+    );
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix> Foam::kineticTheoryModel::divDevRhoReff
+(
+    volVectorField& U
+) const
+{
+    return
+    (
+      - fvm::laplacian(this->rho_*this->nut_, U)
+      - fvc::div
+        (
+            (this->rho_*this->nut_)*dev2(T(fvc::grad(U)))
+          + ((this->rho_*lambda_)*fvc::div(this->phi_))
+           *dimensioned<symmTensor>("I", dimless, symmTensor::I)
+        )
+    );
+}
+
+
 void Foam::kineticTheoryModel::correct()
 {
     // Local references
@@ -343,202 +541,4 @@ void Foam::kineticTheoryModel::correct()
     }
 }
 
-
-/*
-Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::pp() const
-{
-
-    // Particle pressure coefficient
-    // Coefficient in front of Theta (Eq. 3.22, p. 45)
-    volScalarField PsCoeff
-    (
-        granularPressureModel_->granularPressureCoeff
-        (
-            alpha,
-            gs0,
-            rho,
-            e_
-        )
-    );
-
-    // Frictional pressure
-    volScalarField pf
-    (
-        frictionalStressModel_->frictionalPressure
-        (
-            alpha,
-            alphaMinFriction_,
-            alphaMax_
-        )
-    );
-
-    // Return total particle pressure
-    return PsCoeff*Theta_ + pf;
-}
-*/
-
-
-Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::pPrime() const
-{
-    // Local references
-    const volScalarField& alpha = this->alpha_;
-    const volScalarField& rho = phase_.rho();
-
-    return
-    (
-        Theta_
-       *granularPressureModel_->granularPressureCoeffPrime
-        (
-            alpha,
-            radialModel_->g0(alpha, alphaMinFriction_, alphaMax_),
-            radialModel_->g0prime(alpha, alphaMinFriction_, alphaMax_),
-            rho,
-            e_
-        )
-     +  frictionalStressModel_->frictionalPressurePrime
-        (
-            alpha,
-            alphaMinFriction_,
-            alphaMax_
-        )
-    );
-}
-
-
-Foam::tmp<Foam::surfaceScalarField> Foam::kineticTheoryModel::pPrimef() const
-{
-    // Local references
-    const volScalarField& alpha = this->alpha_;
-    const volScalarField& rho = phase_.rho();
-
-    return fvc::interpolate
-    (
-        Theta_
-       *granularPressureModel_->granularPressureCoeffPrime
-        (
-            alpha,
-            radialModel_->g0(alpha, alphaMinFriction_, alphaMax_),
-            radialModel_->g0prime(alpha, alphaMinFriction_, alphaMax_),
-            rho,
-            e_
-        )
-     +  frictionalStressModel_->frictionalPressurePrime
-        (
-            alpha,
-            alphaMinFriction_,
-            alphaMax_
-        )
-    );
-}
-
-
-void Foam::kineticTheoryModel::correctNut()
-{}
-
-
-Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::k() const
-{
-    notImplemented("kineticTheoryModel::k()");
-    return nut_;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::epsilon() const
-{
-    notImplemented("kineticTheoryModel::epsilon()");
-    return nut_;
-}
-
-
-Foam::tmp<Foam::volSymmTensorField> Foam::kineticTheoryModel::R() const
-{
-    return tmp<volSymmTensorField>
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                IOobject::groupName("R", this->U_.group()),
-                this->runTime_.timeName(),
-                this->mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-          - (this->nut_)*dev(twoSymm(fvc::grad(this->U_)))
-          - (lambda_*fvc::div(this->phi_))*symmTensor::I
-        )
-    );
-}
-
-
-Foam::tmp<Foam::volSymmTensorField> Foam::kineticTheoryModel::devRhoReff() const
-{
-    return tmp<volSymmTensorField>
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                IOobject::groupName("devRhoReff", this->U_.group()),
-                this->runTime_.timeName(),
-                this->mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-          - (this->rho_*this->nut_)
-           *dev(twoSymm(fvc::grad(this->U_)))
-          - ((this->rho_*lambda_)*fvc::div(this->phi_))*symmTensor::I
-        )
-    );
-}
-
-
-Foam::tmp<Foam::fvVectorMatrix> Foam::kineticTheoryModel::divDevRhoReff
-(
-    volVectorField& U
-) const
-{
-    return
-    (
-      - fvm::laplacian(this->rho_*this->nut_, U)
-      - fvc::div
-        (
-            (this->rho_*this->nut_)*dev2(T(fvc::grad(U)))
-          + ((this->rho_*lambda_)*fvc::div(this->phi_))
-           *dimensioned<symmTensor>("I", dimless, symmTensor::I)
-        )
-    );
-}
-
-
-bool Foam::kineticTheoryModel::read()
-{
-    if
-    (
-        eddyViscosity
-        <
-            RASModel<PhaseIncompressibleTurbulenceModel<phaseModel> >
-        >::read()
-    )
-    {
-        this->coeffDict().lookup("equilibrium") >> equilibrium_;
-        e_.readIfPresent(this->coeffDict());
-        alphaMax_.readIfPresent(this->coeffDict());
-        alphaMinFriction_.readIfPresent(this->coeffDict());
-
-        viscosityModel_->read();
-        conductivityModel_->read();
-        radialModel_->read();
-        granularPressureModel_->read();
-        frictionalStressModel_->read();
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H
index 981aaf5a68..4f6164682f 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H
@@ -25,6 +25,17 @@ Class
     Foam::kineticTheoryModel
 
 Description
+    Kinetic theory particle phase RAS model
+
+    Reference:
+    \verbatim
+        "Derivation, implementation, and validation of computer simulation
+        models for gas-solid fluidized beds",
+        B.G.M. van Wachem,
+        Ph.D. Thesis, Delft University of Technology, Amsterdam, 2000.
+    \endverbatim
+
+    There are no default model coefficients.
 
 SourceFiles
     kineticTheoryModel.C
@@ -118,6 +129,9 @@ class kineticTheoryModel
 
     // Private Member Functions
 
+        void correctNut()
+        {}
+
         //- Disallow default bitwise copy construct
         kineticTheoryModel(const kineticTheoryModel&);
 
@@ -125,13 +139,6 @@ class kineticTheoryModel
         void operator=(const kineticTheoryModel&);
 
 
-protected:
-
-    // Protected member functions
-
-        virtual void correctNut();
-
-
 public:
 
     //- Runtime type information
@@ -160,6 +167,9 @@ public:
 
     // Member Functions
 
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
         //- Return the effective viscosity
         virtual tmp<volScalarField> nuEff() const
         {
@@ -197,9 +207,6 @@ public:
 
         //- Solve the kinetic theory equations and correct the viscosity
         virtual void correct();
-
-        //- Re-read model coefficients if they have changed
-        virtual bool read();
 };
 
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phaseIncompressibleTurbulenceModels.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phaseIncompressibleTurbulenceModels.C
index aa49f6d1eb..a1817fefb5 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phaseIncompressibleTurbulenceModels.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phaseIncompressibleTurbulenceModels.C
@@ -182,7 +182,10 @@ namespace Foam
 
 #include "LESModel.H"
 #include "Smagorinsky.H"
+#include "SmagorinskyZhang.H"
 #include "kEqn.H"
+#include "NicenoKEqn.H"
+#include "continuousGasKEqn.H"
 
 namespace Foam
 {
@@ -215,6 +218,21 @@ namespace Foam
         );
     }
 
+    namespace LESModels
+    {
+        typedef SmagorinskyZhang<incompressibleTransportTurbulenceModel>
+            incompressibleSmagorinskyZhang;
+
+        defineNamedTemplateTypeNameAndDebug(incompressibleSmagorinskyZhang, 0);
+
+        addToRunTimeSelectionTable
+        (
+            incompressibleLESModel,
+            incompressibleSmagorinskyZhang,
+            dictionary
+        );
+    }
+
     namespace LESModels
     {
         typedef kEqn<incompressibleTransportTurbulenceModel>
@@ -229,6 +247,36 @@ namespace Foam
             dictionary
         );
     }
+
+    namespace LESModels
+    {
+        typedef NicenoKEqn<incompressibleTransportTurbulenceModel>
+            incompressibleNicenoKEqn;
+
+        defineNamedTemplateTypeNameAndDebug(incompressibleNicenoKEqn, 0);
+
+        addToRunTimeSelectionTable
+        (
+            incompressibleLESModel,
+            incompressibleNicenoKEqn,
+            dictionary
+        );
+    }
+
+    namespace LESModels
+    {
+        typedef continuousGasKEqn<incompressibleTransportTurbulenceModel>
+            incompressiblecontinuousGasKEqn;
+
+        defineNamedTemplateTypeNameAndDebug(incompressiblecontinuousGasKEqn, 0);
+
+        addToRunTimeSelectionTable
+        (
+            incompressibleLESModel,
+            incompressiblecontinuousGasKEqn,
+            dictionary
+        );
+    }
 }
 
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
index 6c0d8852db..dda70d3714 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
@@ -65,6 +65,11 @@ Foam::phasePressureModel::phasePressureModel
     )
 {
     this->nut_ == dimensionedScalar("zero", this->nut_.dimensions(), 0.0);
+
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
 }
 
 
@@ -76,35 +81,27 @@ Foam::phasePressureModel::~phasePressureModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::phasePressureModel::correct()
-{}
-
-
-void Foam::phasePressureModel::correctNut()
-{}
-
-
-Foam::tmp<Foam::volScalarField> Foam::phasePressureModel::pPrime() const
+bool Foam::phasePressureModel::read()
 {
-    return
-        g0_
-       *min
-        (
-            exp(preAlphaExp_*(this->alpha_ - alphaMax_)),
-            expMax_
-        );
-}
+    if
+    (
+        eddyViscosity
+        <
+            RASModel<PhaseIncompressibleTurbulenceModel<phaseModel> >
+        >::read()
+    )
+    {
+        this->coeffDict().lookup("alphaMax") >> alphaMax_;
+        this->coeffDict().lookup("preAlphaExp") >> preAlphaExp_;
+        this->coeffDict().lookup("expMax") >> expMax_;
+        g0_.readIfPresent(this->coeffDict());
 
-
-Foam::tmp<Foam::surfaceScalarField> Foam::phasePressureModel::pPrimef() const
-{
-    return
-        g0_
-       *min
-        (
-            exp(preAlphaExp_*(fvc::interpolate(this->alpha_) - alphaMax_)),
-            expMax_
-        );
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
@@ -148,6 +145,30 @@ Foam::tmp<Foam::volSymmTensorField> Foam::phasePressureModel::R() const
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::phasePressureModel::pPrime() const
+{
+    return
+        g0_
+       *min
+        (
+            exp(preAlphaExp_*(this->alpha_ - alphaMax_)),
+            expMax_
+        );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::phasePressureModel::pPrimef() const
+{
+    return
+        g0_
+       *min
+        (
+            exp(preAlphaExp_*(fvc::interpolate(this->alpha_) - alphaMax_)),
+            expMax_
+        );
+}
+
+
 Foam::tmp<Foam::volSymmTensorField> Foam::phasePressureModel::devRhoReff() const
 {
     return tmp<volSymmTensorField>
@@ -190,28 +211,8 @@ Foam::tmp<Foam::fvVectorMatrix> Foam::phasePressureModel::divDevRhoReff
 }
 
 
-bool Foam::phasePressureModel::read()
-{
-    if
-    (
-        eddyViscosity
-        <
-            RASModel<PhaseIncompressibleTurbulenceModel<phaseModel> >
-        >::read()
-    )
-    {
-        this->coeffDict().lookup("alphaMax") >> alphaMax_;
-        this->coeffDict().lookup("preAlphaExp") >> preAlphaExp_;
-        this->coeffDict().lookup("expMax") >> expMax_;
-        g0_.readIfPresent(this->coeffDict());
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
+void Foam::phasePressureModel::correct()
+{}
 
 
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.H
index 2441094639..aebdf09014 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.H
@@ -25,6 +25,23 @@ Class
     Foam::phasePressureModel
 
 Description
+    Particle-particle phase-pressure RAS model
+
+    The derivative of the phase-pressure with respect to the phase-fraction
+    is evaluated as
+
+        g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        phasePressureCoeffs
+        {
+            preAlphaExp     500;
+            expMax          1000;
+            alphaMax        0.62;
+            g0              1000;
+        }
+    \endverbatim
 
 SourceFiles
     phasePressureModel.C
@@ -39,7 +56,6 @@ SourceFiles
 #include "PhaseIncompressibleTurbulenceModel.H"
 #include "phaseModel.H"
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -80,6 +96,9 @@ class phasePressureModel
 
     // Private Member Functions
 
+        void correctNut()
+        {}
+
         //- Disallow default bitwise copy construct
         phasePressureModel(const phasePressureModel&);
 
@@ -87,13 +106,6 @@ class phasePressureModel
         void operator=(const phasePressureModel&);
 
 
-protected:
-
-    // Protected member functions
-
-        virtual void correctNut();
-
-
 public:
 
     //- Runtime type information
@@ -122,6 +134,9 @@ public:
 
     // Member Functions
 
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
         //- Return the effective viscosity
         virtual tmp<volScalarField> nuEff() const
         {
@@ -159,9 +174,6 @@ public:
 
         //- Solve the kinetic theory equations and correct the viscosity
         virtual void correct();
-
-        //- Re-read model coefficients if they have changed
-        virtual bool read();
 };
 
 
diff --git a/src/TurbulenceModels/compressible/compressibleTurbulenceModel.H b/src/TurbulenceModels/compressible/compressibleTurbulenceModel.H
index d762dc5a05..0754085f4a 100644
--- a/src/TurbulenceModels/compressible/compressibleTurbulenceModel.H
+++ b/src/TurbulenceModels/compressible/compressibleTurbulenceModel.H
@@ -96,7 +96,7 @@ public:
     {}
 
 
-    // Member Functions
+    // Member functions
 
         //- Return the effective stress tensor including the laminar stress
         virtual tmp<volSymmTensorField> devRhoReff() const = 0;
diff --git a/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.H b/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.H
index 8766fe454f..a73255b3e0 100644
--- a/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.H
+++ b/src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.H
@@ -28,6 +28,31 @@ Group
     grpRASTurbulence
 
 Description
+    Continuous-phase k-epsilon model including bubble-generated turbulence.
+
+    Reference:
+    \verbatim
+        "The simulation of multidimensional multiphase flows",
+        Lahey R.T.,
+        Nucl. Eng. & Design
+        2005 (235) pp.1043-1060.
+    \endverbatim
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        LaheyKEpsilonCoeffs
+        {
+            Cmu             0.09;
+            C1              1.44;
+            C2              1.92;
+            C3              -0.33;
+            sigmak          1.0;
+            sigmaEps        1.3;
+            Cp              0.25;
+            Cmub            0.6;
+            alphaInversion  0.3;
+        }
+    \endverbatim
 
 SourceFiles
     LaheyKEpsilon.C
@@ -63,6 +88,20 @@ class LaheyKEpsilon
         > *gasTurbulencePtr_;
 
 
+    // Private Member Functions
+
+        //- Return the turbulence model for the gas phase
+        const PhaseIncompressibleTurbulenceModel
+        <
+            typename BasicTurbulenceModel::transportModel
+        >&
+        gasTurbulence() const;
+
+        // Disallow default bitwise copy construct and assignment
+        LaheyKEpsilon(const LaheyKEpsilon&);
+        LaheyKEpsilon& operator=(const LaheyKEpsilon&);
+
+
 protected:
 
     // Protected data
@@ -74,7 +113,7 @@ protected:
             dimensionedScalar Cmub_;
 
 
-    // Protected member functions
+    // Protected Member Functions
 
         virtual void correctNut();
         tmp<volScalarField> bubbleG() const;
@@ -82,6 +121,7 @@ protected:
         virtual tmp<fvScalarMatrix> kSource() const;
         virtual tmp<fvScalarMatrix> epsilonSource() const;
 
+
 public:
 
     typedef typename BasicTurbulenceModel::alphaField alphaField;
@@ -116,15 +156,11 @@ public:
 
     // Member Functions
 
-        //- Return the turbulence model for the gas phase
-        const PhaseIncompressibleTurbulenceModel<transportModel>&
-        gasTurbulence() const;
+        //- Read model coefficients if they have changed
+        virtual bool read();
 
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
-
-        //- Read RASProperties dictionary
-        virtual bool read();
 };
 
 
diff --git a/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.H b/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.H
index ecb9486339..7cf70368ac 100644
--- a/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.H
+++ b/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.H
@@ -28,6 +28,31 @@ Group
     grpRASTurbulence
 
 Description
+    k-epsilon model for the gas-phase in a two-phase system
+    supporting phase-inversion.
+
+    In the limit that the gas-phase fraction approaches zero a contribution from
+    the other phase is blended into the k and epsilon equations up to the
+    phase-fraction of alphaInversion at which point phase-inversion is
+    considered to have occurred and the model reverts to the pure single-phase
+    form.
+
+    This model is unpublished and is provided as a stable numerical framework
+    on which a more physical model may be built.
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        continuousGasKEpsilonCoeffs
+        {
+            Cmu             0.09;
+            C1              1.44;
+            C2              1.92;
+            C3              -0.33;
+            sigmak          1.0;
+            sigmaEps        1.3;
+            alphaInversion  0.7;
+        }
+    \endverbatim
 
 SourceFiles
     continuousGasKEpsilon.C
@@ -62,6 +87,13 @@ class continuousGasKEpsilon
         volScalarField nutEff_;
 
 
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        continuousGasKEpsilon(const continuousGasKEpsilon&);
+        continuousGasKEpsilon& operator=(const continuousGasKEpsilon&);
+
+
 protected:
 
     // Protected data
@@ -71,7 +103,7 @@ protected:
             dimensionedScalar alphaInversion_;
 
 
-    // Protected member functions
+    // Protected Member Functions
 
         virtual void correctNut();
         tmp<volScalarField> phaseTransferCoeff() const;
@@ -113,6 +145,9 @@ public:
 
     // Member Functions
 
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
         //- Return the turbulence model for the liquid phase
         const turbulenceModel& liquidTurbulence() const;
 
@@ -124,9 +159,6 @@ public:
 
         //- Return the Reynolds stress tensor
         virtual tmp<volSymmTensorField> R() const;
-
-        //- Read RASProperties dictionary
-        virtual bool read();
 };
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.H b/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.H
index 6c6bb36d31..3fc24b426c 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.H
+++ b/src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.H
@@ -168,6 +168,10 @@ public:
 
     // Member Functions
 
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+
         // Access
 
             //- Const access to the coefficients dictionary
@@ -216,9 +220,6 @@ public:
 
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
-
-        //- Read LESProperties dictionary
-        virtual bool read();
 };
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/LES/Smagorinsky/Smagorinsky.C b/src/TurbulenceModels/turbulenceModels/LES/Smagorinsky/Smagorinsky.C
index 5ed7bbc96f..da0db918a1 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/Smagorinsky/Smagorinsky.C
+++ b/src/TurbulenceModels/turbulenceModels/LES/Smagorinsky/Smagorinsky.C
@@ -59,21 +59,21 @@ Smagorinsky<BasicTurbulenceModel>::Smagorinsky
         propertiesName
     ),
 
-    ck_
+    Ck_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "ck",
+            "Ck",
             this->coeffDict_,
-            0.02
+            0.094
         )
     ),
 
-    ce_
+    Ce_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "ce",
+            "Ce",
             this->coeffDict_,
             1.048
         )
@@ -94,8 +94,8 @@ bool Smagorinsky<BasicTurbulenceModel>::read()
 {
     if (eddyViscosity<LESModel<BasicTurbulenceModel> >::read())
     {
-        ck_.readIfPresent(this->coeffDict());
-        ce_.readIfPresent(this->coeffDict());
+        Ck_.readIfPresent(this->coeffDict());
+        Ce_.readIfPresent(this->coeffDict());
 
         return true;
     }
@@ -105,6 +105,7 @@ bool Smagorinsky<BasicTurbulenceModel>::read()
     }
 }
 
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> Smagorinsky<BasicTurbulenceModel>::k
 (
@@ -113,9 +114,9 @@ tmp<volScalarField> Smagorinsky<BasicTurbulenceModel>::k
 {
     volSymmTensorField D(symm(gradU));
 
-    volScalarField a(ce_/this->delta());
+    volScalarField a(Ce_/this->delta());
     volScalarField b((2.0/3.0)*tr(D));
-    volScalarField c(2*ck_*this->delta()*(dev(D) && D));
+    volScalarField c(2*Ck_*this->delta()*(dev(D) && D));
 
     return sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a));
 }
@@ -136,7 +137,7 @@ tmp<volScalarField> Smagorinsky<BasicTurbulenceModel>::epsilon() const
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            ce_*k()*sqrt(k())/this->delta()
+            Ce_*k()*sqrt(k())/this->delta()
         )
     );
 }
@@ -147,7 +148,7 @@ void Smagorinsky<BasicTurbulenceModel>::correctNut()
 {
     volScalarField k(this->k(fvc::grad(this->U_)));
 
-    this->nut_ = ck_*this->delta()*sqrt(k);
+    this->nut_ = Ck_*this->delta()*sqrt(k);
     this->nut_.correctBoundaryConditions();
 }
 
diff --git a/src/TurbulenceModels/turbulenceModels/LES/Smagorinsky/Smagorinsky.H b/src/TurbulenceModels/turbulenceModels/LES/Smagorinsky/Smagorinsky.H
index 8467525c85..c74278dee9 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/Smagorinsky/Smagorinsky.H
+++ b/src/TurbulenceModels/turbulenceModels/LES/Smagorinsky/Smagorinsky.H
@@ -28,7 +28,7 @@ Group
     grpLESTurbulence
 
 Description
-    The Smagorinsky Model.
+    The Smagorinsky SGS model.
 
     Algebraic eddy viscosity SGS model founded on the assumption that
     local equilibrium prevails.
@@ -40,8 +40,17 @@ Description
     where
 
         D = symm(grad(U));
-        k from D:B + ce*k^3/2/delta = 0
-        nuSgs = ck*sqrt(k)*delta
+        k from D:B + Ce*k^3/2/delta = 0
+        nuSgs = Ck*sqrt(k)*delta
+    \endverbatim
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        SmagorinskyCoeffs
+        {
+            Ck                  0.094;
+            Ce                  1.048;
+        }
     \endverbatim
 
 SourceFiles
@@ -71,26 +80,22 @@ class Smagorinsky
 :
     public eddyViscosity<LESModel<BasicTurbulenceModel> >
 {
-
-protected:
-
-    // Protected data
-
-        dimensionedScalar ck_;
-        dimensionedScalar ce_;
-
-
     // Private Member Functions
 
-        //- Update sub-grid scale fields
-        void updateSubGridScaleFields(const volTensorField& gradU);
-
         // Disallow default bitwise copy construct and assignment
         Smagorinsky(const Smagorinsky&);
         Smagorinsky& operator=(const Smagorinsky&);
 
 
-    // Protected member functions
+protected:
+
+    // Protected data
+
+        dimensionedScalar Ck_;
+        dimensionedScalar Ce_;
+
+
+    // Protected Member Functions
 
         //- Return SGS kinetic energy
         //  calculated from the given velocity gradient
@@ -133,6 +138,9 @@ public:
 
     // Member Functions
 
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
         //- Return SGS kinetic energy
         virtual tmp<volScalarField> k() const
         {
@@ -144,9 +152,6 @@ public:
 
         //- Correct Eddy-Viscosity and related properties
         virtual void correct();
-
-        //- Read model coefficients if they have changed
-        virtual bool read();
 };
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/LES/kEqn/kEqn.C b/src/TurbulenceModels/turbulenceModels/LES/kEqn/kEqn.C
index 704c8e4007..dd99d718f9 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/kEqn/kEqn.C
+++ b/src/TurbulenceModels/turbulenceModels/LES/kEqn/kEqn.C
@@ -72,21 +72,21 @@ kEqn<BasicTurbulenceModel>::kEqn
         this->mesh_
     ),
 
-    ck_
+    Ck_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "ck",
+            "Ck",
             this->coeffDict_,
             0.094
         )
     ),
 
-    ce_
+    Ce_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "ce",
+            "Ce",
             this->coeffDict_,
             1.048
         )
@@ -107,8 +107,8 @@ bool kEqn<BasicTurbulenceModel>::read()
 {
     if (eddyViscosity<LESModel<BasicTurbulenceModel> >::read())
     {
-        ck_.readIfPresent(this->coeffDict());
-        ce_.readIfPresent(this->coeffDict());
+        Ck_.readIfPresent(this->coeffDict());
+        Ce_.readIfPresent(this->coeffDict());
 
         return true;
     }
@@ -134,7 +134,7 @@ tmp<volScalarField> kEqn<BasicTurbulenceModel>::epsilon() const
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            ce_*k()*sqrt(k())/this->delta()
+            Ce_*k()*sqrt(k())/this->delta()
         )
     );
 }
@@ -143,7 +143,7 @@ tmp<volScalarField> kEqn<BasicTurbulenceModel>::epsilon() const
 template<class BasicTurbulenceModel>
 void kEqn<BasicTurbulenceModel>::correctNut()
 {
-    this->nut_ = ck_*sqrt(k_)*this->delta();
+    this->nut_ = Ck_*sqrt(k_)*this->delta();
     this->nut_.correctBoundaryConditions();
 }
 
@@ -196,7 +196,7 @@ void kEqn<BasicTurbulenceModel>::correct()
      ==
         alpha*rho*G
       - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
-      - fvm::Sp(ce_*alpha*rho*sqrt(k_)/this->delta(), k_)
+      - fvm::Sp(Ce_*alpha*rho*sqrt(k_)/this->delta(), k_)
       + kSource()
     );
 
diff --git a/src/TurbulenceModels/turbulenceModels/LES/kEqn/kEqn.H b/src/TurbulenceModels/turbulenceModels/LES/kEqn/kEqn.H
index ac065c0be7..838655d00a 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/kEqn/kEqn.H
+++ b/src/TurbulenceModels/turbulenceModels/LES/kEqn/kEqn.H
@@ -35,7 +35,7 @@ Description
     \verbatim
         d/dt(rho*k) + div(rho*U*k) - div(rho*nuEff*grad(k))
         =
-        -rho*D:B - ce*rho*k^(3/2)/delta
+        -rho*D:B - Ce*rho*k^(3/2)/delta
 
     and
 
@@ -44,10 +44,19 @@ Description
     where
 
         D = symm(grad(U));
-        nuSgs = ck*sqrt(k)*delta
+        nuSgs = Ck*sqrt(k)*delta
         nuEff = nuSgs + nu
     \endverbatim
 
+    The default model coefficients correspond to the following:
+    \verbatim
+        NicenoKEqnCoeffs
+        {
+            Ck                  0.094;
+            Ce                  1.048;
+        }
+    \endverbatim
+
 SourceFiles
     kEqn.C
 
@@ -75,6 +84,12 @@ class kEqn
 :
     public eddyViscosity<LESModel<BasicTurbulenceModel> >
 {
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        kEqn(const kEqn&);
+        kEqn& operator=(const kEqn&);
+
 
 protected:
 
@@ -82,21 +97,11 @@ protected:
 
         volScalarField k_;
 
-        dimensionedScalar ck_;
-        dimensionedScalar ce_;
+        dimensionedScalar Ck_;
+        dimensionedScalar Ce_;
 
 
-    // Private Member Functions
-
-        //- Update sub-grid scale fields
-        void updateSubGridScaleFields();
-
-        // Disallow default bitwise copy construct and assignment
-        kEqn(const kEqn&);
-        kEqn& operator=(const kEqn&);
-
-
-    // Protected member functions
+    // Protected Member Functions
 
         virtual void correctNut();
         virtual tmp<fvScalarMatrix> kSource() const;
@@ -136,6 +141,9 @@ public:
 
     // Member Functions
 
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
         //- Return SGS kinetic energy
         virtual tmp<volScalarField> k() const
         {
@@ -156,9 +164,6 @@ public:
 
         //- Correct Eddy-Viscosity and related properties
         virtual void correct();
-
-        //- Read model coefficients if they have changed
-        virtual bool read();
 };
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/RASModel/RASModel.H b/src/TurbulenceModels/turbulenceModels/RAS/RASModel/RASModel.H
index 9185de493b..6e1636d0ae 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/RASModel/RASModel.H
+++ b/src/TurbulenceModels/turbulenceModels/RAS/RASModel/RASModel.H
@@ -170,6 +170,10 @@ public:
 
     // Member Functions
 
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+
         // Access
 
             //- Return the lower allowable limit for k (default: SMALL)
@@ -236,9 +240,6 @@ public:
 
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
-
-        //- Read RASProperties dictionary
-        virtual bool read();
 };
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.H b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.H
index bc840a26e3..c30004efc3 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.H
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.H
@@ -70,6 +70,12 @@ class kEpsilon
 :
     public eddyViscosity<RASModel<BasicTurbulenceModel> >
 {
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        kEpsilon(const kEpsilon&);
+        kEpsilon& operator=(const kEpsilon&);
+
 
 protected:
 
@@ -90,7 +96,7 @@ protected:
             volScalarField epsilon_;
 
 
-    // Protected member functions
+    // Protected Member Functions
 
         virtual void correctNut();
         virtual tmp<fvScalarMatrix> kSource() const;
@@ -131,6 +137,9 @@ public:
 
     // Member Functions
 
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
         //- Return the effective diffusivity for k
         tmp<volScalarField> DkEff() const
         {
@@ -171,9 +180,6 @@ public:
 
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
-
-        //- Re-read model coefficients if they have changed
-        virtual bool read();
 };
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/eddyViscosity/eddyViscosity.H b/src/TurbulenceModels/turbulenceModels/eddyViscosity/eddyViscosity.H
index 19fcaaa193..59b0c61c0f 100644
--- a/src/TurbulenceModels/turbulenceModels/eddyViscosity/eddyViscosity.H
+++ b/src/TurbulenceModels/turbulenceModels/eddyViscosity/eddyViscosity.H
@@ -62,7 +62,7 @@ protected:
             volScalarField nut_;
 
 
-    // Protected member functions
+    // Protected Member Functions
 
         virtual void correctNut() = 0;
 
@@ -97,6 +97,9 @@ public:
 
     // Member Functions
 
+        //- Re-read model coefficients if they have changed
+        virtual bool read() = 0;
+
         //- Return the turbulence viscosity
         virtual tmp<volScalarField> nut() const
         {
@@ -123,9 +126,6 @@ public:
 
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct() = 0;
-
-        //- Re-read model coefficients if they have changed
-        virtual bool read() = 0;
 };
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/turbulenceModel.H b/src/TurbulenceModels/turbulenceModels/turbulenceModel.H
index 31385546e4..df10e3b081 100644
--- a/src/TurbulenceModels/turbulenceModels/turbulenceModel.H
+++ b/src/TurbulenceModels/turbulenceModels/turbulenceModel.H
@@ -114,6 +114,9 @@ public:
 
     // Member Functions
 
+        //- Read model coefficients if they have changed
+        virtual bool read() = 0;
+
         const Time& time() const
         {
             return runTime_;
@@ -205,9 +208,6 @@ public:
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct() = 0;
 
-        //- Read LESProperties or RASProperties dictionary
-        virtual bool read() = 0;
-
         //- Default dummy write function
         virtual bool writeData(Ostream&) const
         {