diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C
index ed4acb4da05ae4f21a0e521e5f029459c0ec5f2b..a23f1d1682f7cd8ed2c0736576d0f5dd9d38be57 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C
@@ -40,7 +40,7 @@ namespace RASModels
 {
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 scalar mutRoughWallFunctionFvPatchScalarField::fnRough
 (
@@ -65,6 +65,60 @@ scalar mutRoughWallFunctionFvPatchScalarField::fnRough
 }
 
 
+tmp<scalarField> mutRoughWallFunctionFvPatchScalarField::calcMut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
+    const tmp<volScalarField> tk = rasModel.k();
+    const volScalarField& k = tk();
+    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
+
+    const scalar Cmu25 = pow(Cmu_, 0.25);
+
+    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
+    scalarField& mutw = tmutw();
+
+    forAll(mutw, faceI)
+    {
+        label faceCellI = patch().faceCells()[faceI];
+
+        scalar uStar = Cmu25*sqrt(k[faceCellI]);
+
+        scalar yPlus = uStar*y[faceI]/(muw[faceI]/rhow[faceI]);
+
+        scalar KsPlus = uStar*Ks_[faceI]/(muw[faceI]/rhow[faceI]);
+
+        scalar Edash = E_;
+        scalar yPlusLamNew = yPlusLam_;
+        if (KsPlus > 2.25)
+        {
+            Edash /= fnRough(KsPlus, Cs_[faceI]);
+            yPlusLamNew = rasModel.yPlusLam(kappa_, Edash);
+        }
+
+        if (debug)
+        {
+            Info<< "yPlus = " << yPlus
+                << ", KsPlus = " << KsPlus
+                << ", Edash = " << Edash
+                << ", yPlusLam = " << yPlusLam_
+                << endl;
+        }
+
+        if (yPlus > yPlusLamNew)
+        {
+            mutw[faceI] =
+                muw[faceI]*(yPlus*kappa_/log(max(Edash*yPlus, 1+1e-4)) - 1);
+        }
+    }
+
+    return tmutw;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
@@ -158,61 +212,6 @@ void mutRoughWallFunctionFvPatchScalarField::rmap
 }
 
 
-tmp<scalarField> mutRoughWallFunctionFvPatchScalarField::calcMut() const
-{
-    const label patchI = patch().index();
-
-    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
-    const scalarField& y = rasModel.y()[patchI];
-    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
-    const tmp<volScalarField> tk = rasModel.k();
-    const volScalarField& k = tk();
-    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
-
-    const scalar Cmu25 = pow(Cmu_, 0.25);
-
-    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
-    scalarField& mutw = tmutw();
-
-    forAll(mutw, faceI)
-    {
-        label faceCellI = patch().faceCells()[faceI];
-
-        scalar uStar = Cmu25*sqrt(k[faceCellI]);
-
-        scalar yPlus = uStar*y[faceI]/(muw[faceI]/rhow[faceI]);
-
-        scalar KsPlus = uStar*Ks_[faceI]/(muw[faceI]/rhow[faceI]);
-
-        scalar Edash = E_;
-        scalar yPlusLamNew = yPlusLam;
-        if (KsPlus > 2.25)
-        {
-            Edash /= fnRough(KsPlus, Cs_[faceI]);
-            yPlusLamNew = rasModel.yPlusLam(kappa_, Edash);
-        }
-
-        if (debug)
-        {
-            Info<< "yPlus = " << yPlus
-                << ", KsPlus = " << KsPlus
-                << ", Edash = " << Edash
-                << ", yPlusLam = " << yPlusLam
-                << endl;
-        }
-
-        if (yPlus > yPlusLamNew)
-        {
-            mutw[faceI] =
-                muw[faceI]*(yPlus*kappa_/log(max(Edash*yPlus, 1+1e-4)) - 1);
-        }
-    }
-
-    return tmutw;
-}
-
-
 void mutRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchField<scalar>::write(os);
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
index c1876610959d0d2b96cd9fa1036ba8ce836ea133..3a1d0b7ab01a1a8ed62c6bc9f11a532150f5b2d2 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
@@ -39,6 +39,8 @@ namespace compressible
 namespace RASModels
 {
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
 tmp<scalarField>
 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
 (
@@ -48,7 +50,6 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
     const label patchI = patch().index();
 
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
     const scalarField& y = rasModel.y()[patchI];
     const scalarField& muw = rasModel.mu().boundaryField()[patchI];
     const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI];
@@ -74,7 +75,7 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
                 const scalar Re = rho[facei]*magUpara*y[facei]/muw[facei];
                 const scalar kappaRe = kappa_*Re;
 
-                scalar yp = yPlusLam;
+                scalar yp = yPlusLam_;
                 const scalar ryPlusLam = 1.0/yp;
 
                 int iter = 0;
@@ -146,7 +147,7 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
             const scalar Re = rho[facei]*magUpara*y[facei]/muw[facei];
             const scalar kappaRe = kappa_*Re;
 
-            scalar yp = yPlusLam;
+            scalar yp = yPlusLam_;
             const scalar ryPlusLam = 1.0/yp;
 
             int iter = 0;
@@ -171,6 +172,38 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
 }
 
 
