diff --git a/src/turbulenceModels/compressible/LRR/LRR.C b/src/turbulenceModels/compressible/LRR/LRR.C
index 380baf2d6d..ea9d71543b 100644
--- a/src/turbulenceModels/compressible/LRR/LRR.C
+++ b/src/turbulenceModels/compressible/LRR/LRR.C
@@ -55,17 +55,23 @@ LRR::LRR
 :
     turbulenceModel(typeName, rho, U, phi, thermophysicalModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    Clrr1(turbulenceModelCoeffs_.lookup("Clrr1")),
-    Clrr2(turbulenceModelCoeffs_.lookup("Clrr2")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    Cs(turbulenceModelCoeffs_.lookup("Cs")),
-    Ceps(turbulenceModelCoeffs_.lookup("Ceps")),
-    couplingFactor_(0.0),
-    alphaR(turbulenceModelCoeffs_.lookup("alphaR")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    alphah(turbulenceModelCoeffs_.lookup("alphah")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    Clrr1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clrr1", 1.8)),
+    Clrr2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clrr2", 0.6)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    Cs(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cs", 0.25)),
+    Ceps(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ceps", 0.15)),
+    couplingFactor_
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("couplingFactor", 0.0)
+    ),
+    alphaR(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaR", 1.22)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
+    alphah(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0)),
 
     R_
     (
@@ -121,21 +127,16 @@ LRR::LRR
 {
 #   include "wallViscosityI.H"
 
-    if (turbulenceModelCoeffs_.found("couplingFactor"))
+    if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
     {
-        turbulenceModelCoeffs_.lookup("couplingFactor") >> couplingFactor_;
-
-        if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
-        {
-            FatalErrorIn
-            (
-                "LRR::LRR"
-                "(const volVectorField& U, const surfaceScalarField& phi,"
-                "incompressibleTransportModel& lamTransportModel)"
-            )   << "couplingFactor = " << couplingFactor_
-                << " is not in range 0 - 1"
-                << exit(FatalError);
-        }
+        FatalErrorIn
+        (
+            "LRR::LRR"
+            "(const volVectorField& U, const surfaceScalarField& phi,"
+            "incompressibleTransportModel& lamTransportModel)"
+        )   << "couplingFactor = " << couplingFactor_
+            << " is not in range 0 - 1" << nl
+            << exit(FatalError);
     }
 }
 
@@ -191,24 +192,31 @@ bool LRR::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("Clrr1") >> Clrr1;
-        turbulenceModelCoeffs_.lookup("Clrr2") >> Clrr2;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("Cs") >> Cs;
-        turbulenceModelCoeffs_.lookup("Ceps") >> Ceps;
-        turbulenceModelCoeffs_.lookup("alphaR") >> alphaR;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("alphah") >> alphah;
-
-        turbulenceModelCoeffs_.lookup("couplingFactor") >> couplingFactor_;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        Clrr1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clrr1", 1.8);
+        Clrr2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clrr2", 0.6);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        Cs = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cs", 0.25);
+        Ceps = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ceps", 0.15);
+        alphaR = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaR", 1.22);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.76923
+            );
+        alphah = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0);
+        couplingFactor_ = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "couplingFactor",
+                0.0
+            );
 
         if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
         {
             FatalErrorIn("LRR::read()")
                 << "couplingFactor = " << couplingFactor_
-                << " is not in range 0 - 1"
+                << " is not in range 0 - 1" << nl
                 << exit(FatalError);
         }
 
diff --git a/src/turbulenceModels/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C
index da59204bd1..179f05cfae 100644
--- a/src/turbulenceModels/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C
+++ b/src/turbulenceModels/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C
@@ -57,19 +57,25 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
 :
     turbulenceModel(typeName, rho, U, phi, thermophysicalModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    Clg1(turbulenceModelCoeffs_.lookup("Clg1")),
-    Clg2(turbulenceModelCoeffs_.lookup("Clg2")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    Cs(turbulenceModelCoeffs_.lookup("Cs")),
-    Ceps(turbulenceModelCoeffs_.lookup("Ceps")),
-    C1Ref(turbulenceModelCoeffs_.lookup("C1Ref")),
-    C2Ref(turbulenceModelCoeffs_.lookup("C2Ref")),
-    couplingFactor_(0.0),
-    alphaR(turbulenceModelCoeffs_.lookup("alphaR")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    alphah(turbulenceModelCoeffs_.lookup("alphah")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    Clg1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clg1", 1.8)),
+    Clg2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clg2", 0.6)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    Cs(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cs", 0.25)),
+    Ceps(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ceps", 0.15)),
+    C1Ref(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1Ref", 0.5)),
+    C2Ref(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2Ref", 0.3)),
+    couplingFactor_
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("couplingFactor", 0.0)
+    ),
+    alphaR(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaR", 1.22)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
+    alphah(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0)),
 
     y_(mesh_),
 
@@ -127,21 +133,16 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
 {
 #   include "wallViscosityI.H"
 
-    if (turbulenceModelCoeffs_.found("couplingFactor"))
+    if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
     {
-        turbulenceModelCoeffs_.lookup("couplingFactor") >> couplingFactor_;
-
-        if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
-        {
-            FatalErrorIn
-            (
-                "LaunderGibsonRSTM::LaunderGibsonRSTM"
-                "(const volVectorField& U, const surfaceScalarField& phi,"
-                "incompressibleTransportModel& lamTransportModel)"
-            )   << "couplingFactor = " << couplingFactor_
-                << " is not in range 0 - 1"
-                << exit(FatalError);
-        }
+        FatalErrorIn
+        (
+            "LaunderGibsonRSTM::LaunderGibsonRSTM"
+            "(const volVectorField& U, const surfaceScalarField& phi,"
+            "incompressibleTransportModel& lamTransportModel)"
+        )   << "couplingFactor = " << couplingFactor_
+            << " is not in range 0 - 1" << nl
+            << exit(FatalError);
     }
 }
 
@@ -197,26 +198,34 @@ bool LaunderGibsonRSTM::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("Clg1") >> Clg1;
-        turbulenceModelCoeffs_.lookup("Clg2") >> Clg2;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("Cs") >> Cs;
-        turbulenceModelCoeffs_.lookup("Ceps") >> Ceps;
-        turbulenceModelCoeffs_.lookup("C1Ref") >> C1Ref;
-        turbulenceModelCoeffs_.lookup("C2Ref") >> C2Ref;
-        turbulenceModelCoeffs_.lookup("alphaR") >> alphaR;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("alphah") >> alphah;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        Clg1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clg1", 1.8);
+        Clg2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clg2", 0.6);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        Cs = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cs", 0.25);
+        Ceps = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ceps", 0.15);
+        C1Ref = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1Ref", 0.5);
+        C2Ref = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2Ref", 0.3);
+        alphaR = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaR", 1.22);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.76923
+            );
+        alphah = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0);
 
-        turbulenceModelCoeffs_.lookup("couplingFactor") >> couplingFactor_;
+        couplingFactor_ = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "couplingFactor",
+                0.0
+            );
 
         if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
         {
             FatalErrorIn("LaunderGibsonRSTM::read()")
                 << "couplingFactor = " << couplingFactor_
-                << " is not in range 0 - 1"
+                << " is not in range 0 - 1" << nl
                 << exit(FatalError);
         }
 
