diff --git a/src/TurbulenceModels/phaseCompressible/RAS/kOmegaSSTSato/kOmegaSSTSato.H b/src/TurbulenceModels/phaseCompressible/RAS/kOmegaSSTSato/kOmegaSSTSato.H
index e25ba65fd7..4a95eb1aba 100644
--- a/src/TurbulenceModels/phaseCompressible/RAS/kOmegaSSTSato/kOmegaSSTSato.H
+++ b/src/TurbulenceModels/phaseCompressible/RAS/kOmegaSSTSato/kOmegaSSTSato.H
@@ -150,6 +150,7 @@ protected:
     // Protected data
 
         // Model coefficients
+
             dimensionedScalar Cmub_;
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C b/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C
index dbe579d0cf..de51e92981 100644
--- a/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C
+++ b/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C
@@ -135,13 +135,25 @@ void kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::correctNut()
 
 
 template<class TurbulenceModel, class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::epsilonByk
+tmp<volScalarField::Internal>
+kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::Pk
 (
-    const volScalarField& F1,
-    const volScalarField& F2
+    const volScalarField::Internal& G
 ) const
 {
-    return betaStar_*omega_;
+    return min(G, (c1_*betaStar_)*this->k_()*this->omega_());
+}
+
+
+template<class TurbulenceModel, class BasicTurbulenceModel>
+tmp<volScalarField::Internal>
+kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::epsilonByk
+(
+    const volScalarField::Internal& F1,
+    const volScalarField::Internal& F2
+) const
+{
+    return betaStar_*omega_();
 }
 
 
@@ -178,9 +190,9 @@ kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::omegaSource() const
 template<class TurbulenceModel, class BasicTurbulenceModel>
 tmp<fvScalarMatrix> kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::Qsas
 (
-    const volScalarField& S2,
-    const volScalarField& gamma,
-    const volScalarField& beta
+    const volScalarField::Internal& S2,
+    const volScalarField::Internal& gamma,
+    const volScalarField::Internal& beta
 ) const
 {
     return tmp<fvScalarMatrix>
@@ -419,12 +431,15 @@ void kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::correct()
 
     BasicTurbulenceModel::correct();
 
-    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
+    volScalarField::Internal 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);
+    volScalarField::Internal GbyNu(dev(twoSymm(tgradU()())) && tgradU()());
+    volScalarField::Internal G(this->GName(), nut()*GbyNu);
     tgradU.clear();
 
     // Update omega and G at the wall
@@ -439,8 +454,8 @@ void kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::correct()
     volScalarField F23(this->F23());
 
     {
-        volScalarField gamma(this->gamma(F1));
-        volScalarField beta(this->beta(F1));
+        volScalarField::Internal gamma(this->gamma(F1));
+        volScalarField::Internal beta(this->beta(F1));
 
         // Turbulent frequency equation
         tmp<fvScalarMatrix> omegaEqn
@@ -449,20 +464,21 @@ void kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::correct()
           + fvm::div(alphaRhoPhi, omega_)
           - fvm::laplacian(alpha*rho*DomegaEff(F1), omega_)
          ==
-            alpha*rho*gamma
+            alpha()*rho()*gamma
            *min
             (
                 GbyNu,
-                (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23*sqrt(S2))
+                (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((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_,
+                alpha()*rho()*(F1() - scalar(1))*CDkOmega()/omega_(),
                 omega_
             )
-          + Qsas(S2, gamma, beta)
+          + Qsas(S2(), gamma, beta)
           + omegaSource()
           + fvOptions(alpha, rho, omega_)
         );
@@ -482,9 +498,9 @@ void kOmegaSST<TurbulenceModel, BasicTurbulenceModel>::correct()
       + 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_)
+        alpha()*rho()*Pk(G)
+      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
+      - fvm::Sp(alpha()*rho()*epsilonByk(F1, F23), k_)
       + kSource()
       + fvOptions(alpha, rho, k_)
     );
diff --git a/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.H b/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.H
index 5adf9961ab..72301b6792 100644
--- a/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.H
+++ b/src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.H
@@ -162,12 +162,12 @@ protected:
             volScalarField omega_;
 
 
-    // Private Member Functions
+    // Protected Member Functions
 
-        tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
-        tmp<volScalarField> F2() const;
-        tmp<volScalarField> F3() const;
-        tmp<volScalarField> F23() const;
+        virtual tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
+        virtual tmp<volScalarField> F2() const;
+        virtual tmp<volScalarField> F3() const;
+        virtual tmp<volScalarField> F23() const;
 
         tmp<volScalarField> blend
         (
@@ -179,6 +179,16 @@ protected:
             return F1*(psi1 - psi2) + psi2;
         }
 
+        tmp<volScalarField::Internal> blend
+        (
+            const volScalarField::Internal& 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_);
@@ -189,28 +199,37 @@ protected:
             return blend(F1, alphaOmega1_, alphaOmega2_);
         }
 
