diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
index fdddddfbfe..821eb5adb0 100644
--- a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -121,9 +121,8 @@ PDRkEpsilon::PDRkEpsilon
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateK("k", mesh_)
     ),
-
     epsilon_
     (
         IOobject
@@ -134,9 +133,8 @@ PDRkEpsilon::PDRkEpsilon
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateEpsilon("epsilon", mesh_)
     ),
-
     mut_
     (
         IOobject
@@ -147,9 +145,8 @@ PDRkEpsilon::PDRkEpsilon
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)
+        autoCreateMut("mut", mesh_)
     ),
-
     alphat_
     (
         IOobject
@@ -163,7 +160,9 @@ PDRkEpsilon::PDRkEpsilon
         autoCreateAlphat("alphat", mesh_)
     )
 {
-    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    bound(epsilon_, epsilonMin_);
+
+    mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 
     alphat_ = mut_/Prt_;
@@ -250,7 +249,7 @@ void PDRkEpsilon::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
         mut_.correctBoundaryConditions();
 
         // Re-calculate thermal diffusivity
@@ -302,7 +301,7 @@ void PDRkEpsilon::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -320,7 +319,7 @@ void PDRkEpsilon::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
     // Re-calculate viscosity
     mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
diff --git a/applications/solvers/multiphase/settlingFoam/kEpsilon.H b/applications/solvers/multiphase/settlingFoam/kEpsilon.H
index 522a9afac7..e27594ffed 100644
--- a/applications/solvers/multiphase/settlingFoam/kEpsilon.H
+++ b/applications/solvers/multiphase/settlingFoam/kEpsilon.H
@@ -5,8 +5,8 @@ if (turbulence)
         y.correct();
     }
 
-    dimensionedScalar k0("k0", k.dimensions(), SMALL);
-    dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), SMALL);
+    dimensionedScalar kMin("kMin", k.dimensions(), SMALL);
+    dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL);
 
     volScalarField divU = fvc::div(phi/fvc::interpolate(rho));
 
@@ -15,7 +15,7 @@ if (turbulence)
     tgradU.clear();
 
     volScalarField Gcoef =
-        Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilon0);
+        Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin);
 
     #include "wallFunctions.H"
 
@@ -40,7 +40,7 @@ if (turbulence)
     epsEqn.relax();
     epsEqn.solve();
 
-    bound(epsilon, epsilon0);
+    bound(epsilon, epsilonMin);
 
 
     // Turbulent kinetic energy equation
@@ -60,11 +60,13 @@ if (turbulence)
       - fvm::Sp(rho*epsilon/k, k)
     );
 
-    bound(k, k0);
+    //FIXME: why no kEqn.relax()?
+
+    bound(k, kMin);
 
 
     //- Re-calculate viscosity
-    mut = rho*Cmu*sqr(k)/(epsilon + epsilon0);
+    mut = rho*Cmu*sqr(k)/(epsilon + epsilonMin);
 
     #include "wallViscosity.H"
 }
diff --git a/src/finiteVolume/cfdTools/general/bound/bound.C b/src/finiteVolume/cfdTools/general/bound/bound.C
index b629033754..7054434c0f 100644
--- a/src/finiteVolume/cfdTools/general/bound/bound.C
+++ b/src/finiteVolume/cfdTools/general/bound/bound.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,11 +30,11 @@ License
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
-void Foam::bound(volScalarField& vsf, const dimensionedScalar& vsf0)
+void Foam::bound(volScalarField& vsf, const dimensionedScalar& lowerBound)
 {
-    scalar minVsf = min(vsf).value();
+    const scalar minVsf = min(vsf).value();
 
-    if (minVsf < vsf0.value())
+    if (minVsf < lowerBound.value())
     {
         Info<< "bounding " << vsf.name()
             << ", min: " << minVsf
@@ -47,13 +47,13 @@ void Foam::bound(volScalarField& vsf, const dimensionedScalar& vsf0)
             max
             (
                 vsf.internalField(),
-                fvc::average(max(vsf, vsf0))().internalField()
-                *pos(-vsf.internalField())
+                fvc::average(max(vsf, lowerBound))().internalField()
+              * pos(-vsf.internalField())
             ),
-            vsf0.value()
+            lowerBound.value()
         );
 
-        vsf.boundaryField() = max(vsf.boundaryField(), vsf0.value());
+        vsf.boundaryField() = max(vsf.boundaryField(), lowerBound.value());
     }
 }
 
diff --git a/src/finiteVolume/cfdTools/general/bound/bound.H b/src/finiteVolume/cfdTools/general/bound/bound.H
index 08a58e5550..f301c5a217 100644
--- a/src/finiteVolume/cfdTools/general/bound/bound.H
+++ b/src/finiteVolume/cfdTools/general/bound/bound.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,11 +23,13 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 InNamespace
-    Foam::bound
+    Foam
 
 Description
-    Bound the given scalar field if it has gone unbounded.  Used extensively
-    in RAS and LES turbulence models.
+    Bound the given scalar field if it has gone unbounded.
+
+    Used extensively in RAS and LES turbulence models, but also of use
+    within solvers.
 
 SourceFiles
     bound.C
@@ -49,7 +51,7 @@ namespace Foam
 
 //- Bound the given scalar field if it has gone unbounded.
 //  Used extensively in RAS and LES turbulence models.
-void bound(volScalarField& vsf, const dimensionedScalar& vsf0);
+void bound(volScalarField&, const dimensionedScalar& lowerBound);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C b/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
index 41063cd2cd..cff839b5e3 100644
--- a/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
+++ b/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -132,7 +132,7 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
     }
 
     K = 0.5*tr(B_);
-    bound(K, k0());
+    bound(K, kMin_);
 
     updateSubGridScaleFields(K);
 }