diff --git a/src/turbulenceModels/compressible/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/compressible/LaunderSharmaKE/LaunderSharmaKE.C
index 7e758e9d9d..f5328518ff 100644
--- a/src/turbulenceModels/compressible/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/compressible/LaunderSharmaKE/LaunderSharmaKE.C
@@ -52,7 +52,7 @@ tmp<volScalarField> LaunderSharmaKE::fMu() const
 
 tmp<volScalarField> LaunderSharmaKE::f2() const
 {
-    return 
+    return
         scalar(1)
       - 0.3*exp(-min(sqr(rho_*sqr(k_)/(mu()*epsilon_)), scalar(50.0)));
 }
@@ -71,13 +71,16 @@ LaunderSharmaKE::LaunderSharmaKE
 :
     turbulenceModel(typeName, rho, U, phi, thermophysicalModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    C3(turbulenceModelCoeffs_.lookup("C3")),
-    alphak(turbulenceModelCoeffs_.lookup("alphak")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    alphah(turbulenceModelCoeffs_.lookup("alphah")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    C3(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C3", -0.33)),
+    alphak(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
+    alphah(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0)),
 
     k_
     (
@@ -176,12 +179,16 @@ bool LaunderSharmaKE::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("C3") >> C3;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("alphah") >> alphah;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        C3 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C3", -0.33);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.76923
+            );
+        alphah = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0);
 
         return true;
     }
diff --git a/src/turbulenceModels/compressible/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/compressible/RNGkEpsilon/RNGkEpsilon.C
index c37b3632e8..0abe4b55f0 100644
--- a/src/turbulenceModels/compressible/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/turbulenceModels/compressible/RNGkEpsilon/RNGkEpsilon.C
@@ -55,15 +55,15 @@ RNGkEpsilon::RNGkEpsilon
 :
     turbulenceModel(typeName, rho, U, phi, thermophysicalModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    C3(turbulenceModelCoeffs_.lookup("C3")),
-    alphak(turbulenceModelCoeffs_.lookup("alphak")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    alphah(turbulenceModelCoeffs_.lookup("alphah")),
-    eta0(turbulenceModelCoeffs_.lookup("eta0")),
-    beta(turbulenceModelCoeffs_.lookup("beta")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.0845)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.42)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.68)),
+    C3(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C3", -0.33)),
+    alphak(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.39)),
+    alphaEps(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 1.39)),
+    alphah(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0)),
+    eta0(turbulenceModelCoeffs_.lookupOrDefault<scalar>("eta0", 4.38)),
+    beta(turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta", 0.012)),
 
     k_
     (
@@ -164,15 +164,19 @@ bool RNGkEpsilon::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("C3") >> C3;
-        turbulenceModelCoeffs_.lookup("alphak") >> alphak;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("alphah") >> alphah;
-        turbulenceModelCoeffs_.lookup("eta0") >> eta0;
-        turbulenceModelCoeffs_.lookup("beta") >> beta;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.0845);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.42);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.68);
+        C3 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C3", -0.33);
+        alphak = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.39);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                1.39
+            );
+        alphah = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0);
+        eta0 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("eta0", 4.38);
+        beta = turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta", 0.012);
 
         return true;
     }
@@ -210,7 +214,7 @@ void RNGkEpsilon::correct()
     volScalarField eta = sqrt(mag(S2))*k_/epsilon_;
     volScalarField eta3 = eta*sqr(eta);
 
-    volScalarField R = 
+    volScalarField R =
         ((eta*(-eta/eta0 + scalar(1)))/(beta*eta3 + scalar(1)));
 
 #   include "wallFunctionsI.H"
diff --git a/src/turbulenceModels/compressible/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/compressible/SpalartAllmaras/SpalartAllmaras.C
index a2ccbf80a4..265191af0a 100644
--- a/src/turbulenceModels/compressible/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/compressible/SpalartAllmaras/SpalartAllmaras.C
@@ -115,16 +115,16 @@ SpalartAllmaras::SpalartAllmaras
 :
     turbulenceModel(typeName, rho, U, phi, thermophysicalModel),
 
-    alphaNut(turbulenceModelCoeffs_.lookup("alphaNut")),
-    alphah(turbulenceModelCoeffs_.lookup("alphah")),
+    alphaNut(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaNut", 1.5)),
+    alphah(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0)),
 
-    Cb1(turbulenceModelCoeffs_.lookup("Cb1")),
-    Cb2(turbulenceModelCoeffs_.lookup("Cb2")),
+    Cb1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cb1", 0.1355)),
+    Cb2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cb2", 0.622)),
     Cw1(Cb1/sqr(kappa_) + alphaNut*(1.0 + Cb2)),
-    Cw2(turbulenceModelCoeffs_.lookup("Cw2")),
-    Cw3(turbulenceModelCoeffs_.lookup("Cw3")),
-    Cv1(turbulenceModelCoeffs_.lookup("Cv1")),
-    Cv2(turbulenceModelCoeffs_.lookup("Cv2")),
+    Cw2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cw2", 0.3)),
+    Cw3(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cw3", 2.0)),
+    Cv1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cv1", 7.1)),
+    Cv2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cv2", 5.0)),
 
     nuTilda_
     (
@@ -214,16 +214,20 @@ bool SpalartAllmaras::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("alphaNut") >> alphaNut;
-        turbulenceModelCoeffs_.lookup("alphah") >> alphah;
+        alphaNut = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaNut",
+                1.5
+            );
+        alphah = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0);
 