+tmp<scalarField>
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
+    const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI];
+
+    scalarField magUp = mag(Uw.patchInternalField() - Uw);
+
+    tmp<scalarField> tyPlus = calcYPlus(magUp);
+    scalarField& yPlus = tyPlus();
+
+    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
+    scalarField& mutw = tmutw();
+
+    forAll(yPlus, facei)
+    {
+        if (yPlus[facei] > yPlusLam_)
+        {
+            const scalar Re = rho[facei]*magUp[facei]*y[facei]/muw[facei];
+            mutw[facei] = muw[facei]*(sqr(yPlus[facei])/Re - 1);
+        }
+    }
+
+    return tmutw;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
@@ -247,39 +280,6 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-tmp<scalarField>
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
-{
-    const label patchI = patch().index();
-
-    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
-    const scalarField& y = rasModel.y()[patchI];
-    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
-    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
-    const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI];
-
-    scalarField magUp = mag(Uw.patchInternalField() - Uw);
-
-    tmp<scalarField> tyPlus = calcYPlus(magUp);
-    scalarField& yPlus = tyPlus();
-
-    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
-    scalarField& mutw = tmutw();
-
-    forAll(yPlus, facei)
-    {
-        if (yPlus[facei] > yPlusLam)
-        {
-            const scalar Re = rho[facei]*magUp[facei]*y[facei]/muw[facei];
-            mutw[facei] = muw[facei]*(sqr(yPlus[facei])/Re - 1);
-        }
-    }
-
-    return tmutw;
-}
-
-
 tmp<scalarField>
 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus() const
 {
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
index 052b22664b42d1f92358729e2b8cbde466c1ae41..c2a49ac7510ef48bd1b38b8e95a693328abc1ba0 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
@@ -50,7 +50,6 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
     const label patchI = patch().index();
 
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
     const scalarField& y = rasModel.y()[patchI];
     const fvPatchScalarField& rhow = rasModel.rho().boundaryField()[patchI];
     const fvPatchScalarField& muw = rasModel.mu().boundaryField()[patchI];
@@ -62,7 +61,7 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
     {
         scalar kappaRe = kappa_*magUp[faceI]*y[faceI]/(muw[faceI]/rhow[faceI]);
 
-        scalar yp = yPlusLam;
+        scalar yp = yPlusLam_;
         scalar ryPlusLam = 1.0/yp;
 
         int iter = 0;
@@ -82,6 +81,35 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
 }
 
 
+tmp<scalarField>
+mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcMut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField magUp = mag(Uw.patchInternalField() - Uw);
+    const fvPatchScalarField& muw = rasModel.mu().boundaryField()[patchI];
+
+    tmp<scalarField> tyPlus = calcYPlus(magUp);
+    scalarField& yPlus = tyPlus();
+
+    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
+    scalarField& mutw = tmutw();
+
+    forAll(yPlus, faceI)
+    {
+        if (yPlus[faceI] > yPlusLam_)
+        {
+            mutw[faceI] =
+                muw[faceI]*(yPlus[faceI]*kappa_/log(E_*yPlus[faceI]) - 1.0);
+        }
+    }
+
+    return tmutw;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::
@@ -143,36 +171,6 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-tmp<scalarField>
-mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcMut() const
-{
-    const label patchI = patch().index();
-
-    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
-    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
-    const scalarField magUp = mag(Uw.patchInternalField() - Uw);
-    const fvPatchScalarField& muw = rasModel.mu().boundaryField()[patchI];
-
-    tmp<scalarField> tyPlus = calcYPlus(magUp);
-    scalarField& yPlus = tyPlus();
-
-    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
-    scalarField& mutw = tmutw();
-
-    forAll(yPlus, faceI)
-    {
-        if (yPlus[faceI] > yPlusLam)
-        {
-            mutw[faceI] =
-                muw[faceI]*(yPlus[faceI]*kappa_/log(E_*yPlus[faceI]) - 1.0);
-        }
-    }
-
-    return tmutw;
-}
-
-
 tmp<scalarField>
 mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() const
 {
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasWallFunction/mutSpalartAllmarasWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasWallFunction/mutSpalartAllmarasWallFunctionFvPatchScalarField.C
index 23874146ea1b65f91386026e2892e6a83c9e0190..85e9c7e60a447e83a47ed6ab456a91634e3cc8bd 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasWallFunction/mutSpalartAllmarasWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasWallFunction/mutSpalartAllmarasWallFunctionFvPatchScalarField.C
@@ -105,6 +105,21 @@ tmp<scalarField> mutSpalartAllmarasWallFunctionFvPatchScalarField::calcUTau
 }
 
 
+tmp<scalarField>
+mutSpalartAllmarasWallFunctionFvPatchScalarField::calcMut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField magGradU = mag(Uw.snGrad());
+    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
+    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
+
+    return max(0.0, rhow*sqr(calcUTau(magGradU))/magGradU - muw);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 mutSpalartAllmarasWallFunctionFvPatchScalarField::
@@ -166,21 +181,6 @@ mutSpalartAllmarasWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-tmp<scalarField>
-mutSpalartAllmarasWallFunctionFvPatchScalarField::calcMut() const
-{
-    const label patchI = patch().index();
-
-    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
-    const scalarField magGradU = mag(Uw.snGrad());
-    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
-    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
-
-    return max(0.0, rhow*sqr(calcUTau(magGradU))/magGradU - muw);
-}
-
-
 tmp<scalarField>
 mutSpalartAllmarasWallFunctionFvPatchScalarField::yPlus() const
 {
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.C
index dcd52cc4e1b5b54ba500334d068f7f815ccc9165..10813c9a055a9dc88ba7cc1cbdb46f0b5c6f000e 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.C
@@ -40,7 +40,7 @@ namespace compressible
 namespace RASModels
 {
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 void mutWallFunctionFvPatchScalarField::checkType()
 {
@@ -56,6 +56,63 @@ void mutWallFunctionFvPatchScalarField::checkType()
 }
 
 
+scalar mutWallFunctionFvPatchScalarField::calcYPlusLam
+(
+    const scalar kappa,
+    const scalar E
+) const
+{
+    scalar ypl = 11.0;
+
+    for (int i=0; i<10; i++)
+    {
+        ypl = log(E*ypl)/kappa;
+    }
+
+    return ypl;
+}
+
+
+tmp<scalarField> mutWallFunctionFvPatchScalarField::calcMut() const
+{
+    const label patchI = patch().index();
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
+    const tmp<volScalarField> tk = rasModel.k();
+    const volScalarField& k = tk();
+    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
+
+    const scalar Cmu25 = pow(Cmu_, 0.25);
+
+    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
+    scalarField& mutw = tmutw();
+
+    forAll(mutw, faceI)
+    {
+        label faceCellI = patch().faceCells()[faceI];
+
+        scalar yPlus =
+            Cmu25*y[faceI]*sqrt(k[faceCellI])/(muw[faceI]/rhow[faceI]);
+
+        if (yPlus > yPlusLam_)
+        {
+            mutw[faceI] = muw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1);
+        }
+    }
+
+    return tmutw;
+}
+
+
+void mutWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
+{
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
@@ -67,7 +124,8 @@ mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
     fixedValueFvPatchScalarField(p, iF),
     Cmu_(0.09),
     kappa_(0.41),
-    E_(9.8)
+    E_(9.8),
+    yPlusLam_(calcYPlusLam(kappa_, E_))
 {}
 
 
@@ -82,7 +140,8 @@ mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
     Cmu_(ptf.Cmu_),
     kappa_(ptf.kappa_),
-    E_(ptf.E_)
+    E_(ptf.E_),
+    yPlusLam_(ptf.yPlusLam_)
 {}
 
 
@@ -96,7 +155,8 @@ mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
     fixedValueFvPatchScalarField(p, iF, dict),
     Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
-    E_(dict.lookupOrDefault<scalar>("E", 9.8))
+    E_(dict.lookupOrDefault<scalar>("E", 9.8)),
+    yPlusLam_(calcYPlusLam(kappa_, E_))
 {}
 
 
@@ -108,7 +168,8 @@ mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
     fixedValueFvPatchScalarField(wfpsf),
     Cmu_(wfpsf.Cmu_),
     kappa_(wfpsf.kappa_),
-    E_(wfpsf.E_)
+    E_(wfpsf.E_),
+    yPlusLam_(wfpsf.yPlusLam_)
 {}
 
 
@@ -121,7 +182,8 @@ mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
     fixedValueFvPatchScalarField(wfpsf, iF),
     Cmu_(wfpsf.Cmu_),
     kappa_(wfpsf.kappa_),
-    E_(wfpsf.E_)
+    E_(wfpsf.E_),
+    yPlusLam_(wfpsf.yPlusLam_)
 {}
 
 
@@ -135,40 +197,6 @@ void mutWallFunctionFvPatchScalarField::updateCoeffs()
 }
 
 
-tmp<scalarField> mutWallFunctionFvPatchScalarField::calcMut() const
-{
-    const label patchI = patch().index();
-    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
-    const scalarField& y = rasModel.y()[patchI];
-    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
-    const tmp<volScalarField> tk = rasModel.k();
-    const volScalarField& k = tk();
-    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
-
-    const scalar Cmu25 = pow(Cmu_, 0.25);
-
-    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
-    scalarField& mutw = tmutw();
-
-    forAll(mutw, faceI)
-    {
-        label faceCellI = patch().faceCells()[faceI];
-
-        scalar yPlus =
-            Cmu25*y[faceI]*sqrt(k[faceCellI])
-           /(muw[faceI]/rhow[faceI]);
-
-        if (yPlus > yPlusLam)
-        {
-            mutw[faceI] = muw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1);
-        }
-    }
-
-    return tmutw;
-}
-
-
 tmp<scalarField> mutWallFunctionFvPatchScalarField::yPlus() const
 {
     const label patchI = patch().index();
@@ -194,14 +222,6 @@ void mutWallFunctionFvPatchScalarField::write(Ostream& os) const
 }
 
 
-void mutWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
-{
-    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
-    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
-    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makePatchTypeField(fvPatchScalarField, mutWallFunctionFvPatchScalarField);
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.H
index 3945c345e9452db2a50337ca096fe25c7e2a2f35..6134256e2311e8c57b07f9fdae75328d3f7c3d1b 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.H
@@ -70,12 +70,18 @@ protected:
         //- E coefficient
         scalar E_;
 
+        //- Y+ at the edge of the laminar sublayer
+        scalar yPlusLam_;
+
 
     // Protected member functions
 
         //- Check the type of the patch
         virtual void checkType();
 
+        //- Calculate the Y+ at the edge of the laminar sublayer
+        virtual scalar calcYPlusLam(const scalar kappa, const scalar E) const;
+
         //- Calculate the turbulence viscosity
         virtual tmp<scalarField> calcMut() const;
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C
index e26e652a1aa66bfe497c47b4d661181659fc71c4..9ef3051f82472b40d698372e249417553149680f 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C
@@ -39,7 +39,7 @@ namespace incompressible
 namespace RASModels
 {
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 scalar nutRoughWallFunctionFvPatchScalarField::fnRough
 (
@@ -64,6 +64,68 @@ scalar nutRoughWallFunctionFvPatchScalarField::fnRough
 }
 
 
+tmp<scalarField> nutRoughWallFunctionFvPatchScalarField::calcNut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+    const tmp<volScalarField> tk = rasModel.k();
+    const volScalarField& k = tk();
+    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
+
+    const scalar Cmu25 = pow(Cmu_, 0.25);
+
+    tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
+    scalarField& nutw = tnutw();
+
+    forAll(nutw, faceI)
+    {
+        label faceCellI = patch().faceCells()[faceI];
+
+        scalar uStar = Cmu25*sqrt(k[faceCellI]);
+        scalar yPlus = uStar*y[faceI]/nuw[faceI];
+        scalar KsPlus = uStar*Ks_[faceI]/nuw[faceI];
+
+        scalar Edash = E_;
+
+        if (KsPlus > 2.25)
+        {
+            Edash /= fnRough(KsPlus, Cs_[faceI]);
+        }
+
+        if (yPlus > yPlusLam_)
+        {
+            scalar limitingNutw = max(nutw[faceI], nuw[faceI]);
+
+            // To avoid oscillations limit the change in the wall viscosity
+            // which is particularly important if it temporarily becomes zero
+            nutw[faceI] =
+                max
+                (
+                    min
+                    (
+                        nuw[faceI]
+                       *(yPlus*kappa_/log(max(Edash*yPlus, 1+1e-4)) - 1),
+                        2*limitingNutw
+                    ), 0.5*limitingNutw
+                );
+        }
+
+        if (debug)
+        {
+            Info<< "yPlus = " << yPlus
+                << ", KsPlus = " << KsPlus
+                << ", Edash = " << Edash
+                << ", nutw = " << nutw[faceI]
+                << endl;
+        }
+    }
+
+    return tnutw;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 nutRoughWallFunctionFvPatchScalarField::nutRoughWallFunctionFvPatchScalarField
@@ -157,69 +219,6 @@ void nutRoughWallFunctionFvPatchScalarField::rmap
 }
 
 
-tmp<scalarField> nutRoughWallFunctionFvPatchScalarField::calcNut() const
-{
-    const label patchI = patch().index();
-
-    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
-    const scalarField& y = rasModel.y()[patchI];
-    const tmp<volScalarField> tk = rasModel.k();
-    const volScalarField& k = tk();
-    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
-
-    const scalar Cmu25 = pow(Cmu_, 0.25);
-
-    tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
-    scalarField& nutw = tnutw();
-
-    forAll(nutw, faceI)
-    {
-        label faceCellI = patch().faceCells()[faceI];
-
-        scalar uStar = Cmu25*sqrt(k[faceCellI]);
-        scalar yPlus = uStar*y[faceI]/nuw[faceI];
-        scalar KsPlus = uStar*Ks_[faceI]/nuw[faceI];
-
-        scalar Edash = E_;
-
-        if (KsPlus > 2.25)
-        {
-            Edash /= fnRough(KsPlus, Cs_[faceI]);
-        }
-
-        if (yPlus > yPlusLam)
-        {
-            scalar limitingNutw = max(nutw[faceI], nuw[faceI]);
-
-            // To avoid oscillations limit the change in the wall viscosity
-            // which is particularly important if it temporarily becomes zero
-            nutw[faceI] =
-                max
-                (
-                    min
-                    (
-                        nuw[faceI]
-                       *(yPlus*kappa_/log(max(Edash*yPlus, 1+1e-4)) - 1),
-                        2*limitingNutw
-                    ), 0.5*limitingNutw
-                );
-        }
-
-        if (debug)
-        {
-            Info<< "yPlus = " << yPlus
-                << ", KsPlus = " << KsPlus
-                << ", Edash = " << Edash
-                << ", nutw = " << nutw[faceI]
-                << endl;
-        }
-    }
-
-    return tnutw;
-}
-
-
 void nutRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchField<scalar>::write(os);
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
index e2de2dae71e66b405f17e27aba5404ae6f076a80..760559f6798a2863069b2776ffa925ec0a91a09f 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
@@ -39,7 +39,39 @@ namespace incompressible
 namespace RASModels
 {
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+tmp<scalarField>
+nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcNut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
+
+    // The flow velocity at the adjacent cell centre
+    const scalarField magUp = mag(Uw.patchInternalField() - Uw);
+
+    tmp<scalarField> tyPlus = calcYPlus(magUp);
+    scalarField& yPlus = tyPlus();
+
+    tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
+    scalarField& nutw = tnutw();
+
+    forAll(yPlus, facei)
+    {
+        if (yPlus[facei] > yPlusLam_)
+        {
+            const scalar Re = magUp[facei]*y[facei]/nuw[facei];
+            nutw[facei] = nuw[facei]*(sqr(yPlus[facei])/Re - 1);
+        }
+    }
+
+    return tnutw;
+}
+
 
 tmp<scalarField>
 nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
@@ -50,7 +82,6 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
     const label patchI = patch().index();
 
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
     const scalarField& y = rasModel.y()[patchI];
     const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
 
@@ -75,7 +106,7 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
                 const scalar Re = magUpara*y[facei]/nuw[facei];
                 const scalar kappaRe = kappa_*Re;
 
-                scalar yp = yPlusLam;
+                scalar yp = yPlusLam_;
                 const scalar ryPlusLam = 1.0/yp;
 
                 int iter = 0;
@@ -133,7 +164,7 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
                  && yp > VSMALL
                 );
 
-                yPlus[facei] = yp;
+                yPlus[facei] = max(0.0, yp);
             }
         }
     }