diff --git a/src/turbulenceModels/compressible/LES/LESModel/LESModel.C b/src/turbulenceModels/compressible/LES/LESModel/LESModel.C
index 6ac295dea6..daeec0a05e 100644
--- a/src/turbulenceModels/compressible/LES/LESModel/LESModel.C
+++ b/src/turbulenceModels/compressible/LES/LESModel/LESModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,11 +79,11 @@ LESModel::LESModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    k0_("k0", dimVelocity*dimVelocity, SMALL),
+    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
 
     delta_(LESdelta::New("delta", U.mesh(), *this))
 {
-    readIfPresent("k0", k0_);
+    readIfPresent("kMin", kMin_);
 
     // Force the construction of the mesh deltaCoeffs which may be needed
     // for the construction of the derived models and BCs
@@ -167,7 +167,7 @@ bool LESModel::read()
             coeffDict_ <<= *dictPtr;
         }
 
-        readIfPresent("k0", k0_);
+        readIfPresent("kMin", kMin_);
 
         delta_().read(*this);
 
diff --git a/src/turbulenceModels/compressible/LES/LESModel/LESModel.H b/src/turbulenceModels/compressible/LES/LESModel/LESModel.H
index 2cdde6c61e..6b10b930cc 100644
--- a/src/turbulenceModels/compressible/LES/LESModel/LESModel.H
+++ b/src/turbulenceModels/compressible/LES/LESModel/LESModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -83,7 +83,7 @@ protected:
         Switch printCoeffs_;
         dictionary coeffDict_;
 
-        dimensionedScalar k0_;
+        dimensionedScalar kMin_;
 
         autoPtr<LESdelta> delta_;
 
@@ -170,16 +170,16 @@ public:
                 return coeffDict_;
             }
 
-            //- Return the value of k0 which k is not allowed to be less than
-            const dimensionedScalar& k0() const
+            //- Return the lower allowable limit for k (default: SMALL)
+            const dimensionedScalar& kMin() const
             {
-                return k0_;
+                return kMin_;
             }
 
-            //- Allow k0 to be changed
-            dimensionedScalar& k0()
+            //- Allow kMin to be changed
+            dimensionedScalar& kMin()
             {
-                return k0_;
+                return kMin_;
             }
 
             //- Access function to filter width
diff --git a/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C b/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C
index edd0f4bcf0..68233e0f28 100644
--- a/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C
+++ b/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,10 +74,10 @@ dimensionedScalar dynOneEqEddy::ce_(const volSymmTensorField& D) const
         pow(KK + filter_(k_), 1.5)/(2*delta()) - filter_(pow(k_, 1.5))/delta();
 
     volScalarField ee =
-        2*delta()*ck_(D)*
-        (
+        2*delta()*ck_(D)
+       *(
             filter_(sqrt(k_)*magSqr(D))
-            - 2*sqrt(KK + filter_(k_))*magSqr(filter_(D))
+          - 2*sqrt(KK + filter_(k_))*magSqr(filter_(D))
         );
 
     return average(ee*mm)/average(mm*mm);
@@ -129,7 +129,11 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
       - fvm::Sp(ce_(D)*rho()*sqrt(k_)/delta(), k_)
     );
 
-    bound(k_, dimensionedScalar("0", k_.dimensions(), 1.0e-10));
+    //FIXME: why not this?
+    // kEqn.relax();
+    // kEqn.solve();
+
+    bound(k_, kMin_);
 
     updateSubGridScaleFields(D);
 }
diff --git a/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C b/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C
index 77b12af60e..696e9b8ddf 100644
--- a/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C
+++ b/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -117,7 +117,7 @@ void lowReOneEqEddy::correct(const tmp<volTensorField>& tgradU)
       - fvm::Sp(ce_*rho()*sqrt(k_)/delta(), k_)
     );
 
-    bound(k_, k0());
+    bound(k_, kMin_);
 
     updateSubGridScaleFields();
 }
diff --git a/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C b/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C
index c6212dd74c..dd032dad79 100644
--- a/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C
+++ b/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,7 +107,7 @@ void oneEqEddy::correct(const tmp<volTensorField>& tgradU)
     kEqn.relax();
     kEqn.solve();
 
-    bound(k_, k0());
+    bound(k_, kMin_);
 
     updateSubGridScaleFields();
 }
diff --git a/src/turbulenceModels/compressible/RAS/LRR/LRR.C b/src/turbulenceModels/compressible/RAS/LRR/LRR.C
index e2d925c3bb..06e23ecc92 100644
--- a/src/turbulenceModels/compressible/RAS/LRR/LRR.C
+++ b/src/turbulenceModels/compressible/RAS/LRR/LRR.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -229,7 +229,9 @@ LRR::LRR
             << exit(FatalError);
     }
 
-    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    bound(epsilon_, epsilonMin_);
+
+    mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 
     alphat_ = mut_/Prt_;
@@ -324,7 +326,7 @@ void LRR::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
         mut_.correctBoundaryConditions();
 
         // Re-calculate thermal diffusivity
@@ -359,7 +361,7 @@ void LRR::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Reynolds stress equation
@@ -406,15 +408,15 @@ void LRR::correct()
             R_.dimensions(),
             symmTensor
             (
-                k0_.value(), -GREAT, -GREAT,
-                k0_.value(), -GREAT,
-                k0_.value()
+                kMin_.value(), -GREAT, -GREAT,
+                kMin_.value(), -GREAT,
+                kMin_.value()
             )
         )
     );
 
     k_ = 0.5*tr(R_);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
index 66c54ec222..e612cebecd 100644
--- a/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
+++ b/src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -258,7 +258,9 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             << exit(FatalError);
     }
 
-    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    bound(epsilon_, epsilonMin_);
+
+    mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 
     alphat_ = mut_/Prt_;
@@ -357,7 +359,7 @@ void LaunderGibsonRSTM::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
         mut_.correctBoundaryConditions();
 
         // Re-calculate thermal diffusivity
@@ -397,7 +399,7 @@ void LaunderGibsonRSTM::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Reynolds stress equation
@@ -453,15 +455,15 @@ void LaunderGibsonRSTM::correct()
             R_.dimensions(),
             symmTensor
             (
-                k0_.value(), -GREAT, -GREAT,
-                k0_.value(), -GREAT,
-                k0_.value()
+                kMin_.value(), -GREAT, -GREAT,
+                kMin_.value(), -GREAT,
+                kMin_.value()
             )
         )
     );
 
     k_ == 0.5*tr(R_);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate turbulent viscosity