-        turbulenceModelCoeffs_.lookup("Cb1") >> Cb1;
-        turbulenceModelCoeffs_.lookup("Cb2") >> Cb2;
+        Cb1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cb1", 0.1355);
+        Cb2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cb2", 0.622);
         Cw1 = Cb1/sqr(kappa_) + alphaNut*(1.0 + Cb2);
-        turbulenceModelCoeffs_.lookup("Cw2") >> Cw2;
-        turbulenceModelCoeffs_.lookup("Cw3") >> Cw3;
-        turbulenceModelCoeffs_.lookup("Cv1") >> Cv1;
-        turbulenceModelCoeffs_.lookup("Cv2") >> Cv2;
+        Cw2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cw2", 0.3);
+        Cw3 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cw3", 2.0);
+        Cv1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cv1", 7.1);
+        Cv2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cv2", 5.0);
 
         return true;
     }
diff --git a/src/turbulenceModels/compressible/kEpsilon/kEpsilon.C b/src/turbulenceModels/compressible/kEpsilon/kEpsilon.C
index a4438c23ee..02e863550c 100644
--- a/src/turbulenceModels/compressible/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/compressible/kEpsilon/kEpsilon.C
@@ -54,13 +54,16 @@ kEpsilon::kEpsilon
 :
     turbulenceModel(typeName, rho, U, phi, thermophysicalModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    C3(turbulenceModelCoeffs_.lookup("C3")),
-    alphak(turbulenceModelCoeffs_.lookup("alphak")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    alphah(turbulenceModelCoeffs_.lookup("alphah")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    C3(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C3", -0.33)),
+    alphak(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
+    alphah(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0)),
 
     k_
     (
@@ -162,13 +165,17 @@ bool kEpsilon::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("C3") >> C3;
-        turbulenceModelCoeffs_.lookup("alphak") >> alphak;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("alphah") >> alphah;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        C3 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C3", -0.33);
+        alphak = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.76923
+            );
+        alphah = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0);
 
         return true;
     }
@@ -232,7 +239,7 @@ void kEpsilon::correct()
       + fvm::div(phi_, k_)
       - fvm::laplacian(DkEff(), k_)
      ==
-        G 
+        G
       - fvm::SuSp((2.0/3.0)*rho_*divU, k_)
       - fvm::Sp(rho_*epsilon_/k_, k_)
     );
diff --git a/src/turbulenceModels/compressible/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/compressible/kOmegaSST/kOmegaSST.C
index 4956b162be..8592f98003 100644
--- a/src/turbulenceModels/compressible/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/compressible/kOmegaSST/kOmegaSST.C
@@ -97,23 +97,32 @@ kOmegaSST::kOmegaSST
 :
     turbulenceModel(typeName, rho, U, phi, thermophysicalModel),
 
-    alphaK1(turbulenceModelCoeffs_.lookup("alphaK1")),
-    alphaK2(turbulenceModelCoeffs_.lookup("alphaK2")),
-    alphaOmega1(turbulenceModelCoeffs_.lookup("alphaOmega1")),
-    alphaOmega2(turbulenceModelCoeffs_.lookup("alphaOmega2")),
-    alphah(turbulenceModelCoeffs_.lookup("alphah")),
-    gamma1(turbulenceModelCoeffs_.lookup("gamma1")),
-    gamma2(turbulenceModelCoeffs_.lookup("gamma2")),
-    beta1(turbulenceModelCoeffs_.lookup("beta1")),
-    beta2(turbulenceModelCoeffs_.lookup("beta2")),
-    betaStar(turbulenceModelCoeffs_.lookup("betaStar")),
-    a1(turbulenceModelCoeffs_.lookup("a1")),
-    c1(turbulenceModelCoeffs_.lookup("c1")),
+    alphaK1
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaK1", 0.85034)
+    ),
+    alphaK2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaK2", 1.0)),
+    alphaOmega1
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaOmega1", 0.5)
+    ),
+    alphaOmega2
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaOmega2", 0.85616)
+    ),
+    alphah(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0)),
+    gamma1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("gamma1", 0.5532)),
+    gamma2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("gamma2", 0.4403)),
+    beta1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta1", 0.075)),
+    beta2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta2", 0.0828)),
+    betaStar(turbulenceModelCoeffs_.lookupOrDefault<scalar>("betaStar", 0.09)),
+    a1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("a1", 0.31)),
+    c1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("c1", 10.0)),
 
     omega0_("omega0", dimless/dimTime, SMALL),
     omegaSmall_("omegaSmall", dimless/dimTime, SMALL),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
 
     y_(mesh_),
 
@@ -216,19 +225,43 @@ bool kOmegaSST::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("alphaK1") >> alphaK1;
-        turbulenceModelCoeffs_.lookup("alphaK2") >> alphaK2;
-        turbulenceModelCoeffs_.lookup("alphaOmega1") >> alphaOmega1;
-        turbulenceModelCoeffs_.lookup("alphaOmega2") >> alphaOmega2;
-        turbulenceModelCoeffs_.lookup("alphah") >> alphah;
-        turbulenceModelCoeffs_.lookup("gamma1") >> gamma1;
-        turbulenceModelCoeffs_.lookup("gamma2") >> gamma2;
-        turbulenceModelCoeffs_.lookup("beta1") >> beta1;
-        turbulenceModelCoeffs_.lookup("beta2") >> beta2;
-        turbulenceModelCoeffs_.lookup("betaStar") >> betaStar;
-        turbulenceModelCoeffs_.lookup("a1") >> a1;
-        turbulenceModelCoeffs_.lookup("c1") >> c1;
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
+        alphaK1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaK1",
+                0.85034
+            );
+        alphaK2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaK2", 1.0);
+        alphaOmega1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaOmega1",
+                0.5
+            );
+        alphaOmega2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaOmega2",
+                0.85616
+            );
+        alphah = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0);
+        gamma1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "gamma1",
+                0.5532
+            );
+        gamma2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "gamma2",
+                0.4403
+            );
+        beta1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta1", 0.075);
+        beta2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta2", 0.0828);
+        betaStar = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "betaStar",
+                0.09
+            );
+        a1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("a1", 0.31);
+        c1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("c1", 10.0);
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
 
         return true;
     }
