diff --git a/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C b/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
index cff839b5e3714b3dcc3e6a315f50ccba72ff1b1b..609fa986f183e5db1b702c804f9c6ba7d2a3ac2e 100644
--- a/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
+++ b/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
@@ -106,7 +106,7 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
 
     volScalarField K = 0.5*tr(B_);
 
-    solve
+    tmp<fvSymmTensorMatrix> BEqn
     (
         fvm::ddt(rho(), B_)
       + fvm::div(phi(), B_)
@@ -118,6 +118,8 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
       - (2*ce_ - 0.667*cm_)*I*rho()*epsilon()
     );
 
+    BEqn().relax();
+    BEqn().solve();
 
     // Bounding the component kinetic energies
 
diff --git a/src/turbulenceModels/compressible/LES/GenEddyVisc/GenEddyVisc.C b/src/turbulenceModels/compressible/LES/GenEddyVisc/GenEddyVisc.C
index 055c76214e814589ba1dd8c461a0d6f6c53810d3..366808ad08ab9cbe1de929c3e3932df81ea33c3b 100644
--- a/src/turbulenceModels/compressible/LES/GenEddyVisc/GenEddyVisc.C
+++ b/src/turbulenceModels/compressible/LES/GenEddyVisc/GenEddyVisc.C
@@ -108,7 +108,9 @@ GenEddyVisc::GenEddyVisc
         ),
         mesh_
     )
-{}
+{
+    bound(k_, kMin_);
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
diff --git a/src/turbulenceModels/compressible/LES/LESModel/LESModel.C b/src/turbulenceModels/compressible/LES/LESModel/LESModel.C
index daeec0a05eee6d48dba7c8850e8da2d6ac7a55ea..556ee0ca51554326b4d66a0987bd4875aa5faa61 100644
--- a/src/turbulenceModels/compressible/LES/LESModel/LESModel.C
+++ b/src/turbulenceModels/compressible/LES/LESModel/LESModel.C
@@ -79,7 +79,7 @@ LESModel::LESModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
+    kMin_("kMin", sqr(dimVelocity), SMALL),
 
     delta_(LESdelta::New("delta", U.mesh(), *this))
 {
@@ -131,8 +131,12 @@ autoPtr<LESModel> LESModel::New
     {
         FatalErrorIn
         (
-            "LESModel::New(const volVectorField& U, const "
-            "surfaceScalarField& phi, const basicThermo&)"
+            "LESModel::New"
+            "("
+                "const volVectorField&, "
+                "const surfaceScalarField&, "
+                "const basicThermo&"
+            ")"
         )   << "Unknown LESModel type " << modelName
             << endl << endl
             << "Valid LESModel types are :" << endl
@@ -167,7 +171,7 @@ bool LESModel::read()
             coeffDict_ <<= *dictPtr;
         }
 
-        readIfPresent("kMin", kMin_);
+        kMin_.readIfPresent(*this);
 
         delta_().read(*this);
 
diff --git a/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C
index 22f688e7d5715a941ed8898c03c47b621a9af260..b054f2eebac1a0d79bf32699cd6f3f8a1f5b6591 100644
--- a/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C
@@ -302,7 +302,7 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& tgradU)
     volScalarField Stilda =
         fv3()*::sqrt(2.0)*mag(skew(gradU)) + fv2()*nuTilda_/sqr(kappa_*dTilda_);
 
-    solve
+    tmp<fvScalarMatrix> nuTildaEqn
     (
         fvm::ddt(rho(), nuTilda_)
       + fvm::div(phi(), nuTilda_)
@@ -318,6 +318,9 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& tgradU)
       - fvm::Sp(rho()*Cw1_*fw(Stilda)*nuTilda_/sqr(dTilda_), nuTilda_)
     );
 
+    nuTildaEqn().relax();
+    nuTildaEqn().solve();
+
     bound(nuTilda_, dimensionedScalar("zero", nuTilda_.dimensions(), 0.0));
     nuTilda_.correctBoundaryConditions();
 
diff --git a/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C b/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C
index 68233e0f28bdf4c035de4ac39f1027d57a0070d1..08587e11e958e48966a41e27e91c6b50d9204537 100644
--- a/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C
+++ b/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C
@@ -118,7 +118,7 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
     volScalarField divU = fvc::div(phi()/fvc::interpolate(rho()));
     volScalarField G = 2*muSgs_*(gradU && D);
 