diff --git a/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
index f1dc4e9203..149ed0b289 100644
--- a/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -187,7 +187,9 @@ LaunderSharmaKE::LaunderSharmaKE
         autoCreateAlphat("alphat", mesh_)
     )
 {
-    mut_ = rho_*Cmu_*fMu()*sqr(k_)/(epsilon_ + epsilonSmall_);
+    bound(epsilon_, epsilonMin_);
+
+    mut_ = rho_*Cmu_*fMu()*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 
     alphat_ = mut_/Prt_;
@@ -275,7 +277,7 @@ void LaunderSharmaKE::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ == rho_*Cmu_*fMu()*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_ == rho_*Cmu_*fMu()*sqr(k_)/epsilon_;
 
         // Re-calculate thermal diffusivity
         alphat_ = mut_/Prt_;
@@ -320,7 +322,7 @@ void LaunderSharmaKE::correct()
 
     epsEqn().relax();
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -338,11 +340,11 @@ void LaunderSharmaKE::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
-    mut_ == Cmu_*fMu()*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_ == Cmu_*fMu()*rho_*sqr(k_)/epsilon_;
 
 
     // Re-calculate thermal diffusivity
diff --git a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C
index e5aa916702..5ac546d69c 100644
--- a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C
+++ b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.C
@@ -81,11 +81,9 @@ RASModel::RASModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    k0_("k0", dimVelocity*dimVelocity, SMALL),
-    epsilon0_("epsilon0", k0_.dimensions()/dimTime, SMALL),
-    epsilonSmall_("epsilonSmall", epsilon0_.dimensions(), SMALL),
-    omega0_("omega0", dimless/dimTime, SMALL),
-    omegaSmall_("omegaSmall", omega0_.dimensions(), SMALL),
+    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
+    epsilonMin_("epsilonMin", kMin_.dimensions()/dimTime, SMALL),
+    omegaMin_("omegaMin", dimless/dimTime, SMALL),
 
     y_(mesh_)
 {
@@ -219,11 +217,9 @@ bool RASModel::read()
             coeffDict_ <<= *dictPtr;
         }
 
-        k0_.readIfPresent(*this);
-        epsilon0_.readIfPresent(*this);
-        epsilonSmall_.readIfPresent(*this);
-        omega0_.readIfPresent(*this);
-        omegaSmall_.readIfPresent(*this);
+        kMin_.readIfPresent(*this);
+        epsilonMin_.readIfPresent(*this);
+        omegaMin_.readIfPresent(*this);
 
         return true;
     }
diff --git a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.H b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.H
index d5a12902cc..9bddec38aa 100644
--- a/src/turbulenceModels/compressible/RAS/RASModel/RASModel.H
+++ b/src/turbulenceModels/compressible/RAS/RASModel/RASModel.H
@@ -93,19 +93,13 @@ protected:
         scalar yPlusLam_;
 
         //- Lower limit of k
-        dimensionedScalar k0_;
+        dimensionedScalar kMin_;
 
         //- Lower limit of epsilon
-        dimensionedScalar epsilon0_;
-
-        //- Small epsilon value used to avoid divide by zero
-        dimensionedScalar epsilonSmall_;
+        dimensionedScalar epsilonMin_;
 
         //- Lower limit for omega
-        dimensionedScalar omega0_;
-
-        //- Small omega value used to avoid divide by zero
-        dimensionedScalar omegaSmall_;
+        dimensionedScalar omegaMin_;
 
         //- Near wall distance boundary field
         nearWallDist y_;
@@ -185,66 +179,40 @@ public:
 
         // Access
 
-            //- Return the lower allowable limit for k
-            const dimensionedScalar& k0() const
+            //- Return the lower allowable limit for k (default: SMALL)
+            const dimensionedScalar& kMin() const
             {
-                return k0_;
+                return kMin_;
             }
 
-            //- Return the lower allowable limit for epsilon
-            const dimensionedScalar& epsilon0() const
+            //- Return the lower allowable limit for epsilon (default: SMALL)
+            const dimensionedScalar& epsilonMin() const
             {
-                return epsilon0_;
+                return epsilonMin_;
             }
 
-            //- Return the value of epsilonSmall which is added to epsilon when
-            //  calculating mut
-            const dimensionedScalar& epsilonSmall() const
+            //- Return the lower allowable limit for omega (default: SMALL)
+            const dimensionedScalar& omegaMin() const
             {
-                return epsilonSmall_;
+                return omegaMin_;
             }
 
-            //- Return the lower allowable limit for omega
-            const dimensionedScalar& omega0() const
+            //- Allow kMin to be changed
+            dimensionedScalar& kMin()
             {
-                return omega0_;
+                return kMin_;
             }
 
-            //- Return the value of omegaSmall which is added to omega when
-            //  calculating mut
-            const dimensionedScalar& omegaSmall() const
+            //- Allow epsilonMin to be changed
+            dimensionedScalar& epsilonMin()
             {
-                return omegaSmall_;
+                return epsilonMin_;
             }
 
-            //- Allow k0 to be changed
-            dimensionedScalar& k0()
+            //- Allow omegaMin to be changed
+            dimensionedScalar& omegaMin()
             {
-                return k0_;
-            }
-
-            //- Allow epsilon0 to be changed
-            dimensionedScalar& epsilon0()
-            {
-                return epsilon0_;
-            }
-
-            //- Allow epsilonSmall to be changed
-            dimensionedScalar& epsilonSmall()
-            {
-                return epsilonSmall_;
-            }
-
-            //- Allow omega0 to be changed
-            dimensionedScalar& omega0()
-            {
-                return omega0_;
-            }
-
-            //- Allow omegaSmall to be changed
-            dimensionedScalar& omegaSmall()
-            {
-                return omegaSmall_;
+                return omegaMin_;
             }
 
             //- Return the near wall distances