@@ -244,7 +277,7 @@ void kOmegaSST::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ = 
+        mut_ =
             a1*rho_*k_/max(a1*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_)))));
 #       include "kOmegaWallViscosityI.H"
         return;
diff --git a/src/turbulenceModels/compressible/realizableKE/realizableKE.C b/src/turbulenceModels/compressible/realizableKE/realizableKE.C
index 421fc4db4a..85ebd9a583 100644
--- a/src/turbulenceModels/compressible/realizableKE/realizableKE.C
+++ b/src/turbulenceModels/compressible/realizableKE/realizableKE.C
@@ -54,7 +54,7 @@ tmp<volScalarField> realizableKE::rCmu
     tmp<volSymmTensorField> tS = dev(symm(gradU));
     const volSymmTensorField& S = tS();
 
-    volScalarField W = 
+    volScalarField W =
         (2*sqrt(2.0))*((S&S)&&S)
        /(
             magS*S2
@@ -63,7 +63,7 @@ tmp<volScalarField> realizableKE::rCmu
 
     tS.clear();
 
-    volScalarField phis = 
+    volScalarField phis =
         (1.0/3.0)*acos(min(max(sqrt(6.0)*W, -scalar(1)), scalar(1)));
     volScalarField As = sqrt(6.0)*cos(phis);
     volScalarField Us = sqrt(S2/2.0 + magSqr(skew(gradU)));
@@ -94,13 +94,16 @@ realizableKE::realizableKE
 )
 :
     turbulenceModel(typeName, rho, U, phi, thermophysicalModel),
-    
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    A0(turbulenceModelCoeffs_.lookup("A0")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    alphak(turbulenceModelCoeffs_.lookup("alphak")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    alphah(turbulenceModelCoeffs_.lookup("alphah")),
+
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    A0(turbulenceModelCoeffs_.lookupOrDefault<scalar>("A0", 4.0)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.9)),
+    alphak(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.833333)
+    ),
+    alphah(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0)),
 
     k_
     (
@@ -203,12 +206,16 @@ bool realizableKE::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("A0") >> A0;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("alphak") >> alphak;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("alphah") >> alphah;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        A0 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("A0", 4.0);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.9);
+        alphak = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.833333
+            );
+        alphah = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphah", 1.0);
 
         return true;
     }
@@ -240,10 +247,10 @@ void realizableKE::correct()
     volTensorField gradU = fvc::grad(U_);
     volScalarField S2 = 2*magSqr(dev(symm(gradU)));
     volScalarField magS = sqrt(S2);
-    
+
     volScalarField eta = magS*k_/epsilon_;
     volScalarField C1 = max(eta/(scalar(5) + eta), scalar(0.43));
-    
+
     volScalarField G = mut_*(gradU && dev(twoSymm(gradU)));
 
 #   include "wallFunctionsI.H"
@@ -286,7 +293,7 @@ void realizableKE::correct()
     kEqn().relax();
     solve(kEqn);
     bound(k_, k0_);
-    
+
     // Re-calculate viscosity
     mut_ = rCmu(gradU, S2, magS)*rho_*sqr(k_)/epsilon_;
 
diff --git a/src/turbulenceModels/incompressible/LRR/LRR.C b/src/turbulenceModels/incompressible/LRR/LRR.C
index 1c37f360f7..b672516b1f 100644
--- a/src/turbulenceModels/incompressible/LRR/LRR.C
+++ b/src/turbulenceModels/incompressible/LRR/LRR.C
@@ -52,15 +52,21 @@ LRR::LRR
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    Clrr1(turbulenceModelCoeffs_.lookup("Clrr1")),
-    Clrr2(turbulenceModelCoeffs_.lookup("Clrr2")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    Cs(turbulenceModelCoeffs_.lookup("Cs")),
-    Ceps(turbulenceModelCoeffs_.lookup("Ceps")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    couplingFactor_(0.0),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    Clrr1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clrr1", 1.8)),
+    Clrr2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clrr2", 0.6)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    Cs(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cs", 0.25)),
+    Ceps(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ceps", 0.15)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
+    couplingFactor_
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("couplingFactor", 0.0)
+    ),
 
     R_
     (
@@ -105,21 +111,16 @@ LRR::LRR
 {
 #   include "wallViscosityI.H"
 
-    if (turbulenceModelCoeffs_.found("couplingFactor"))
+    if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
     {
-        turbulenceModelCoeffs_.lookup("couplingFactor") >> couplingFactor_;
-
-        if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
-        {
-            FatalErrorIn
-            (
-                "LRR::LRR"
-                "(const volVectorField& U, const surfaceScalarField& phi,"
-                "transportModel& lamTransportModel)"
-            )   << "couplingFactor = " << couplingFactor_
-                << " is not in range 0 - 1"
-                << exit(FatalError);
-        }
+        FatalErrorIn
+        (
+            "LRR::LRR"
+            "(const volVectorField& U, const surfaceScalarField& phi,"
+            "transportModel& lamTransportModel)"
+        )   << "couplingFactor = " << couplingFactor_
+            << " is not in range 0 - 1" << nl
+            << exit(FatalError);
     }
 }
 
@@ -173,16 +174,24 @@ bool LRR::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("Clrr1") >> Clrr1;
-        turbulenceModelCoeffs_.lookup("Clrr2") >> Clrr2;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("Cs") >> Cs;
-        turbulenceModelCoeffs_.lookup("Ceps") >> Ceps;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        Clrr1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clrr1", 1.8);
+        Clrr2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clrr2", 0.6);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        Cs = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cs", 0.25);
+        Ceps = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ceps", 0.15);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.76923
+            );
 