@@ -146,7 +177,7 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
             const scalar Re = magUpara*y[facei]/nuw[facei];
             const scalar kappaRe = kappa_*Re;
 
-            scalar yp = yPlusLam;
+            scalar yp = yPlusLam_;
             const scalar ryPlusLam = 1.0/yp;
 
             int iter = 0;
@@ -159,7 +190,7 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
 
             } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.0001 && ++iter < 10);
 
-            yPlus[facei] = yp;
+            yPlus[facei] = max(0.0, yp);
         }
     }
 
@@ -243,39 +274,6 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-tmp<scalarField>
-nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcNut() const
-{
-    const label patchI = patch().index();
-
-    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
-    const scalarField& y = rasModel.y()[patchI];
-    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
-    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
-
-    // The flow velocity at the adjacent cell centre
-    const scalarField magUp = mag(Uw.patchInternalField() - Uw);
-
-    tmp<scalarField> tyPlus = calcYPlus(magUp);
-    scalarField& yPlus = tyPlus();
-
-    tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
-    scalarField& nutw = tnutw();
-
-    forAll(yPlus, facei)
-    {
-        if (yPlus[facei] > yPlusLam)
-        {
-            const scalar Re = magUp[facei]*y[facei]/nuw[facei];
-            nutw[facei] = nuw[facei]*(sqr(yPlus[facei])/Re - 1);
-        }
-    }
-
-    return tnutw;
-}
-
-
 tmp<scalarField>
 nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus() const
 {
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
index 7bea1742fb8e66a6c552e30330a77e6fe48e824e..e3c3e556711e8c174e2694f4b81f5ae38cee8d31 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
@@ -41,6 +41,35 @@ namespace RASModels
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+tmp<scalarField>
+nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcNut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField magUp = mag(Uw.patchInternalField() - Uw);
+    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
+
+    tmp<scalarField> tyPlus = calcYPlus(magUp);
+    scalarField& yPlus = tyPlus();
+
+    tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
+    scalarField& nutw = tnutw();
+
+    forAll(yPlus, facei)
+    {
+        if (yPlus[facei] > yPlusLam_)
+        {
+            nutw[facei] =
+                nuw[facei]*(yPlus[facei]*kappa_/log(E_*yPlus[facei]) - 1.0);
+        }
+    }
+
+    return tnutw;
+}
+
+
 tmp<scalarField>
 nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
 (
@@ -50,7 +79,6 @@ nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
     const label patchI = patch().index();
 
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
     const scalarField& y = rasModel.y()[patchI];
     const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
 
@@ -61,7 +89,7 @@ nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
     {
         scalar kappaRe = kappa_*magUp[facei]*y[facei]/nuw[facei];
 
-        scalar yp = yPlusLam;
+        scalar yp = yPlusLam_;
         scalar ryPlusLam = 1.0/yp;
 
         int iter = 0;
@@ -142,36 +170,6 @@ nutSpalartAllmarasStandardWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-tmp<scalarField>
-nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcNut() const
-{
-    const label patchI = patch().index();
-
-    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
-    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
-    const scalarField magUp = mag(Uw.patchInternalField() - Uw);
-    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
-
-    tmp<scalarField> tyPlus = calcYPlus(magUp);
-    scalarField& yPlus = tyPlus();
-
-    tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
-    scalarField& nutw = tnutw();
-
-    forAll(yPlus, facei)
-    {
-        if (yPlus[facei] > yPlusLam)
-        {
-            nutw[facei] =
-                nuw[facei]*(yPlus[facei]*kappa_/log(E_*yPlus[facei]) - 1.0);
-        }
-    }
-
-    return tnutw;
-}
-
-
 tmp<scalarField>
 nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() const
 {
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasWallFunction/nutSpalartAllmarasWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasWallFunction/nutSpalartAllmarasWallFunctionFvPatchScalarField.C
index 99187abb64279d758ebf32bd09db3303ed82d89a..a2d4b746a827b090735968c9e3332227ee11dcf2 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasWallFunction/nutSpalartAllmarasWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasWallFunction/nutSpalartAllmarasWallFunctionFvPatchScalarField.C
@@ -39,7 +39,21 @@ namespace incompressible
 namespace RASModels
 {
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+tmp<scalarField>
+nutSpalartAllmarasWallFunctionFvPatchScalarField::calcNut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField magGradU = mag(Uw.snGrad());
+    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
+
+    return max(0.0, sqr(calcUTau(magGradU))/magGradU - nuw);
+}
+
 
 tmp<scalarField> nutSpalartAllmarasWallFunctionFvPatchScalarField::calcUTau
 (
@@ -159,20 +173,6 @@ nutSpalartAllmarasWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-tmp<scalarField>
-nutSpalartAllmarasWallFunctionFvPatchScalarField::calcNut() const
-{
-    const label patchI = patch().index();
-
-    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
-    const scalarField magGradU = mag(Uw.snGrad());
-    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
-
-    return max(0.0, sqr(calcUTau(magGradU))/magGradU - nuw);
-}
-
-
 tmp<scalarField>
 nutSpalartAllmarasWallFunctionFvPatchScalarField::yPlus() const
 {
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.C
index 0fdefb2b51431dbf3116acf394e0e8c94db94b1d..1417f3b64c7c43131f27d9bf1c893a8b3ea246e7 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.C
@@ -56,6 +56,62 @@ void nutWallFunctionFvPatchScalarField::checkType()
 }
 
 
+scalar nutWallFunctionFvPatchScalarField::calcYPlusLam
+(
+    const scalar kappa,
+    const scalar E
+) const
+{
+    scalar ypl = 11.0;
+
+    for (int i=0; i<10; i++)
+    {
+        ypl = log(E*ypl)/kappa;
+    }
+
+    return ypl;
+}
+
+
+tmp<scalarField> nutWallFunctionFvPatchScalarField::calcNut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+    const tmp<volScalarField> tk = rasModel.k();
+    const volScalarField& k = tk();
+    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
+
+    const scalar Cmu25 = pow(Cmu_, 0.25);
+
+    tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
+    scalarField& nutw = tnutw();
+
+    forAll(nutw, faceI)
+    {
+        label faceCellI = patch().faceCells()[faceI];
+
+        scalar yPlus = Cmu25*y[faceI]*sqrt(k[faceCellI])/nuw[faceI];
+
+        if (yPlus > yPlusLam_)
+        {
+            nutw[faceI] = nuw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1.0);
+        }
+    }
+
+    return tnutw;
+}
+
+
+void nutWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
+{
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField
@@ -67,7 +123,8 @@ nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField
     fixedValueFvPatchScalarField(p, iF),
     Cmu_(0.09),
     kappa_(0.41),
-    E_(9.8)
+    E_(9.8),
+    yPlusLam_(calcYPlusLam(kappa_, E_))
 {
     checkType();
 }
@@ -84,7 +141,8 @@ nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField
     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
     Cmu_(ptf.Cmu_),
     kappa_(ptf.kappa_),
-    E_(ptf.E_)
+    E_(ptf.E_),
+    yPlusLam_(ptf.yPlusLam_)
 {
     checkType();
 }
@@ -100,7 +158,8 @@ nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField
     fixedValueFvPatchScalarField(p, iF, dict),
     Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
-    E_(dict.lookupOrDefault<scalar>("E", 9.8))
+    E_(dict.lookupOrDefault<scalar>("E", 9.8)),
+    yPlusLam_(calcYPlusLam(kappa_, E_))
 {
     checkType();
 }
@@ -114,7 +173,8 @@ nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField
     fixedValueFvPatchScalarField(wfpsf),
     Cmu_(wfpsf.Cmu_),
     kappa_(wfpsf.kappa_),
-    E_(wfpsf.E_)
+    E_(wfpsf.E_),
+    yPlusLam_(wfpsf.yPlusLam_)
 {
     checkType();
 }
@@ -129,7 +189,8 @@ nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField
     fixedValueFvPatchScalarField(wfpsf, iF),
     Cmu_(wfpsf.Cmu_),
     kappa_(wfpsf.kappa_),
-    E_(wfpsf.E_)
+    E_(wfpsf.E_),
+    yPlusLam_(wfpsf.yPlusLam_)
 {
     checkType();
 }
@@ -145,38 +206,6 @@ void nutWallFunctionFvPatchScalarField::updateCoeffs()
 }
 
 
-tmp<scalarField> nutWallFunctionFvPatchScalarField::calcNut() const
-{
-    const label patchI = patch().index();
-
-    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
-    const scalarField& y = rasModel.y()[patchI];
-    const tmp<volScalarField> tk = rasModel.k();
-    const volScalarField& k = tk();
-    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
-
-    const scalar Cmu25 = pow(Cmu_, 0.25);
-
-    tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
-    scalarField& nutw = tnutw();
-
-    forAll(nutw, faceI)
-    {
-        label faceCellI = patch().faceCells()[faceI];
-
-        scalar yPlus = Cmu25*y[faceI]*sqrt(k[faceCellI])/nuw[faceI];
-
-        if (yPlus > yPlusLam)
-        {
-            nutw[faceI] = nuw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1.0);
-        }
-    }
-
-    return tnutw;
-}
-
-
 tmp<scalarField> nutWallFunctionFvPatchScalarField::yPlus() const
 {
     const label patchI = patch().index();
@@ -201,14 +230,6 @@ void nutWallFunctionFvPatchScalarField::write(Ostream& os) const
 }
 
 
-void nutWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
-{
-    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
-    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
-    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makePatchTypeField(fvPatchScalarField, nutWallFunctionFvPatchScalarField);
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.H
index cb68c36f1f285a5ec4e9467d2078ca7350ce0bb4..490d14cd0108ea55e8c9b416724b354ee69240c7 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.H
@@ -70,12 +70,18 @@ protected:
         //- E coefficient
         scalar E_;
 
+        //- Y+ at the edge of the laminar sublayer
+        scalar yPlusLam_;
+
 
     // Protected member functions
 
         //- Check the type of the patch
         virtual void checkType();
 
+        //- Calculate the Y+ at the edge of the laminar sublayer
+        virtual scalar calcYPlusLam(const scalar kappa, const scalar E) const;
+
         //- Calculate the turbulence viscosity
         virtual tmp<scalarField> calcNut() const;