diff --git a/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C
index 2fc2bdab29..fd3955a573 100644
--- a/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -186,7 +186,9 @@ RNGkEpsilon::RNGkEpsilon
         autoCreateAlphat("alphat", mesh_)
     )
 {
-    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    bound(epsilon_, epsilonMin_);
+
+    mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 
     alphat_ = mut_/Prt_;
@@ -276,7 +278,7 @@ void RNGkEpsilon::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
         mut_.correctBoundaryConditions();
 
         // Re-calculate thermal diffusivity
@@ -327,7 +329,7 @@ void RNGkEpsilon::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -339,12 +341,12 @@ void RNGkEpsilon::correct()
       - fvm::laplacian(DkEff(), k_)
       ==
         G - fvm::SuSp(2.0/3.0*rho_*divU, k_)
-      - fvm::Sp(rho_*(epsilon_)/k_, k_)
+      - fvm::Sp(rho_*epsilon_/k_, k_)
     );
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C b/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
index e0c04f1e44..43e63367ce 100644
--- a/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -168,7 +168,9 @@ kEpsilon::kEpsilon
         autoCreateAlphat("alphat", mesh_)
     )
 {
-    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    bound(epsilon_, epsilonMin_);
+
+    mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 
     alphat_ = mut_/Prt_;
@@ -257,7 +259,7 @@ void kEpsilon::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
         mut_.correctBoundaryConditions();
 
         // Re-calculate thermal diffusivity
@@ -300,7 +302,7 @@ void kEpsilon::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -318,7 +320,7 @@ void kEpsilon::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
index 23e9f0569f..6f9da54c97 100644
--- a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
@@ -258,12 +258,14 @@ kOmegaSST::kOmegaSST
         autoCreateAlphat("alphat", mesh_)
     )
 {
+    bound(omega_, omegaMin_);
+
     mut_ =
     (
         a1_*rho_*k_
       / max
         (
-            a1_*(omega_ + omegaSmall_),
+            a1_*omega_,
             F2()*sqrt(magSqr(symm(fvc::grad(U_))))
         )
     );
@@ -422,7 +424,7 @@ void kOmegaSST::correct()
     omegaEqn().boundaryManipulate(omega_.boundaryField());
 
     solve(omegaEqn);
-    bound(omega_, omega0_);
+    bound(omega_, omegaMin_);
 
     // Turbulent kinetic energy equation
     tmp<fvScalarMatrix> kEqn
@@ -438,7 +440,7 @@ void kOmegaSST::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/compressible/RAS/realizableKE/realizableKE.C b/src/turbulenceModels/compressible/RAS/realizableKE/realizableKE.C
index a79eff7b8a..c03e603acf 100644
--- a/src/turbulenceModels/compressible/RAS/realizableKE/realizableKE.C
+++ b/src/turbulenceModels/compressible/RAS/realizableKE/realizableKE.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,7 +69,7 @@ tmp<volScalarField> realizableKE::rCmu
     volScalarField As = sqrt(6.0)*cos(phis);
     volScalarField Us = sqrt(S2/2.0 + magSqr(skew(gradU)));
 
-    return 1.0/(A0_ + As*Us*k_/(epsilon_ + epsilonSmall_));
+    return 1.0/(A0_ + As*Us*k_/epsilon_);
 }
 
 
@@ -200,10 +200,10 @@ realizableKE::realizableKE
         autoCreateAlphat("alphat", mesh_)
     )
 {
-    bound(k_, k0_);
-    bound(epsilon_, epsilon0_);
+    bound(k_, kMin_);
+    bound(epsilon_, epsilonMin_);
 
-    mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 
     alphat_ = mut_/Prt_;
@@ -341,7 +341,7 @@ void realizableKE::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -353,12 +353,12 @@ void realizableKE::correct()
       - fvm::laplacian(DkEff(), k_)
      ==
         G - fvm::SuSp(2.0/3.0*rho_*divU, k_)
-      - fvm::Sp(rho_*(epsilon_)/k_, k_)
+      - fvm::Sp(rho_*epsilon_/k_, k_)
     );
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
     // Re-calculate viscosity
     mut_ = rCmu(gradU, S2, magS)*rho_*sqr(k_)/epsilon_;
diff --git a/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C b/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
index dd35b035ec..a5a0c562a1 100644
--- a/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
+++ b/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -122,15 +122,15 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
     forAll(B_, celli)
     {
         B_[celli].component(symmTensor::XX) =
-            max(B_[celli].component(symmTensor::XX), k0().value());
+            max(B_[celli].component(symmTensor::XX), kMin_.value());
         B_[celli].component(symmTensor::YY) =
-            max(B_[celli].component(symmTensor::YY), k0().value());
+            max(B_[celli].component(symmTensor::YY), kMin_.value());
         B_[celli].component(symmTensor::ZZ) =
-            max(B_[celli].component(symmTensor::ZZ), k0().value());
+            max(B_[celli].component(symmTensor::ZZ), kMin_.value());
     }
 
     K = 0.5*tr(B_);
-    bound(K, k0());
+    bound(K, kMin_);
 
     updateSubGridScaleFields(K);
 }
diff --git a/src/turbulenceModels/incompressible/LES/LESModel/LESModel.C b/src/turbulenceModels/incompressible/LES/LESModel/LESModel.C
index fd7791d9f3..d32fd1647a 100644
--- a/src/turbulenceModels/incompressible/LES/LESModel/LESModel.C
+++ b/src/turbulenceModels/incompressible/LES/LESModel/LESModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,10 +78,10 @@ LESModel::LESModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    k0_("k0", dimVelocity*dimVelocity, SMALL),
+    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
     delta_(LESdelta::New("delta", U.mesh(), *this))
 {
-    readIfPresent("k0", k0_);
+    readIfPresent("kMin", kMin_);
 
     // Force the construction of the mesh deltaCoeffs which may be needed
     // for the construction of the derived models and BCs
@@ -167,7 +167,7 @@ bool LESModel::read()
 
         delta_().read(*this);
 
-        readIfPresent("k0", k0_);
+        readIfPresent("kMin", kMin_);
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H b/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H
index 49fe52c454..8b4b43ceda 100644
--- a/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H
+++ b/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -83,7 +83,7 @@ protected:
         Switch printCoeffs_;
         dictionary coeffDict_;
 
-        dimensionedScalar k0_;
+        dimensionedScalar kMin_;
 
         autoPtr<LESdelta> delta_;
 
@@ -171,16 +171,16 @@ public:
             return delta_();
         }
 