-    solve
+    tmp<fvScalarMatrix> kEqn
     (
         fvm::ddt(rho(), k_)
       + fvm::div(phi(), k_)
@@ -129,9 +129,8 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
       - fvm::Sp(ce_(D)*rho()*sqrt(k_)/delta(), k_)
     );
 
-    //FIXME: why not this?
-    // kEqn.relax();
-    // kEqn.solve();
+    kEqn().relax();
+    kEqn().solve();
 
     bound(k_, kMin_);
 
diff --git a/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C b/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C
index 696e9b8ddf738599ada4e90efeaafd046910acbc..4ab154f209d3b3c31ec1eeda67ceb6f096051d2e 100644
--- a/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C
+++ b/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C
@@ -106,7 +106,7 @@ void lowReOneEqEddy::correct(const tmp<volTensorField>& tgradU)
     volScalarField divU = fvc::div(phi()/fvc::interpolate(rho()));
     volScalarField G = 2*muSgs_*(gradU && dev(symm(gradU)));
 
-    solve
+    tmp<fvScalarMatrix> kEqn
     (
         fvm::ddt(rho(), k_)
       + fvm::div(phi(), k_)
@@ -117,6 +117,9 @@ void lowReOneEqEddy::correct(const tmp<volTensorField>& tgradU)
       - fvm::Sp(ce_*rho()*sqrt(k_)/delta(), k_)
     );
 
+    kEqn().relax();
+    kEqn().solve();
+
     bound(k_, kMin_);
 
     updateSubGridScaleFields();
diff --git a/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C b/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C
index dd032dad791a6dcfdca51ef08db37abf16776d73..b2a01a7ac59f8ce6e8cef57a089932953201ef75 100644
--- a/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C
+++ b/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C
@@ -93,7 +93,7 @@ void oneEqEddy::correct(const tmp<volTensorField>& tgradU)
     volScalarField divU = fvc::div(phi()/fvc::interpolate(rho()));
     volScalarField G = 2*muSgs_*(gradU && dev(symm(gradU)));
 
-    fvScalarMatrix kEqn
+    tmp<fvScalarMatrix> kEqn
     (
         fvm::ddt(rho(), k_)
       + fvm::div(phi(), k_)
@@ -104,8 +104,8 @@ void oneEqEddy::correct(const tmp<volTensorField>& tgradU)
       - fvm::Sp(ce_*rho()*sqrt(k_)/delta(), k_)
     );
 
-    kEqn.relax();
-    kEqn.solve();
+    kEqn().relax();
+    kEqn().solve();
 
     bound(k_, kMin_);
 