-        turbulenceModelCoeffs_.lookup("couplingFactor") >> couplingFactor_;
+        couplingFactor_ = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "couplingFactor",
+                0.0
+            );
 
         if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
         {
diff --git a/src/turbulenceModels/incompressible/LamBremhorstKE/LamBremhorstKE.C b/src/turbulenceModels/incompressible/LamBremhorstKE/LamBremhorstKE.C
index bb506dc999..d89072281f 100644
--- a/src/turbulenceModels/incompressible/LamBremhorstKE/LamBremhorstKE.C
+++ b/src/turbulenceModels/incompressible/LamBremhorstKE/LamBremhorstKE.C
@@ -51,10 +51,13 @@ LamBremhorstKE::LamBremhorstKE
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
 
     k_
     (
@@ -153,10 +156,14 @@ bool LamBremhorstKE::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+        (
+            "alphaEps",
+            0.76923
+        );
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/incompressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C
index b01febba7b..27a16dc4c8 100644
--- a/src/turbulenceModels/incompressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C
+++ b/src/turbulenceModels/incompressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C
@@ -52,18 +52,24 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    Clg1(turbulenceModelCoeffs_.lookup("Clg1")),
-    Clg2(turbulenceModelCoeffs_.lookup("Clg2")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    Cs(turbulenceModelCoeffs_.lookup("Cs")),
-    Ceps(turbulenceModelCoeffs_.lookup("Ceps")),
-    alphaR(turbulenceModelCoeffs_.lookup("alphaR")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    C1Ref(turbulenceModelCoeffs_.lookup("C1Ref")),
-    C2Ref(turbulenceModelCoeffs_.lookup("C2Ref")),
-    couplingFactor_(0.0),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    Clg1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clg1", 1.8)),
+    Clg2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clg2", 0.6)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    Cs(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cs", 0.25)),
+    Ceps(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ceps", 0.15)),
+    alphaR(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaR", 1.22)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
+    C1Ref(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1Ref", 0.5)),
+    C2Ref(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2Ref", 0.3)),
+    couplingFactor_
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("couplingFactor", 0.0)
+    ),
 
     yr_(mesh_),
 
@@ -110,21 +116,16 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
 {
 #   include "wallViscosityI.H"
 
-    if (turbulenceModelCoeffs_.found("couplingFactor"))
+    if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
     {
-        turbulenceModelCoeffs_.lookup("couplingFactor") >> couplingFactor_;
-
-        if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
-        {
-            FatalErrorIn
-            (
-                "LaunderGibsonRSTM::LaunderGibsonRSTM"
-                "(const volVectorField& U, const surfaceScalarField& phi,"
-                "transportModel& lamTransportModel)"
-            )   << "couplingFactor = " << couplingFactor_
-                << " is not in range 0 - 1"
-                << exit(FatalError);
-        }
+        FatalErrorIn
+        (
+            "LaunderGibsonRSTM::LaunderGibsonRSTM"
+            "(const volVectorField& U, const surfaceScalarField& phi,"
+            "transportModel& lamTransportModel)"
+        )   << "couplingFactor = " << couplingFactor_
+            << " is not in range 0 - 1" << nl
+            << exit(FatalError);
     }
 }
 
@@ -178,19 +179,27 @@ bool LaunderGibsonRSTM::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("Clg1") >> Clg1;
-        turbulenceModelCoeffs_.lookup("Clg2") >> Clg2;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("Cs") >> Cs;
-        turbulenceModelCoeffs_.lookup("Ceps") >> Ceps;
-        turbulenceModelCoeffs_.lookup("alphaR") >> alphaR;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("C1Ref") >> C1Ref;
-        turbulenceModelCoeffs_.lookup("C2Ref") >> C2Ref;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        Clg1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clg1", 1.8);
+        Clg2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Clg2", 0.6);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        Cs = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cs", 0.25);
+        Ceps = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ceps", 0.15);
+        alphaR = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaR", 1.22);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.76923
+            );
+        C1Ref = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1Ref", 0.5);
+        C2Ref = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2Ref", 0.3);
 
-        turbulenceModelCoeffs_.lookup("couplingFactor") >> couplingFactor_;
+        couplingFactor_ = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "couplingFactor",
+                0.0
+            );
 
         if (couplingFactor_ < 0.0 || couplingFactor_ > 1.0)
         {
diff --git a/src/turbulenceModels/incompressible/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/incompressible/LaunderSharmaKE/LaunderSharmaKE.C
index f6d08178b3..015e92a031 100644
--- a/src/turbulenceModels/incompressible/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/incompressible/LaunderSharmaKE/LaunderSharmaKE.C
@@ -66,10 +66,13 @@ LaunderSharmaKE::LaunderSharmaKE
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
 
     k_
     (
@@ -158,10 +161,14 @@ bool LaunderSharmaKE::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.76923
+            );
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/LienCubicKE/LienCubicKE.C b/src/turbulenceModels/incompressible/LienCubicKE/LienCubicKE.C
index a64913bda7..7de1b1d374 100644
--- a/src/turbulenceModels/incompressible/LienCubicKE/LienCubicKE.C
+++ b/src/turbulenceModels/incompressible/LienCubicKE/LienCubicKE.C
@@ -52,16 +52,19 @@ LienCubicKE::LienCubicKE
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    alphak(turbulenceModelCoeffs_.lookup("alphak")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    A1(turbulenceModelCoeffs_.lookup("A1")),
-    A2(turbulenceModelCoeffs_.lookup("A2")),
-    Ctau1(turbulenceModelCoeffs_.lookup("Ctau1")),
-    Ctau2(turbulenceModelCoeffs_.lookup("Ctau2")),
-    Ctau3(turbulenceModelCoeffs_.lookup("Ctau3")),
-    alphaKsi(turbulenceModelCoeffs_.lookup("alphaKsi")),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    alphak(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
+    A1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("A1", 1.25)),
+    A2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("A2", 1000.0)),
+    Ctau1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau1", -4.0)),
+    Ctau2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau2", 13.0)),
+    Ctau3(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau3", -2.0)),
+    alphaKsi(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaKsi", 0.9)),
 
     k_
     (
@@ -195,16 +198,24 @@ bool LienCubicKE::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("alphak") >> alphak;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("A1") >> A1;
-        turbulenceModelCoeffs_.lookup("A2") >> A2;
-        turbulenceModelCoeffs_.lookup("Ctau1") >> Ctau1;
-        turbulenceModelCoeffs_.lookup("Ctau2") >> Ctau2;
-        turbulenceModelCoeffs_.lookup("Ctau3") >> Ctau3;
-        turbulenceModelCoeffs_.lookup("alphaKsi") >> alphaKsi;
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        alphak = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.76923
+            );
+        A1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("A1", 1.25);
+        A2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("A2", 1000.0);
+        Ctau1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau1", -4.0);
+        Ctau2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau2", 13.0);
+        Ctau3 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau3", -2.0);
+        alphaKsi = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaKsi",
+                0.9
+            );
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/LienCubicKELowRe/LienCubicKELowRe.C b/src/turbulenceModels/incompressible/LienCubicKELowRe/LienCubicKELowRe.C
index b72ee3cf14..2138565629 100644
--- a/src/turbulenceModels/incompressible/LienCubicKELowRe/LienCubicKELowRe.C
+++ b/src/turbulenceModels/incompressible/LienCubicKELowRe/LienCubicKELowRe.C
@@ -52,20 +52,23 @@ LienCubicKELowRe::LienCubicKELowRe
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    alphak(turbulenceModelCoeffs_.lookup("alphak")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    A1(turbulenceModelCoeffs_.lookup("A1")),
-    A2(turbulenceModelCoeffs_.lookup("A2")),
-    Ctau1(turbulenceModelCoeffs_.lookup("Ctau1")),
-    Ctau2(turbulenceModelCoeffs_.lookup("Ctau2")),
-    Ctau3(turbulenceModelCoeffs_.lookup("Ctau3")),
-    alphaKsi(turbulenceModelCoeffs_.lookup("alphaKsi")),
-    CmuWall(turbulenceModelCoeffs_.lookup("Cmu")),
-    Am(turbulenceModelCoeffs_.lookup("Am")),
-    Aepsilon(turbulenceModelCoeffs_.lookup("Aepsilon")),
-    Amu(turbulenceModelCoeffs_.lookup("Amu")),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    alphak(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
+    A1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("A1", 1.25)),
+    A2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("A2", 1000.0)),
+    Ctau1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau1", -4.0)),
+    Ctau2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau2", 13.0)),
+    Ctau3(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau3", -2.0)),
+    alphaKsi(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaKsi", 0.9)),
+    CmuWall(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    Am(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Am", 0.016)),
+    Aepsilon(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Aepsilon", 0.263)),
+    Amu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Amu", 0.00222)),
 
     k_
     (
@@ -220,20 +223,32 @@ bool LienCubicKELowRe::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("alphak") >> alphak;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("A1") >> A1;
-        turbulenceModelCoeffs_.lookup("A2") >> A2;
-        turbulenceModelCoeffs_.lookup("Ctau1") >> Ctau1;
-        turbulenceModelCoeffs_.lookup("Ctau2") >> Ctau2;
-        turbulenceModelCoeffs_.lookup("Ctau3") >> Ctau3;
-        turbulenceModelCoeffs_.lookup("alphaKsi") >> alphaKsi;
-        turbulenceModelCoeffs_.lookup("Cmu") >> CmuWall;
-        turbulenceModelCoeffs_.lookup("Am") >> Am;
-        turbulenceModelCoeffs_.lookup("Aepsilon") >> Aepsilon;
-        turbulenceModelCoeffs_.lookup("Amu") >> Amu;
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        alphak = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.76923
+            );
+        A1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("A1", 1.25);
+        A2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("A2", 1000.0);
+        Ctau1= turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau1", -4.0);
+        Ctau2= turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau2", 13.0);
+        Ctau3 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau3", -2.0);
+        alphaKsi = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaKsi",
+                0.9
+            );
+        CmuWall= turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        Am = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Am", 0.016);
+        Aepsilon = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "Aepsilon",
+                0.263
+            );
+        Amu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Amu", 0.00222);
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/LienLeschzinerLowRe/LienLeschzinerLowRe.C b/src/turbulenceModels/incompressible/LienLeschzinerLowRe/LienLeschzinerLowRe.C
index 1dbba5a001..7a8df3e45d 100644
--- a/src/turbulenceModels/incompressible/LienLeschzinerLowRe/LienLeschzinerLowRe.C
+++ b/src/turbulenceModels/incompressible/LienLeschzinerLowRe/LienLeschzinerLowRe.C
@@ -52,14 +52,17 @@ LienLeschzinerLowRe::LienLeschzinerLowRe
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    alphak(turbulenceModelCoeffs_.lookup("alphak")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    Am(turbulenceModelCoeffs_.lookup("Am")),
-    Aepsilon(turbulenceModelCoeffs_.lookup("Aepsilon")),
-    Amu(turbulenceModelCoeffs_.lookup("Amu")),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    alphak(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    Am(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Am", 0.016)),
+    Aepsilon(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Aepsilon", 0.263)),
+    Amu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Amu", 0.00222)),
 
     k_
     (
@@ -160,14 +163,22 @@ bool LienLeschzinerLowRe::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("alphak") >> alphak;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("Am") >> Am;
-        turbulenceModelCoeffs_.lookup("Aepsilon") >> Aepsilon;
-        turbulenceModelCoeffs_.lookup("Amu") >> Amu;
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        alphak = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.76923
+            );
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        Am = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Am", 0.016);
+        Aepsilon = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "Aepsilon",
+                0.263
+            );
+        Amu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Amu", 0.00222);
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/NonlinearKEShih/NonlinearKEShih.C b/src/turbulenceModels/incompressible/NonlinearKEShih/NonlinearKEShih.C
index 6d6d904c62..4047412037 100644
--- a/src/turbulenceModels/incompressible/NonlinearKEShih/NonlinearKEShih.C
+++ b/src/turbulenceModels/incompressible/NonlinearKEShih/NonlinearKEShih.C
@@ -52,16 +52,19 @@ NonlinearKEShih::NonlinearKEShih
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    alphak(turbulenceModelCoeffs_.lookup("alphak")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    A1(turbulenceModelCoeffs_.lookup("A1")),
-    A2(turbulenceModelCoeffs_.lookup("A2")),
-    Ctau1(turbulenceModelCoeffs_.lookup("Ctau1")),
-    Ctau2(turbulenceModelCoeffs_.lookup("Ctau2")),
-    Ctau3(turbulenceModelCoeffs_.lookup("Ctau3")),
-    alphaKsi(turbulenceModelCoeffs_.lookup("alphaKsi")),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    alphak(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
+    A1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("A1", 1.25)),
+    A2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("A2", 1000.0)),
+    Ctau1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau1", -4.0)),
+    Ctau2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau2", 13.0)),
+    Ctau3(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau3", -2.0)),
+    alphaKsi(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaKsi", 0.9)),
 
     k_
     (
@@ -179,16 +182,24 @@ bool NonlinearKEShih::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("alphak") >> alphak;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("A1") >> A1;
-        turbulenceModelCoeffs_.lookup("A2") >> A2;
-        turbulenceModelCoeffs_.lookup("Ctau1") >> Ctau1;
-        turbulenceModelCoeffs_.lookup("Ctau2") >> Ctau2;
-        turbulenceModelCoeffs_.lookup("Ctau3") >> Ctau3;
-        turbulenceModelCoeffs_.lookup("alphaKsi") >> alphaKsi;
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        alphak = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.76923
+            );
+        A1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("A1", 1.25);
+        A2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("A2", 1000.0);
+        Ctau1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau1", -4.0);
+        Ctau2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau2", 13.0);
+        Ctau3 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Ctau3", -2.0);
+        alphaKsi = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaKsi",
+                0.9
+            );
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/QZeta/QZeta.C b/src/turbulenceModels/incompressible/QZeta/QZeta.C
index dcafe71f59..b4428c7934 100644
--- a/src/turbulenceModels/incompressible/QZeta/QZeta.C
+++ b/src/turbulenceModels/incompressible/QZeta/QZeta.C
@@ -51,7 +51,7 @@ tmp<volScalarField> QZeta::fMu() const
     }
     else
     {
-        return 
+        return
             exp(-6.0/sqr(scalar(1) + Rt/50.0))
            *(scalar(1) + 3.0*exp(-Rt/10.0));
     }
@@ -77,11 +77,17 @@ QZeta::QZeta
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    alphaZeta(turbulenceModelCoeffs_.lookup("alphaZeta")),
-    Anisotropic(turbulenceModelCoeffs_.lookup("anisotropic")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    alphaZeta
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaZeta", 0.76923)
+    ),
+    Anisotropic
+    (
+         turbulenceModelCoeffs_.lookupOrDefault<Switch>("anisotropic", false)
+    ),
 
     k_
     (
@@ -198,11 +204,19 @@ bool QZeta::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("alphaZeta") >> alphaZeta;
-        turbulenceModelCoeffs_.lookup("anisotropic") >> Anisotropic;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.92);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.44);
+        alphaZeta = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaZeta",
+                0.76923
+            );
+        Anisotropic = turbulenceModelCoeffs_.lookupOrDefault<Switch>
+            (
+                "anisotropic",
+                false
+            );
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/incompressible/RNGkEpsilon/RNGkEpsilon.C
index 05fd38e9f0..cfe2d3d6b2 100644
--- a/src/turbulenceModels/incompressible/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/turbulenceModels/incompressible/RNGkEpsilon/RNGkEpsilon.C
@@ -51,13 +51,13 @@ RNGkEpsilon::RNGkEpsilon
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    alphak(turbulenceModelCoeffs_.lookup("alphak")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
-    eta0(turbulenceModelCoeffs_.lookup("eta0")),
-    beta(turbulenceModelCoeffs_.lookup("beta")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.0845)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.42)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.68)),
+    alphak(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.39)),
+    alphaEps(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 1.39)),
+    eta0(turbulenceModelCoeffs_.lookupOrDefault<scalar>("eta0", 4.38)),
+    beta(turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta", 0.012)),
 
     k_
     (
@@ -149,13 +149,17 @@ bool RNGkEpsilon::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("alphak") >> alphak;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
-        turbulenceModelCoeffs_.lookup("eta0") >> eta0;
-        turbulenceModelCoeffs_.lookup("beta") >> beta;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.0845);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.42);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.68);
+        alphak = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.39);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                1.39
+            );
+        eta0 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("eta0", 4.38);
+        beta = turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta", 0.012);
 
         return true;
     }