-        //- Return the value of k0 which k is not allowed to be less than
-        const dimensionedScalar& k0() const
+        //- Return the lower allowable limit for k (default: SMALL)
+        const dimensionedScalar& kMin() const
         {
-            return k0_;
+            return kMin_;
         }
 
-        //- Allow k0 to be changed
-        dimensionedScalar& k0()
+        //- Allow kMin to be changed
+        dimensionedScalar& kMin()
         {
-            return k0_;
+            return kMin_;
         }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C b/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C
index 58368cc587..10c0116b6a 100644
--- a/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C
+++ b/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -132,15 +132,15 @@ void LRRDiffStress::correct(const tmp<volTensorField>& tgradU)
     forAll(B_, celli)
     {
         B_[celli].component(symmTensor::XX) =
-            max(B_[celli].component(symmTensor::XX), k0().value());
+            max(B_[celli].component(symmTensor::XX), kMin_.value());
         B_[celli].component(symmTensor::YY) =
-            max(B_[celli].component(symmTensor::YY), k0().value());
+            max(B_[celli].component(symmTensor::YY), kMin_.value());
         B_[celli].component(symmTensor::ZZ) =
-            max(B_[celli].component(symmTensor::ZZ), k0().value());
+            max(B_[celli].component(symmTensor::ZZ), kMin_.value());
     }
 
     K = 0.5*tr(B_);
-    bound(K, k0());
+    bound(K, kMin_);
 
     updateSubGridScaleFields(K);
 }
diff --git a/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C b/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
index 4aadf02bf2..08f9c62fe7 100644
--- a/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
+++ b/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,10 +80,10 @@ dimensionedScalar dynOneEqEddy::ce(const volSymmTensorField& D) const
         pow(KK + filter_(k_), 1.5)/(2*delta()) - filter_(pow(k_, 1.5))/delta();
 
     volScalarField ee =
-        2*delta()*ck(D)
+      2*delta()*ck(D)
        *(
-           filter_(sqrt(k_)*magSqr(D))
-         - 2*sqrt(KK + filter_(k_))*magSqr(filter_(D))
+            filter_(sqrt(k_)*magSqr(D))
+          - 2*sqrt(KK + filter_(k_))*magSqr(filter_(D))
         );
 
     dimensionedScalar mmmm = average(magSqr(mm));
@@ -135,8 +135,10 @@ dynOneEqEddy::dynOneEqEddy
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void dynOneEqEddy::correct(const tmp<volTensorField>& gradU)
+void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
 {
+    const volTensorField& gradU = tgradU();
+
     GenEddyVisc::correct(gradU);
 
     volSymmTensorField D = symm(gradU);
@@ -156,7 +158,7 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& gradU)
     kEqn.relax();
     kEqn.solve();
 
-    bound(k_, k0());
+    bound(k_, kMin_);
 
     updateSubGridScaleFields(D);
 }
diff --git a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
index a755d4e56e..86cd8227d1 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -263,8 +263,7 @@ kOmegaSSTSAS::kOmegaSSTSAS
         )
     ),
 
-    omega0_("omega0", dimless/dimTime, SMALL),
-    omegaSmall_("omegaSmall", dimless/dimTime, SMALL),
+    omegaMin_("omegaMin", dimless/dimTime, SMALL),
     y_(mesh_),
     Cmu_
     (
@@ -324,6 +323,8 @@ kOmegaSSTSAS::kOmegaSSTSAS
         mesh_
     )
 {
+    bound(omega_, omegaMin_);
+
     updateSubGridScaleFields(magSqr(symm(fvc::grad(U))));
 
     printCoeffs();
@@ -346,9 +347,8 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
 
     volVectorField gradK = fvc::grad(k_);
     volVectorField gradOmega = fvc::grad(omega_);
-    volScalarField L = sqrt(k_)/(pow025(Cmu_)*(omega_ + omegaSmall_));
-    volScalarField CDkOmega =
-        (2.0*alphaOmega2_)*(gradK & gradOmega)/(omega_ + omegaSmall_);
+    volScalarField L = sqrt(k_)/(pow025(Cmu_)*omega_);
+    volScalarField CDkOmega = (2.0*alphaOmega2_)*(gradK & gradOmega)/omega_;
     volScalarField F1 = this->F1(CDkOmega);
     volScalarField G = nuSgs_*2.0*S2;
 
@@ -368,14 +368,12 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
         kEqn.relax();
         kEqn.solve();
     }
-    bound(k_, k0());
+    bound(k_, kMin_);
 
     volScalarField grad_omega_k = max
     (
-        magSqr(gradOmega)/
-        sqr(omega_ + omegaSmall_),
-        magSqr(gradK)/
-        sqr(k_ + k0())
+        magSqr(gradOmega)/sqr(omega_),
+        magSqr(gradK)/sqr(k_)
     );
 
     // Turbulent frequency equation
@@ -397,7 +395,7 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
           + FSAS_
            *max
             (
-                dimensionedScalar("zero",dimensionSet(0, 0 , -2, 0, 0),0. ),
+                dimensionedScalar("zero",dimensionSet(0, 0, -2, 0, 0), 0.0),
                 zetaTilda2_*kappa_*S2*(L/Lvk2(S2))
               - 2.0/alphaPhi_*k_*grad_omega_k
             )
@@ -406,7 +404,7 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
         omegaEqn.relax();
         omegaEqn.solve();
     }
-    bound(omega_, omega0_);
+    bound(omega_, omegaMin_);
 
     updateSubGridScaleFields(S2);
 }
diff --git a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
index 89b9a88190..4c2ac8d55b 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,8 +103,7 @@ protected:
             dimensionedScalar zetaTilda2_;
             dimensionedScalar FSAS_;
 