diff --git a/src/turbulenceModels/compressible/RAS/LRR/LRR.C b/src/turbulenceModels/compressible/RAS/LRR/LRR.C
index 06e23ecc9214a0131596a49870bba7838b9fd897..c567a4df2f54dca589ccdf057783848a0e84c32d 100644
--- a/src/turbulenceModels/compressible/RAS/LRR/LRR.C
+++ b/src/turbulenceModels/compressible/RAS/LRR/LRR.C
@@ -229,6 +229,7 @@ LRR::LRR
             << exit(FatalError);
     }
 
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
diff --git a/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
index e612cebecd5cf6bb521aecc0f55043f7e65518e8..8f019890db89736f62bd10e516d8ac0c7a4d8924 100644
--- a/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
+++ b/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
@@ -258,6 +258,7 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             << exit(FatalError);
     }
 
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
diff --git a/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
index 149ed0b28939e179b0583ca03e8e0dec0770bd2e..93f30ec930d644d6cf269d33444a520fc728958b 100644
--- a/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
@@ -187,6 +187,7 @@ LaunderSharmaKE::LaunderSharmaKE
         autoCreateAlphat("alphat", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     mut_ = rho_*Cmu_*fMu()*sqr(k_)/epsilon_;
diff --git a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C
index 5ac546d69ce50e3d9a13b1a80540a63f1421018f..2c35c3cec0fd88d0f744f379d43aa4003b439e45 100644
--- a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C
+++ b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C
@@ -81,7 +81,7 @@ RASModel::RASModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
+    kMin_("kMin", sqr(dimVelocity), SMALL),
     epsilonMin_("epsilonMin", kMin_.dimensions()/dimTime, SMALL),
     omegaMin_("omegaMin", dimless/dimTime, SMALL),
 
@@ -133,9 +133,13 @@ autoPtr<RASModel> RASModel::New
     {
         FatalErrorIn
         (
-            "RASModel::New(const volScalarField&, "
-            "const volVectorField&, const surfaceScalarField&, "
-            "basicThermo&)"
+            "RASModel::New"
+            "("
+                "const volScalarField&, "
+                "const volVectorField&, "
+                "const surfaceScalarField&, "
+                "basicThermo&"
+            ")"
         )   << "Unknown RASModel type " << modelName
             << endl << endl
             << "Valid RASModel types are :" << endl
diff --git a/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C
index fd3955a5730983a4bdb1caed80a107fc3686d8d1..3159265c96b9ba59eb761e94b4d706300c83ef46 100644
--- a/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C
@@ -186,6 +186,7 @@ RNGkEpsilon::RNGkEpsilon
         autoCreateAlphat("alphat", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
diff --git a/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C b/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
index 43e63367cefdadb69edc1b3231b0ee2b6f153380..973c802fd9e3da5080177e655e2c29d84a30aab6 100644
--- a/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
@@ -168,6 +168,7 @@ kEpsilon::kEpsilon
         autoCreateAlphat("alphat", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
diff --git a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
index 6f9da54c9746d37fb2b67c494c27aac843a844a7..123de2099dbad5eb0a903542ff30f6c4670ac9c8 100644
--- a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
@@ -258,6 +258,7 @@ kOmegaSST::kOmegaSST
         autoCreateAlphat("alphat", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(omega_, omegaMin_);
 
     mut_ =
diff --git a/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C b/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
index a5a0c562a13d69cd015a388f161dd4757ce3b8fe..06b7beedfbdeacf48741184fcb7799a86bb8397d 100644
--- a/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
+++ b/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
@@ -102,7 +102,7 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
     volScalarField K = 0.5*tr(B_);
     volScalarField Epsilon = 2*nuEff()*magSqr(D);
 
-    fvSymmTensorMatrix BEqn
+    tmp<fvSymmTensorMatrix> BEqn
     (
         fvm::ddt(B_)
       + fvm::div(phi(), B_)
@@ -114,8 +114,8 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
       - (2*ce_ - 0.667*cm_)*I*Epsilon
     );
 
-    BEqn.relax();
-    BEqn.solve();
+    BEqn().relax();
+    BEqn().solve();
 
     // Bounding the component kinetic energies
 
diff --git a/src/turbulenceModels/incompressible/LES/LESModel/LESModel.C b/src/turbulenceModels/incompressible/LES/LESModel/LESModel.C
index d32fd1647ac2549706d66352a8562fed7a6152bf..34966eb22c64375dafbcc0864e4553a25d2de8e6 100644
--- a/src/turbulenceModels/incompressible/LES/LESModel/LESModel.C
+++ b/src/turbulenceModels/incompressible/LES/LESModel/LESModel.C
@@ -78,7 +78,7 @@ LESModel::LESModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
+    kMin_("kMin", sqr(dimVelocity), SMALL),
     delta_(LESdelta::New("delta", U.mesh(), *this))
 {
     readIfPresent("kMin", kMin_);
@@ -128,8 +128,12 @@ autoPtr<LESModel> LESModel::New
     {
         FatalErrorIn
         (
-            "LESModel::New(const volVectorField& U, const "
-            "surfaceScalarField& phi, transportModel&)"
+            "LESModel::New"
+            "("
+                "const volVectorField&, "
+                "const surfaceScalarField& ,"
+                "transportModel&"
+            ")"
         )   << "Unknown LESModel type " << modelName
             << endl << endl
             << "Valid LESModel types are :" << endl
@@ -167,7 +171,7 @@ bool LESModel::read()
 
         delta_().read(*this);
 
-        readIfPresent("kMin", kMin_);
+        kMin_.readIfPresent(*this);
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C b/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C
index 10c0116b6a3f332d441cc1ae0ebd5fa07cf42186..5e31e281dc2f3e96fc9313f73d7acc2db6c19450 100644
--- a/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C
+++ b/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C
@@ -90,7 +90,10 @@ LRRDiffStress::LRRDiffStress
         )
     )
 {
-    updateSubGridScaleFields(0.5*tr(B_));
+    volScalarField K = 0.5*tr(B_);
+    bound(K, kMin_);
+
+    updateSubGridScaleFields(K);
 
     printCoeffs();
 }
@@ -111,7 +114,7 @@ void LRRDiffStress::correct(const tmp<volTensorField>& tgradU)
     volScalarField K = 0.5*tr(B_);
     volScalarField Epsilon = 2*nuEff()*magSqr(D);
 
-    fvSymmTensorMatrix BEqn
+    tmp<fvSymmTensorMatrix> BEqn
     (
         fvm::ddt(B_)
       + fvm::div(phi(), B_)
@@ -124,8 +127,8 @@ void LRRDiffStress::correct(const tmp<volTensorField>& tgradU)
       - (0.667 - 2*c1_)*I*pow(K, 1.5)/delta()
     );
 
-    BEqn.relax();
-    BEqn.solve();
+    BEqn().relax();
+    BEqn().solve();
 
     // Bounding the component kinetic energies
 
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
index 0932bb524eeed23598afcb244f8fe0639390b062..431656cc6fe03e4437135052e58cfec85584344a 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
@@ -295,7 +295,7 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
     const volScalarField dTilda = this->dTilda(S);
     const volScalarField STilda = this->STilda(S, dTilda);
 
-    fvScalarMatrix nuTildaEqn
+    tmp<fvScalarMatrix> nuTildaEqn
     (
         fvm::ddt(nuTilda_)
       + fvm::div(phi(), nuTilda_)
@@ -311,8 +311,8 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
       - fvm::Sp(Cw1_*fw(STilda, dTilda)*nuTilda_/sqr(dTilda), nuTilda_)
     );
 
-    nuTildaEqn.relax();
-    nuTildaEqn.solve();
+    nuTildaEqn().relax();
+    nuTildaEqn().solve();
 
     bound(nuTilda_, dimensionedScalar("zero", nuTilda_.dimensions(), 0.0));
     nuTilda_.correctBoundaryConditions();
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C
index 22f9b638e21dadc30e3a0f027fa61d073290304f..48dac551fa5fd40e2322db493db95befb4b7fb8c 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C
@@ -176,7 +176,6 @@ SpalartAllmarasIDDES::SpalartAllmarasIDDES
             1.63
         )
     )
-
 {}
 
 
diff --git a/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C b/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
index 08f9c62fe7be657ac2f44d402d2c1e3c71c2edbc..ab777728354d6f5a868a12084417159b5ecfa81b 100644
--- a/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
+++ b/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
@@ -127,6 +127,8 @@ dynOneEqEddy::dynOneEqEddy
     filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
     filter_(filterPtr_())
 {
+    bound(k_, kMin_);
+
     updateSubGridScaleFields(symm(fvc::grad(U)));
 
     printCoeffs();
@@ -145,7 +147,7 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
 
     volScalarField P = 2.0*nuSgs_*magSqr(D);
 
-    fvScalarMatrix kEqn
+    tmp<fvScalarMatrix> kEqn
     (
        fvm::ddt(k_)
      + fvm::div(phi(), k_)
@@ -155,8 +157,8 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
      - fvm::Sp(ce(D)*sqrt(k_)/delta(), k_)
     );
 
-    kEqn.relax();
-    kEqn.solve();
+    kEqn().relax();
+    kEqn().solve();
 
     bound(k_, kMin_);
 
diff --git a/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.C b/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.C
index fedd18ee7d169d5fcb7f2741dc1e05f82f78db2a..f2bf70721feceff64543e5e97c477e77f7d8e647 100644
--- a/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.C
+++ b/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.C
@@ -118,6 +118,8 @@ dynSmagorinsky::dynSmagorinsky
     filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
     filter_(filterPtr_())
 {
+    bound(k_,  kMin_);
+
     updateSubGridScaleFields(dev(symm(fvc::grad(U))));
 
     printCoeffs();
@@ -133,6 +135,7 @@ void dynSmagorinsky::correct(const tmp<volTensorField>& gradU)
     volSymmTensorField D = dev(symm(gradU));
 
     k_ = cI(D)*sqr(delta())*magSqr(D);
+    bound(k_,  kMin_);
 
     updateSubGridScaleFields(D);
 }
diff --git a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
index 86cd8227d10ff732b21c801ea95aa118cd548111..708c090351d36c772f185ec0068e10dedc2f91be 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
@@ -323,6 +323,7 @@ kOmegaSSTSAS::kOmegaSSTSAS
         mesh_
     )
 {
+    bound(k_, kMin_);
     bound(omega_, omegaMin_);
 
     updateSubGridScaleFields(magSqr(symm(fvc::grad(U))));
diff --git a/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C b/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C
index 3decb7080e9be7fbdb7e1f13e64bcb82e2d67ee7..93db2120edd700866b2340144d20ab0ac6daac1f 100644
--- a/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C
+++ b/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C
@@ -119,6 +119,8 @@ locDynOneEqEddy::locDynOneEqEddy
     filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
     filter_(filterPtr_())
 {
+    bound(k_, kMin_);
+
     volScalarField KK = 0.5*(filter_(magSqr(U)) - magSqr(filter_(U)));
     updateSubGridScaleFields(symm(fvc::grad(U)), KK);
 
@@ -139,7 +141,7 @@ void locDynOneEqEddy::correct(const tmp<volTensorField>& gradU)
 
     volScalarField P = 2.0*nuSgs_*magSqr(D);
 
-    fvScalarMatrix kEqn
+    tmp<fvScalarMatrix> kEqn
     (
        fvm::ddt(k_)
      + fvm::div(phi(), k_)
@@ -149,8 +151,8 @@ void locDynOneEqEddy::correct(const tmp<volTensorField>& gradU)
      - fvm::Sp(ce(D, KK)*sqrt(k_)/delta(), k_)
     );
 
-    kEqn.relax();
-    kEqn.solve();
+    kEqn().relax();
+    kEqn().solve();
 
     bound(k_, kMin_);
 
diff --git a/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C b/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C
index 96188e59e4cbc1747e75e7c9c6e2a5ecba8e369d..89ad1892e35ad07e335b44fbfa55b74c5dcff2b6 100644
--- a/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C
+++ b/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C
@@ -86,6 +86,8 @@ oneEqEddy::oneEqEddy
         )
     )
 {
+    bound(k_, kMin_);
+
     updateSubGridScaleFields();
 
     printCoeffs();
@@ -100,7 +102,7 @@ void oneEqEddy::correct(const tmp<volTensorField>& gradU)
 
     volScalarField G = 2.0*nuSgs_*magSqr(symm(gradU));
 
-    fvScalarMatrix kEqn
+    tmp<fvScalarMatrix> kEqn
     (
        fvm::ddt(k_)
      + fvm::div(phi(), k_)
@@ -110,8 +112,8 @@ void oneEqEddy::correct(const tmp<volTensorField>& gradU)
      - fvm::Sp(ce_*sqrt(k_)/delta(), k_)
     );
 
-    kEqn.relax();
-    kEqn.solve();
+    kEqn().relax();
+    kEqn().solve();
 
     bound(k_, kMin_);
 
diff --git a/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.C b/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.C
index 4a16af3576c37166c229ea34e97bc78e2805cda9..f6ab218082d99d3c6be18e74c82d0f53c897a946 100644
--- a/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.C
+++ b/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.C
@@ -72,7 +72,6 @@ spectEddyVisc::spectEddyVisc
     LESModel(typeName, U, phi, transport),
     GenEddyVisc(U, phi, transport),
 
-
     cB_
     (
         dimensioned<scalar>::lookupOrAddToDict
diff --git a/src/turbulenceModels/incompressible/RAS/LRR/LRR.C b/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
index 037da5b430aff25fce6f6d0c6eac469b8164bf17..ae79bb4389bfa00896e6b28f9885f4d058a02c43 100644
--- a/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
+++ b/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
@@ -198,6 +198,7 @@ LRR::LRR
             << exit(FatalError);
     }
 
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*sqr(k_)/epsilon_;
diff --git a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
index b6d20aea2d41f29f303e81da6c3df95af3be8adf..5f036cff106db453c031fa064bbee469c5123988 100644
--- a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
@@ -140,6 +140,7 @@ LamBremhorstKE::LamBremhorstKE
         autoCreateLowReNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*fMu_*sqr(k_)/epsilon_;
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
index 98001c6a68811bc8a753ec6079cfab0dc8505391..41f696222463a969dbdcc6720015557d2594f0fa 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
@@ -224,6 +224,7 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*sqr(k_)/epsilon_;
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
index 6f840c57d38d1d3a880ef854a40b04e0ae7af44b..e62cbf2578f52ce83bd29ca3f4e5c4ddebb05b79 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
@@ -146,6 +146,7 @@ LaunderSharmaKE::LaunderSharmaKE
         autoCreateLowReNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilonTilda_, epsilonMin_);
 
     nut_ = Cmu_*fMu()*sqr(k_)/epsilonTilda_;
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
index f01805f115606666ca4da22ca493277199a5e69b..d91f974fbe8f2a7f210747d659ee355e2d93d3c5 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
@@ -170,11 +170,17 @@ LienCubicKE::LienCubicKE
         autoCreateEpsilon("epsilon", mesh_)
     ),
 
-    //FIXME - epsilon is not bounded
-
     gradU_(fvc::grad(U)),
-    eta_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))),
-    ksi_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))),
+    eta_
+    (
+        k_/(epsilon_ + epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))
+    ),
+    ksi_
+    (
+        k_/(epsilon_ + epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))
+    ),
     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
     fEta_(A2_ + pow(eta_, 3.0)),
 