@@ -182,7 +186,7 @@ void RNGkEpsilon::correct()
     volScalarField G = nut_*S2;
 
     volScalarField eta = sqrt(S2)*k_/epsilon_;
-    volScalarField R = 
+    volScalarField R =
         ((eta*(scalar(1) - eta/eta0))/(scalar(1) + beta*eta*sqr(eta)));
 
 #   include "wallFunctionsI.H"
diff --git a/src/turbulenceModels/incompressible/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/incompressible/SpalartAllmaras/SpalartAllmaras.C
index 3523ba2a48..d2130de943 100644
--- a/src/turbulenceModels/incompressible/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/SpalartAllmaras/SpalartAllmaras.C
@@ -112,15 +112,15 @@ SpalartAllmaras::SpalartAllmaras
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    alphaNut(turbulenceModelCoeffs_.lookup("alphaNut")),
+    alphaNut(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaNut", 1.5)),
 
-    Cb1(turbulenceModelCoeffs_.lookup("Cb1")),
-    Cb2(turbulenceModelCoeffs_.lookup("Cb2")),
+    Cb1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cb1", 0.1355)),
+    Cb2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cb2", 0.622)),
     Cw1(Cb1/sqr(kappa_) + alphaNut*(1.0 + Cb2)),
-    Cw2(turbulenceModelCoeffs_.lookup("Cw2")),
-    Cw3(turbulenceModelCoeffs_.lookup("Cw3")),
-    Cv1(turbulenceModelCoeffs_.lookup("Cv1")),
-    Cv2(turbulenceModelCoeffs_.lookup("Cv2")),
+    Cw2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cw2", 0.3)),
+    Cw3(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cw3", 2.0)),
+    Cv1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cv1", 7.1)),
+    Cv2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cv2", 5.0)),
 
     nuTilda_
     (
@@ -257,15 +257,19 @@ bool SpalartAllmaras::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("alphaNut") >> alphaNut;
+        alphaNut = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaNut",
+                1.5
+            );
 
-        turbulenceModelCoeffs_.lookup("Cb1") >> Cb1;
-        turbulenceModelCoeffs_.lookup("Cb2") >> Cb2;
+        Cb1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cb1", 0.1355);
+        Cb2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cb2", 0.622);
         Cw1 = Cb1/sqr(kappa_) + alphaNut*(1.0 + Cb2);
-        turbulenceModelCoeffs_.lookup("Cw2") >> Cw2;
-        turbulenceModelCoeffs_.lookup("Cw3") >> Cw3;
-        turbulenceModelCoeffs_.lookup("Cv1") >> Cv1;
-        turbulenceModelCoeffs_.lookup("Cv2") >> Cv2;
+        Cw2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cw2", 0.3);
+        Cw3 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cw3", 2.0);
+        Cv1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cv1", 7.1);
+        Cv2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cv2", 5.0);
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/kEpsilon/kEpsilon.C b/src/turbulenceModels/incompressible/kEpsilon/kEpsilon.C
index 901887c54a..d18d1bf470 100644
--- a/src/turbulenceModels/incompressible/kEpsilon/kEpsilon.C
+++ b/src/turbulenceModels/incompressible/kEpsilon/kEpsilon.C
@@ -52,10 +52,13 @@ kEpsilon::kEpsilon
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    C1(turbulenceModelCoeffs_.lookup("C1")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    C1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92)),
+    alphaEps
+    (
+         turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923)
+    ),
 
     k_
     (
@@ -146,10 +149,11 @@ bool kEpsilon::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("C1") >> C1;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        C1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C1", 1.44);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.92);
+        alphaEps =
+            turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.76923);
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/incompressible/kOmegaSST/kOmegaSST.C
index 70c5e70ebe..0057a1a1e5 100644
--- a/src/turbulenceModels/incompressible/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/incompressible/kOmegaSST/kOmegaSST.C
@@ -94,22 +94,31 @@ kOmegaSST::kOmegaSST
 :
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
-    alphaK1(turbulenceModelCoeffs_.lookup("alphaK1")),
-    alphaK2(turbulenceModelCoeffs_.lookup("alphaK2")),
-    alphaOmega1(turbulenceModelCoeffs_.lookup("alphaOmega1")),
-    alphaOmega2(turbulenceModelCoeffs_.lookup("alphaOmega2")),
-    gamma1(turbulenceModelCoeffs_.lookup("gamma1")),
-    gamma2(turbulenceModelCoeffs_.lookup("gamma2")),
-    beta1(turbulenceModelCoeffs_.lookup("beta1")),
-    beta2(turbulenceModelCoeffs_.lookup("beta2")),
-    betaStar(turbulenceModelCoeffs_.lookup("betaStar")),
-    a1(turbulenceModelCoeffs_.lookup("a1")),
-    c1(turbulenceModelCoeffs_.lookup("c1")),
+    alphaK1
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaK1", 0.85034)
+    ),
+    alphaK2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaK2", 1.0)),
+    alphaOmega1
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaOmega1", 0.5)
+    ),
+    alphaOmega2
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaOmega2", 0.85616)
+    ),
+    gamma1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("gamma1", 0.5532)),
+    gamma2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("gamma2", 0.4403)),
+    beta1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta1", 0.075)),
+    beta2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta2", 0.0828)),
+    betaStar(turbulenceModelCoeffs_.lookupOrDefault<scalar>("betaStar", 0.09)),
+    a1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("a1", 0.31)),
+    c1(turbulenceModelCoeffs_.lookupOrDefault<scalar>("c1", 10.0)),
 
     omega0_("omega0", dimless/dimTime, SMALL),
     omegaSmall_("omegaSmall", dimless/dimTime, SMALL),
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
 
     y_(mesh_),
 