-            dimensionedScalar omega0_;
-            dimensionedScalar omegaSmall_;
+            dimensionedScalar omegaMin_;
 
             wallDist y_;
             dimensionedScalar Cmu_;
diff --git a/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C b/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C
index f463e3f208..3decb7080e 100644
--- a/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C
+++ b/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -152,7 +152,7 @@ void locDynOneEqEddy::correct(const tmp<volTensorField>& gradU)
     kEqn.relax();
     kEqn.solve();
 
-    bound(k_, k0());
+    bound(k_, kMin_);
 
     updateSubGridScaleFields(D, KK);
 }
diff --git a/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C b/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C
index 28a0a1ed5a..96188e59e4 100644
--- a/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C
+++ b/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,7 +113,7 @@ void oneEqEddy::correct(const tmp<volTensorField>& gradU)
     kEqn.relax();
     kEqn.solve();
 
-    bound(k_, k0());
+    bound(k_, kMin_);
 
     updateSubGridScaleFields();
 }
diff --git a/src/turbulenceModels/incompressible/RAS/LRR/LRR.C b/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
index 08d1b4a39a..037da5b430 100644
--- a/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
+++ b/src/turbulenceModels/incompressible/RAS/LRR/LRR.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -186,9 +186,6 @@ LRR::LRR
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
-    nut_.correctBoundaryConditions();
-
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
     {
         FatalErrorIn
@@ -201,6 +198,11 @@ LRR::LRR
             << exit(FatalError);
     }
 
+    bound(epsilon_, epsilonMin_);
+
+    nut_ = Cmu_*sqr(k_)/epsilon_;
+    nut_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
@@ -320,7 +322,7 @@ void LRR::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Reynolds stress equation
@@ -368,15 +370,15 @@ void LRR::correct()
             R_.dimensions(),
             symmTensor
             (
-                k0_.value(), -GREAT, -GREAT,
-                k0_.value(), -GREAT,
-                k0_.value()
+                kMin_.value(), -GREAT, -GREAT,
+                kMin_.value(), -GREAT,
+                kMin_.value()
             )
         )
     );
 
     k_ = 0.5*tr(R_);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
index 1edcb224b5..b6d20aea2d 100644
--- a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -140,7 +140,9 @@ LamBremhorstKE::LamBremhorstKE
         autoCreateLowReNut("nut", mesh_)
     )
 {
-    nut_ = Cmu_*fMu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    bound(epsilon_, epsilonMin_);
+
+    nut_ = Cmu_*fMu_*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -260,7 +262,7 @@ void LamBremhorstKE::correct()
 
     epsEqn().relax();
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -276,7 +278,7 @@ void LamBremhorstKE::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
index 33939dfc15..98001c6a68 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -224,7 +224,9 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    bound(epsilon_, epsilonMin_);
+
+    nut_ = Cmu_*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
@@ -362,7 +364,7 @@ void LaunderGibsonRSTM::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Reynolds stress equation
@@ -419,15 +421,15 @@ void LaunderGibsonRSTM::correct()
             R_.dimensions(),
             symmTensor
             (
-                k0_.value(), -GREAT, -GREAT,
-                             k0_.value(), -GREAT,
-                                          k0_.value()
+                kMin_.value(), -GREAT, -GREAT,
+                kMin_.value(), -GREAT,
+                kMin_.value()
             )
         )
     );
 
     k_ == 0.5*tr(R_);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate turbulent viscosity
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
index 3bbdcb6a47..6f840c57d3 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -146,7 +146,9 @@ LaunderSharmaKE::LaunderSharmaKE
         autoCreateLowReNut("nut", mesh_)
     )
 {
-    nut_ = Cmu_*fMu()*sqr(k_)/(epsilonTilda_ + epsilonSmall_);
+    bound(epsilonTilda_, epsilonMin_);
+
+    nut_ = Cmu_*fMu()*sqr(k_)/epsilonTilda_;
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -256,7 +258,7 @@ void LaunderSharmaKE::correct()
 
     epsEqn().relax();
     solve(epsEqn);
-    bound(epsilonTilda_, epsilon0_);
+    bound(epsilonTilda_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -272,7 +274,7 @@ void LaunderSharmaKE::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
index c786116bd9..f01805f115 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -170,6 +170,8 @@ 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())))),
@@ -226,7 +228,9 @@ LienCubicKE::LienCubicKE
         )
     )
 {
-    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_) + C5viscosity_;
+    bound(epsilon_, epsilonMin_);
+
+    nut_ = Cmu_*sqr(k_)/epsilon_ + C5viscosity_;
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -350,7 +354,7 @@ void LienCubicKE::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -367,7 +371,7 @@ void LienCubicKE::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
index 886b984f57..8c7bbc7c7b 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
@@ -280,10 +280,12 @@ LienCubicKELowRe::LienCubicKELowRe
         )
     )
 {
+    bound(epsilon_, epsilonMin_);
+
     nut_ = Cmu_
       * (scalar(1) - exp(-Am_*yStar_))
       / (scalar(1) - exp(-Aepsilon_*yStar_) + SMALL)
-      * sqr(k_)/(epsilon_ + epsilonSmall_)
+      * sqr(k_)/epsilon_
         // cubic term C5, implicit part
       + max
         (
@@ -432,7 +434,7 @@ void LienCubicKELowRe::correct()
 #   include "wallDissipationI.H"
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -449,7 +451,7 @@ void LienCubicKELowRe::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
index 6e47e1d926..a2a807676d 100644
--- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -180,9 +180,11 @@ LienLeschzinerLowRe::LienLeschzinerLowRe
         autoCreateLowReNut("nut", mesh_)
     )
 {
+    bound(epsilon_, epsilonMin_);
+
     nut_ = Cmu_*(scalar(1) - exp(-Am_*yStar_))
        /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL)*sqr(k_)
-       /(epsilon_ + epsilonSmall_);
+       /(epsilon_);
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -322,7 +324,7 @@ void LienLeschzinerLowRe::correct()
     #include "wallDissipationI.H"
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -339,7 +341,7 @@ void LienLeschzinerLowRe::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
index f6d778b3a6..01324af90a 100644
--- a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
+++ b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -144,6 +144,7 @@ NonlinearKEShih::NonlinearKEShih
         )
     ),
 