@@ -228,6 +234,7 @@ LienCubicKE::LienCubicKE
         )
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*sqr(k_)/epsilon_ + C5viscosity_;
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
index 8c7bbc7c7b09b0848b1ca03bd1f665a90757c60d..046f5a6fa1b339ea2968223560624faadaffba22 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
@@ -220,14 +220,22 @@ LienCubicKELowRe::LienCubicKELowRe
     y_(mesh_),
 
     gradU_(fvc::grad(U)),
-    eta_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))),
-    ksi_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))),
+    eta_
+    (
+        k_/(epsilon_ + epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))
+    ),
+    ksi_
+    (
+        k_/(epsilon_ + epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))
+    ),
     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
-    fEta_(A2_ + pow(eta_, 3.0)),
+    fEta_(A2_ + pow3(eta_)),
 
     C5viscosity_
     (
-        -2.0*pow(Cmu_, 3.0)*pow(k_, 4.0)/pow(epsilon_, 3.0)
+        -2.0*pow3(Cmu_)*pow4(k_)/pow3(epsilon_ + epsilonMin_)
        *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()))
     ),
 
@@ -252,7 +260,7 @@ LienCubicKELowRe::LienCubicKELowRe
         symm
         (
             // quadratic terms
-            pow(k_, 3.0)/sqr(epsilon_)
+            pow3(k_)/sqr(epsilon_ + epsilonMin_)
            *(
                 Ctau1_/fEta_
                *(
@@ -263,8 +271,8 @@ LienCubicKELowRe::LienCubicKELowRe
               + Ctau3_/fEta_*(gradU_.T() & gradU_)
             )
             // cubic term C4
-          - 20.0*pow(k_, 4.0)/pow(epsilon_, 3.0)
-           *pow(Cmu_, 3.0)
+          - 20.0*pow4(k_)/pow3(epsilon_ + epsilonMin_)
+           *pow3(Cmu_)
            *(
                 ((gradU_ & gradU_) & gradU_.T())
               + ((gradU_ & gradU_.T()) & gradU_.T())
@@ -280,6 +288,7 @@ LienCubicKELowRe::LienCubicKELowRe
         )
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_
@@ -430,8 +439,8 @@ void LienCubicKELowRe::correct()
 
     epsEqn().relax();
 
-#   include "LienCubicKELowReSetWallDissipation.H"
-#   include "wallDissipationI.H"
+    #include "LienCubicKELowReSetWallDissipation.H"
+    #include "wallDissipationI.H"
 
     solve(epsEqn);
     bound(epsilon_, epsilonMin_);
diff --git a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
index a2a807676dd8bca6f65b5bd04acff34993d56242..72a455d945fb04a3f509e193315683e0d5b9ab76 100644
--- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
@@ -180,6 +180,7 @@ LienLeschzinerLowRe::LienLeschzinerLowRe
         autoCreateLowReNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*(scalar(1) - exp(-Am_*yStar_))
diff --git a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
index 01324af90a21f5f9a0d2e513e5018560ae071fe6..c885d87f163ec5cce7ad3e4a790aa2ee4d8faec5 100644
--- a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
+++ b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
@@ -144,7 +144,6 @@ NonlinearKEShih::NonlinearKEShih
         )
     ),
 