@@ -202,18 +211,42 @@ bool kOmegaSST::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("alphaK1") >> alphaK1;
-        turbulenceModelCoeffs_.lookup("alphaK2") >> alphaK2;
-        turbulenceModelCoeffs_.lookup("alphaOmega1") >> alphaOmega1;
-        turbulenceModelCoeffs_.lookup("alphaOmega2") >> alphaOmega2;
-        turbulenceModelCoeffs_.lookup("gamma1") >> gamma1;
-        turbulenceModelCoeffs_.lookup("gamma2") >> gamma2;
-        turbulenceModelCoeffs_.lookup("beta1") >> beta1;
-        turbulenceModelCoeffs_.lookup("beta2") >> beta2;
-        turbulenceModelCoeffs_.lookup("betaStar") >> betaStar;
-        turbulenceModelCoeffs_.lookup("a1") >> a1;
-        turbulenceModelCoeffs_.lookup("c1") >> c1;
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
+        alphaK1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaK1",
+                0.85034
+            );
+        alphaK2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaK2",
+                1.0
+            );
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaOmega1", 0.5);
+        alphaOmega2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaOmega2",
+                0.85616
+            );
+        gamma1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "gamma1",
+                0.5532
+            );
+        gamma2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "gamma2",
+                0.4403
+            );
+        beta1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta1", 0.075);
+        beta2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("beta2", 0.0828);
+        betaStar = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "betaStar",
+                0.09
+            );
+        a1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("a1", 0.31);
+        c1 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("c1", 10.0);
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
 
         return true;
     }