+    //FIXME: should be named 'kappa_' or 'kappa'?
     kappa_
     (
         dimensioned<scalar>::lookupOrAddToDict
@@ -174,6 +175,7 @@ NonlinearKEShih::NonlinearKEShih
             IOobject::AUTO_WRITE
         ),
         mesh_
+        //FIXME: what about autoCreateK("k", mesh_)
     ),
 
     epsilon_
@@ -187,15 +189,32 @@ 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())))),
     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
     fEta_(A2_ + pow(eta_, 3.0)),
 
-    nut_("nut", Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_)),
+    // 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_
     (
@@ -215,6 +234,11 @@ NonlinearKEShih::NonlinearKEShih
         )
     )
 {
+    bound(epsilon_, epsilonMin_);
+
+    //FIXME: could use this
+    // nut_ = Cmu_*sqr(k_)/epsilon_;
+
     #include "wallNonlinearViscosityI.H"
 
     printCoeffs();
@@ -343,7 +367,7 @@ void NonlinearKEShih::correct()
     #include "wallDissipationI.H"
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -360,7 +384,7 @@ void NonlinearKEShih::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C b/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C
index 8944f94c23..f45c1ee0ff 100644
--- a/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C
+++ b/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C
@@ -80,11 +80,9 @@ RASModel::RASModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    k0_("k0", dimVelocity*dimVelocity, SMALL),
-    epsilon0_("epsilon0", k0_.dimensions()/dimTime, SMALL),
-    epsilonSmall_("epsilonSmall", epsilon0_.dimensions(), SMALL),
-    omega0_("omega0", dimless/dimTime, SMALL),
-    omegaSmall_("omegaSmall", omega0_.dimensions(), SMALL),
+    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
+    epsilonMin_("epsilonMin", kMin_.dimensions()/dimTime, SMALL),
+    omegaMin_("omegaMin", dimless/dimTime, SMALL),
 
     y_(mesh_)
 {
@@ -213,11 +211,9 @@ bool RASModel::read()
             coeffDict_ <<= *dictPtr;
         }
 
-        k0_.readIfPresent(*this);
-        epsilon0_.readIfPresent(*this);
-        epsilonSmall_.readIfPresent(*this);
-        omega0_.readIfPresent(*this);
-        omegaSmall_.readIfPresent(*this);
+        kMin_.readIfPresent(*this);
+        epsilonMin_.readIfPresent(*this);
+        omegaMin_.readIfPresent(*this);
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.H b/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.H
index a8eab758f0..337f044f0a 100644
--- a/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.H
+++ b/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.H
@@ -87,19 +87,13 @@ protected:
         dictionary coeffDict_;
 
         //- Lower limit of k
-        dimensionedScalar k0_;
+        dimensionedScalar kMin_;
 
         //- Lower limit of epsilon
-        dimensionedScalar epsilon0_;
-
-        //- Small epsilon value used to avoid divide by zero
-        dimensionedScalar epsilonSmall_;
+        dimensionedScalar epsilonMin_;
 
         //- Lower limit for omega
-        dimensionedScalar omega0_;
-
-        //- Small omega value used to avoid divide by zero
-        dimensionedScalar omegaSmall_;
+        dimensionedScalar omegaMin_;
 
         //- Near wall distance boundary field
         nearWallDist y_;
@@ -176,66 +170,40 @@ public:
 
         // Access
 
-            //- Return lower allowable limit for k
-            const dimensionedScalar& k0() const
+            //- Return lower allowable limit for k (default: SMALL)
+            const dimensionedScalar& kMin() const
             {
-                return k0_;
+                return kMin_;
             }
 
-            //- Return the lower allowable limit for epsilon
-            const dimensionedScalar& epsilon0() const
+            //- Return the lower allowable limit for epsilon (default: SMALL)
+            const dimensionedScalar& epsilonMin() const
             {
-                return epsilon0_;
+                return epsilonMin_;
             }
 
-            //- Return the value of epsilonSmall which is added to epsilon when
-            //  calculating nut
-            const dimensionedScalar& epsilonSmall() const
+            //- Return the lower allowable limit for omega (default: SMALL)
+            const dimensionedScalar& omegaMin() const
             {
-                return epsilonSmall_;
+                return omegaMin_;
             }
 
-            //- Return the lower allowable limit for omega0
-            const dimensionedScalar& omega0() const
+            //- Allow kMin to be changed
+            dimensionedScalar& kMin()
             {
-                return omega0_;
+                return kMin_;
             }
 
-            //- Return the value of omegaSmall which is added to omega when
-            //  calculating nut
-            const dimensionedScalar& omegaSmall() const
+            //- Allow epsilonMin to be changed
+            dimensionedScalar& epsilonMin()
             {
-                return omegaSmall_;
+                return epsilonMin_;
             }
 
-            //- Allow k0 to be changed
-            dimensionedScalar& k0()
+            //- Allow omegaMin to be changed
+            dimensionedScalar& omegaMin()
             {
-                return k0_;
-            }
-
-            //- Allow epsilon0 to be changed
-            dimensionedScalar& epsilon0()
-            {
-                return epsilon0_;
-            }
-
-            //- Allow epsilonSmall to be changed
-            dimensionedScalar& epsilonSmall()
-            {
-                return epsilonSmall_;
-            }
-
-            //- Allow omega0 to be changed
-            dimensionedScalar& omega0()
-            {
-                return omega0_;
-            }
-
-            //- Allow omegaSmall to be changed
-            dimensionedScalar& omegaSmall()
-            {
-                return omegaSmall_;
+                return omegaMin_;
             }
 
             //- Return the near wall distances
diff --git a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
index aa7e537001..1591608717 100644
--- a/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -155,7 +155,9 @@ RNGkEpsilon::RNGkEpsilon
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    bound(epsilon_, epsilonMin_);
+
+    nut_ = Cmu_*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -274,7 +276,7 @@ void RNGkEpsilon::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -290,7 +292,7 @@ void RNGkEpsilon::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
index 0861863673..8c41bf4d27 100644
--- a/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -128,7 +128,9 @@ kEpsilon::kEpsilon
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    bound(epsilon_, epsilonMin_);
+
+    nut_ = Cmu_*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -237,7 +239,7 @@ void kEpsilon::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -254,7 +256,7 @@ void kEpsilon::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
index 64b095d27b..04f505376b 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -137,7 +137,9 @@ kOmega::kOmega
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ = k_/(omega_ + omegaSmall_);
+    bound(omega_, omegaMin_);
+
+    nut_ = k_/omega_;
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -246,7 +248,7 @@ void kOmega::correct()
     omegaEqn().boundaryManipulate(omega_.boundaryField());
 
     solve(omegaEqn);
-    bound(omega_, omega0_);
+    bound(omega_, omegaMin_);
 
 
     // Turbulent kinetic energy equation
@@ -263,11 +265,11 @@ void kOmega::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
-    nut_ = k_/(omega_ + omegaSmall_);
+    nut_ = k_/omega_;
     nut_.correctBoundaryConditions();
 }
 