-    //FIXME: should be named 'kappa_' or 'kappa'?
     kappa_
     (
         dimensioned<scalar>::lookupOrAddToDict
@@ -175,7 +174,6 @@ NonlinearKEShih::NonlinearKEShih
             IOobject::AUTO_WRITE
         ),
         mesh_
-        //FIXME: what about autoCreateK("k", mesh_)
     ),
 
     epsilon_
@@ -189,39 +187,30 @@ NonlinearKEShih::NonlinearKEShih
             IOobject::AUTO_WRITE
         ),
         mesh_
-        //FIXME: what about autoCreateK("epsilon", mesh_)
     ),
 
-    //FIXME: epsilon is not bounded
-
     gradU_(fvc::grad(U)),
-    eta_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))),
-    ksi_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))),
+    eta_
+    (
+        k_/(epsilon_ + epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))
+    ),
+    ksi_
+    (
+        k_/(epsilon_+ epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))
+    ),
     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
     fEta_(A2_ + pow(eta_, 3.0)),
 
-    // FIXME: epsilon is not bounded
     nut_("nut", Cmu_*sqr(k_)/(epsilon_ + epsilonMin_)),
-    // FIXME: why not use the following?
-    // nut_
-    // (
-    //     IOobject
-    //     (
-    //         "nut",
-    //         runTime_.timeName(),
-    //         mesh_,
-    //         IOobject::NO_READ,
-    //         IOobject::AUTO_WRITE
-    //     ),
-    //     autoCreateNut("nut", mesh_)
-    // ),
 
     nonlinearStress_
     (
         "nonlinearStress",
         symm
         (
-            pow(k_, 3.0)/sqr(epsilon_)
+            pow3(k_)/sqr(epsilon_ + epsilonMin_)
            *(
                 Ctau1_/fEta_
                *(
@@ -234,11 +223,9 @@ NonlinearKEShih::NonlinearKEShih
         )
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
-    //FIXME: could use this
-    // nut_ = Cmu_*sqr(k_)/epsilon_;
-
     #include "wallNonlinearViscosityI.H"
 
     printCoeffs();
diff --git a/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C b/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C
index f45c1ee0ff13fe3cf4c54cce36c75393336b83fc..c7d97e08e9f340da09d9541f373e14338e1331ce 100644
--- a/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C
+++ b/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C
@@ -80,7 +80,7 @@ RASModel::RASModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
+    kMin_("kMin", sqr(dimVelocity), SMALL),
     epsilonMin_("epsilonMin", kMin_.dimensions()/dimTime, SMALL),
     omegaMin_("omegaMin", dimless/dimTime, SMALL),
 
@@ -131,8 +131,12 @@ autoPtr<RASModel> RASModel::New
     {
         FatalErrorIn
         (
-            "RASModel::New(const volVectorField&, "
-            "const surfaceScalarField&, transportModel&)"
+            "RASModel::New"
+            "("
+                "const volVectorField&, "
+                "const surfaceScalarField&, "
+                "transportModel&"
+            ")"
         )   << "Unknown RASModel type " << modelName
             << endl << endl
             << "Valid RASModel types are :" << endl
diff --git a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
index 1591608717616dea676ac56fd4340dd71ebf20ae..6c8d874be3850fbc13f4db48f9f9b6dfff810b02 100644
--- a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
@@ -155,6 +155,7 @@ RNGkEpsilon::RNGkEpsilon
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*sqr(k_)/epsilon_;
diff --git a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
index 8c41bf4d2786f5f2bc5da29e466611ccadec0187..53f2e041494f0233d6994b2880a04519fc6e9354 100644
--- a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
@@ -128,6 +128,7 @@ kEpsilon::kEpsilon
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*sqr(k_)/epsilon_;
diff --git a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
index 04f505376b9e36715c4e0dde6d95ca11a93ceaa0..77225d5edeb55b3174f3c7ed9d4e881d9279e8f5 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
@@ -137,6 +137,7 @@ kOmega::kOmega
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(omega_, omegaMin_);
 
     nut_ = k_/omega_;
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
index 4332fc3706fac8a15d6db7d9f6d62152030dd993..4d37b287528d9cfcf70db6a173984864a454c411 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
@@ -236,6 +236,7 @@ kOmegaSST::kOmegaSST
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(omega_, omegaMin_);
 
     nut_ =
diff --git a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
index d72fd30f3e7a6d439971c6fb64d53999919b7816..422f8b1463a40c8dd680015d34d3bdf91525b935 100644
--- a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
+++ b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
@@ -196,6 +196,7 @@ qZeta::qZeta
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
     bound(q_, qMin_);
     bound(zeta_, zetaMin_);