-        tmp<volScalarField> beta(const volScalarField& F1) const
+        tmp<volScalarField::Internal> beta
+        (
+            const volScalarField::Internal& F1
+        ) const
         {
             return blend(F1, beta1_, beta2_);
         }
 
-        tmp<volScalarField> gamma(const volScalarField& F1) const
+        tmp<volScalarField::Internal> gamma
+        (
+            const volScalarField::Internal& 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
+        //- Return k production rate
+        virtual tmp<volScalarField::Internal> Pk
         (
-            const volScalarField& F1,
-            const volScalarField& F2
+            const volScalarField::Internal& G
+        ) const;
+
+        //- Return epsilon/k which for standard RAS is betaStar*omega
+        virtual tmp<volScalarField::Internal> epsilonByk
+        (
+            const volScalarField::Internal& F1,
+            const volScalarField::Internal& F2
         ) const;
 
         virtual tmp<fvScalarMatrix> kSource() const;
@@ -219,9 +238,9 @@ protected:
 
         virtual tmp<fvScalarMatrix> Qsas
         (
-            const volScalarField& S2,
-            const volScalarField& gamma,
-            const volScalarField& beta
+            const volScalarField::Internal& S2,
+            const volScalarField::Internal& gamma,
+            const volScalarField::Internal& beta
         ) const;
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.C b/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.C
index 0ad21dd8cb..8714800784 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.C
+++ b/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.C
@@ -35,27 +35,27 @@ namespace LESModels
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::Lt() const
+tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::Lt() const
 {
-    return sqrt(this->k_)/(this->betaStar_*this->omega_);
+    return sqrt(this->k_())/(this->betaStar_*this->omega_());
 }
 
 
 template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::FDES
+tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::FDES
 (
-    const volScalarField& F1,
-    const volScalarField& F2
+    const volScalarField::Internal& F1,
+    const volScalarField::Internal& F2
 ) const
 {
     switch (FSST_)
     {
         case 0:
-            return max(Lt()/(CDES_*this->delta()), scalar(1));
+            return max(Lt()/(CDES_*this->delta()()), scalar(1));
         case 1:
-            return max(Lt()*(1 - F1)/(CDES_*this->delta()), scalar(1));
+            return max(Lt()*(1 - F1)/(CDES_*this->delta()()), scalar(1));
         case 2:
-            return max(Lt()*(1 - F2)/(CDES_*this->delta()), scalar(1));
+            return max(Lt()*(1 - F2)/(CDES_*this->delta()()), scalar(1));
         default:
             FatalErrorInFunction
                 << "Incorrect FSST = " << FSST_ << ", should be 0, 1 or 2"
@@ -66,13 +66,13 @@ tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::FDES
 
 
 template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::epsilonByk
+tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::epsilonByk
 (
-    const volScalarField& F1,
-    const volScalarField& F2
+    const volScalarField::Internal& F1,
+    const volScalarField::Internal& F2
 ) const
 {
-    return this->betaStar_*this->omega_*FDES(F1, F2);
+    return this->betaStar_*this->omega_()*FDES(F1, F2);
 }
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.H b/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.H
index 213de48d4e..3d4769d296 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.H
+++ b/src/TurbulenceModels/turbulenceModels/LES/kOmegaSSTDES/kOmegaSSTDES.H
@@ -101,21 +101,21 @@ protected:
     // Protected Member Functions
 
         //- Return the turbulent length-scale
-        tmp<volScalarField> Lt() const;
+        tmp<volScalarField::Internal> Lt() const;
 
         //- The DES dissipation-rate multiplier with options zonal filtering
         //  based on either F1 or F2
-        virtual tmp<volScalarField> FDES
+        virtual tmp<volScalarField::Internal> FDES
         (
-            const volScalarField& F1,
-            const volScalarField& F2
+            const volScalarField::Internal& F1,
+            const volScalarField::Internal& F2
         ) const;
 
         //- Return epsilon/k which for standard RAS is betaStar*omega
-        virtual tmp<volScalarField> epsilonByk
+        virtual tmp<volScalarField::Internal> epsilonByk
         (
-            const volScalarField& F1,
-            const volScalarField& F2
+            const volScalarField::Internal& F1,
+            const volScalarField::Internal& F2
         ) const;
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTSAS/kOmegaSSTSAS.C b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTSAS/kOmegaSSTSAS.C
index 08b3b5af43..5d5ce66dbb 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTSAS/kOmegaSSTSAS.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTSAS/kOmegaSSTSAS.C
@@ -37,23 +37,23 @@ namespace RASModels
 template<class BasicTurbulenceModel>
 tmp<fvScalarMatrix> kOmegaSSTSAS<BasicTurbulenceModel>::Qsas
 (
-    const volScalarField& S2,
-    const volScalarField& gamma,
-    const volScalarField& beta
+    const volScalarField::Internal& S2,
+    const volScalarField::Internal& gamma,
+    const volScalarField::Internal& beta
 ) const
 {
-    volScalarField L
+    volScalarField::Internal L
     (
-        sqrt(this->k_)/(pow025(this->betaStar_)*this->omega_)
+        sqrt(this->k_())/(pow025(this->betaStar_)*this->omega_())
     );
 
-    volScalarField Lvk
+    volScalarField::Internal Lvk
     (
         max
         (
             kappa_*sqrt(S2)
            /(
-                mag(fvc::laplacian(this->U_))
+                mag(fvc::laplacian(this->U_))()()
               + dimensionedScalar
                 (
                     "ROOTVSMALL",
@@ -61,29 +61,29 @@ tmp<fvScalarMatrix> kOmegaSSTSAS<BasicTurbulenceModel>::Qsas
                     ROOTVSMALL
                 )
             ),
-            Cs_*sqrt(kappa_*zeta2_/(beta/this->betaStar_ - gamma))*delta()
+            Cs_*sqrt(kappa_*zeta2_/(beta/this->betaStar_ - gamma))*delta()()
         )
     );
 
     return fvm::Su
     (
-        this->alpha_*this->rho_
+        this->alpha_()*this->rho_()
        *min
         (
             max
             (
                 zeta2_*kappa_*S2*sqr(L/Lvk)
-              - (2*C_/sigmaPhi_)*this->k_
+              - (2*C_/sigmaPhi_)*this->k_()
                *max
                 (
-                    magSqr(fvc::grad(this->omega_))/sqr(this->omega_),
-                    magSqr(fvc::grad(this->k_))/sqr(this->k_)
+                    magSqr(fvc::grad(this->omega_)()())/sqr(this->omega_()),
+                    magSqr(fvc::grad(this->k_)()())/sqr(this->k_())
                 ),
                 dimensionedScalar("0", dimensionSet(0, 0, -2, 0, 0), 0)
             ),
             // Limit SAS production of omega for numerical stability,
             // particularly during start-up
-            this->omega_/(0.1*this->omega_.time().deltaT())
+            this->omega_()/(0.1*this->omega_.time().deltaT())
         ),
         this->omega_
     );
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTSAS/kOmegaSSTSAS.H b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTSAS/kOmegaSSTSAS.H
index 52c918c1c7..0840a26450 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTSAS/kOmegaSSTSAS.H
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTSAS/kOmegaSSTSAS.H
@@ -131,9 +131,9 @@ protected:
         //- SAS omega source
         virtual tmp<fvScalarMatrix> Qsas
         (
-            const volScalarField& S2,
-            const volScalarField& gamma,
-            const volScalarField& beta
+            const volScalarField::Internal& S2,
+            const volScalarField::Internal& gamma,
+            const volScalarField::Internal& beta
         ) const;