diff --git a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
index 8323a23011..4332fc3706 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
@@ -236,12 +236,14 @@ kOmegaSST::kOmegaSST
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(omega_, omegaMin_);
+
     nut_ =
     (
         a1_*k_
       / max
         (
-            a1_*(omega_ + omegaSmall_),
+            a1_*omega_,
             F2()*mag(symm(fvc::grad(U_)))
         )
     );
@@ -376,7 +378,7 @@ void kOmegaSST::correct()
     omegaEqn().boundaryManipulate(omega_.boundaryField());
 
     solve(omegaEqn);
-    bound(omega_, omega0_);
+    bound(omega_, omegaMin_);
 
     // Turbulent kinetic energy equation
     tmp<fvScalarMatrix> kEqn
@@ -392,7 +394,7 @@ void kOmegaSST::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity
diff --git a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
index 9b16f1c621..d72fd30f3e 100644
--- a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
+++ b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -126,6 +126,9 @@ qZeta::qZeta
         )
     ),
 
+    qMin_("qMin", dimVelocity, SMALL),
+    zetaMin_("zetaMin", dimVelocity/dimTime, SMALL),
+
     k_
     (
         IOobject
@@ -193,7 +196,11 @@ qZeta::qZeta
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ = Cmu_*fMu()*sqr(k_)/(epsilon_ + epsilonSmall_);
+    bound(epsilon_, epsilonMin_);
+    bound(q_, qMin_);
+    bound(zeta_, zetaMin_);
+
+    nut_ = Cmu_*fMu()*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -263,6 +270,9 @@ bool qZeta::read()
         sigmaZeta_.readIfPresent(coeffDict());
         anisotropic_.readIfPresent("anisotropic", coeffDict());
 
+        qMin_.readIfPresent(*this);
+        zetaMin_.readIfPresent(*this);
+
         return true;
     }
     else
@@ -302,7 +312,7 @@ void qZeta::correct()
 
     zetaEqn().relax();
     solve(zetaEqn);
-    bound(zeta_, epsilon0_/(2*sqrt(k0_)));
+    bound(zeta_, zetaMin_);
 
 
     // q equation
@@ -318,7 +328,7 @@ void qZeta::correct()
 
     qEqn().relax();
     solve(qEqn);
-    bound(q_, sqrt(k0_));
+    bound(q_, qMin_);
 
 
     // Re-calculate k and epsilon
diff --git a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H
index 6300868d2d..4c396c13a1 100644
--- a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H
+++ b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,6 +66,11 @@ class qZeta
             dimensionedScalar sigmaZeta_;
             Switch anisotropic_;
 
+        //- Lower limit of q
+        dimensionedScalar qMin_;
+
+        //- Lower limit of zeta
+        dimensionedScalar zetaMin_;
 
         // Fields
 
@@ -107,6 +112,33 @@ public:
 
     // Member Functions
 
+       // Access
+
+            //- Return lower allowable limit for q (default: SMALL)
+            const dimensionedScalar& qMin() const
+            {
+                return qMin_;
+            }
+
+            //- Return the lower allowable limit for zeta (default: SMALL)
+            const dimensionedScalar& zetaMin() const
+            {
+                return zetaMin_;
+            }
+
+            //- Allow qMin to be changed
+            dimensionedScalar& qMin()
+            {
+                return qMin_;
+            }
+
+            //- Allow zetaMin to be changed
+            dimensionedScalar& zetaMin()
+            {
+                return zetaMin_;
+            }
+
+
         //- Return the turbulence viscosity
         virtual tmp<volScalarField> nut() const
         {
diff --git a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
index 5072f9b1cf..5a53c6b0fc 100644
--- a/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
+++ b/src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,7 +69,7 @@ tmp<volScalarField> realizableKE::rCmu
     volScalarField As = sqrt(6.0)*cos(phis);
     volScalarField Us = sqrt(S2/2.0 + magSqr(skew(gradU)));
 
-    return 1.0/(A0_ + As*Us*k_/(epsilon_ + epsilonSmall_));
+    return 1.0/(A0_ + As*Us*k_/epsilon_);
 }
 
 
@@ -178,10 +178,10 @@ realizableKE::realizableKE
         autoCreateNut("nut", mesh_)
     )
 {
-    bound(k_, k0_);
-    bound(epsilon_, epsilon0_);
+    bound(k_, kMin_);
+    bound(epsilon_, epsilonMin_);
 
-    nut_ = rCmu(fvc::grad(U_))*sqr(k_)/(epsilon_ + epsilonSmall_);
+    nut_ = rCmu(fvc::grad(U_))*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -303,7 +303,7 @@ void realizableKE::correct()
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
-    bound(epsilon_, epsilon0_);
+    bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
@@ -319,7 +319,7 @@ void realizableKE::correct()
 
     kEqn().relax();
     solve(kEqn);
-    bound(k_, k0_);
+    bound(k_, kMin_);
 
 
     // Re-calculate viscosity