diff --git a/src/turbulenceModels/incompressible/realizableKE/realizableKE.C b/src/turbulenceModels/incompressible/realizableKE/realizableKE.C
index 45f82cb3a7..9bfe0ed084 100644
--- a/src/turbulenceModels/incompressible/realizableKE/realizableKE.C
+++ b/src/turbulenceModels/incompressible/realizableKE/realizableKE.C
@@ -93,11 +93,14 @@ realizableKE::realizableKE
     turbulenceModel(typeName, U, phi, lamTransportModel),
 
 
-    Cmu(turbulenceModelCoeffs_.lookup("Cmu")),
-    A0(turbulenceModelCoeffs_.lookup("A0")),
-    C2(turbulenceModelCoeffs_.lookup("C2")),
-    alphak(turbulenceModelCoeffs_.lookup("alphak")),
-    alphaEps(turbulenceModelCoeffs_.lookup("alphaEps")),
+    Cmu(turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09)),
+    A0(turbulenceModelCoeffs_.lookupOrDefault<scalar>("A0", 4.0)),
+    C2(turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.9)),
+    alphak(turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0)),
+    alphaEps
+    (
+        turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphaEps", 0.833333)
+    ),
 
     k_
     (
@@ -190,11 +193,15 @@ bool realizableKE::read()
 {
     if (turbulenceModel::read())
     {
-        turbulenceModelCoeffs_.lookup("Cmu") >> Cmu;
-        turbulenceModelCoeffs_.lookup("A0") >> A0;
-        turbulenceModelCoeffs_.lookup("C2") >> C2;
-        turbulenceModelCoeffs_.lookup("alphak") >> alphak;
-        turbulenceModelCoeffs_.lookup("alphaEps") >> alphaEps;
+        Cmu = turbulenceModelCoeffs_.lookupOrDefault<scalar>("Cmu", 0.09);
+        A0 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("A0", 4.0);
+        C2 = turbulenceModelCoeffs_.lookupOrDefault<scalar>("C2", 1.9);
+        alphak = turbulenceModelCoeffs_.lookupOrDefault<scalar>("alphak", 1.0);
+        alphaEps = turbulenceModelCoeffs_.lookupOrDefault<scalar>
+            (
+                "alphaEps",
+                0.833333
+            );
 
         return true;
     }