diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
index eaaafc99d4..5adde282fa 100644
--- a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
@@ -50,7 +50,7 @@ PDRkEpsilon::PDRkEpsilon
     const volScalarField& rho,
     const volVectorField& U,
     const surfaceScalarField& phi,
-    basicThermo& thermophysicalModel
+    const basicThermo& thermophysicalModel
 )
 :
     RASModel(typeName, rho, U, phi, thermophysicalModel),
diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H
index 0866a28a8e..7a8f7ae54d 100644
--- a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H
@@ -93,7 +93,7 @@ public:
             const volScalarField& rho,
             const volVectorField& U,
             const surfaceScalarField& phi,
-            basicThermo& thermophysicalModel
+            const basicThermo& thermophysicalModel
         );
 
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
index ee7f2943bd..b6bcfffa3f 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
@@ -129,9 +129,20 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
     {
         Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
             << "updateCoeffs() :"
+            << " patch:" << patch().name()
             << " walltemperature "
-            << " min:" << gMin(*this)
-            << " max:" << gMax(*this)
+            << " min:"
+            <<  returnReduce
+                (
+                    (this->size() > 0 ? min(*this) : VGREAT),
+                    minOp<scalar>()
+                )
+            << " max:"
+            <<  returnReduce
+                (
+                    (this->size() > 0 ? max(*this) : -VGREAT),
+                    maxOp<scalar>()
+                )
             << " avg:" << gAverage(*this)
             << endl;
     }
@@ -163,7 +174,9 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
         label nTotSize = returnReduce(this->size(), sumOp<label>());
 
         Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-            << "updateCoeffs() : Out of " << nTotSize
+            << "updateCoeffs() :"
+            << " patch:" << patch().name()
+            << " out of:" << nTotSize
             << " fixedBC:" << nFixed
             << " gradient:" << nTotSize-nFixed << endl;
     }
@@ -213,9 +226,22 @@ void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::evaluate
 
         if (debug)
         {
-            Info<< "Setting master and slave to wall temperature "
-                << " min:" << gMin(*this)
-                << " max:" << gMax(*this)
+            Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
+                << "updateCoeffs() :"
+                << " patch:" << patch().name()
+                << " setting master and slave to wall temperature "
+                << " min:"
+                <<  returnReduce
+                    (
+                        (this->size() > 0 ? min(*this) : VGREAT),
+                        minOp<scalar>()
+                    )
+                << " max:"
+                <<  returnReduce
+                    (
+                        (this->size() > 0 ? max(*this) : -VGREAT),
+                        maxOp<scalar>()
+                    )
                 << " avg:" << gAverage(*this)
                 << endl;
         }
diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index e3cb83ad38..8aac1f6b7e 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -488,17 +488,18 @@ labelList getNonRegionCells(const labelList& cellRegion, const label regionI)
 }
 
 
-// Get per region-region interface the sizes.
-// If sumParallel does merge.
-EdgeMap<label> getInterfaceSizes
+// Get per region-region interface the sizes. If sumParallel sums sizes.
+// Returns interfaces as straight list for looping in identical order.
+void getInterfaceSizes
 (
     const polyMesh& mesh,
     const labelList& cellRegion,
-    const bool sumParallel
+    const bool sumParallel,
+
+    edgeList& interfaces,
+    EdgeMap<label>& interfaceSizes
 )
 {
-    EdgeMap<label> interfaceSizes;
-
     forAll(mesh.faceNeighbour(), faceI)
     {
         label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
@@ -585,7 +586,12 @@ EdgeMap<label> getInterfaceSizes
         }
     }
 
-    return interfaceSizes;
+    // Make sure all processors have interfaces in same order
+    interfaces = interfaceSizes.toc();
+    if (sumParallel)
+    {
+        Pstream::scatter(interfaces);
+    }
 }
 
 
@@ -705,11 +711,7 @@ autoPtr<mapPolyMesh> createRegionMesh
 
         if (otherRegion != -1)
         {
-            edge interface
-            (
-                min(regionI, otherRegion),
-                max(regionI, otherRegion)
-            );
+            edge interface(regionI, otherRegion);
 
             // Find the patch.
             if (regionI < otherRegion)
@@ -848,6 +850,7 @@ void createAndWriteRegion
 
 
     const polyBoundaryMesh& newPatches = newMesh().boundaryMesh();
+    newPatches.checkParallelSync(true);
 
     // Delete empty patches
     // ~~~~~~~~~~~~~~~~~~~~
@@ -863,13 +866,12 @@ void createAndWriteRegion
     {
         const polyPatch& pp = newPatches[patchI];
 
-        if
-        (
-            !isA<processorPolyPatch>(pp)
-         && returnReduce(pp.size(), sumOp<label>()) > 0
-        )
+        if (!isA<processorPolyPatch>(pp))
         {
-            oldToNew[patchI] = newI++;
+            if (returnReduce(pp.size(), sumOp<label>()) > 0)
+            {
+                oldToNew[patchI] = newI++;
+            }
         }
     }
 
@@ -983,10 +985,15 @@ void createAndWriteRegion
 }
 
 
+// Create for every region-region interface with non-zero size two patches.
+// First one is for minimumregion to maximumregion.
+// Note that patches get created in same order on all processors (if parallel)
+// since looping over synchronised 'interfaces'.
 EdgeMap<label> addRegionPatches
 (
     fvMesh& mesh,
     const regionSplit& cellRegion,
+    const edgeList& interfaces,
     const EdgeMap<label>& interfaceSizes,
     const wordList& regionNames
 )
@@ -998,15 +1005,12 @@ EdgeMap<label> addRegionPatches
 
     EdgeMap<label> interfaceToPatch(cellRegion.nRegions());
 
-    // Keep start of added patches for later.
-    label minAddedPatchI = labelMax;
-
-    forAllConstIter(EdgeMap<label>, interfaceSizes, iter)
+    forAll(interfaces, interI)
     {
-        if (iter() > 0)
-        {
-            const edge& e = iter.key();
+        const edge& e = interfaces[interI];
 
+        if (interfaceSizes[e] > 0)
+        {
             label patchI = addPatch
             (
                 mesh,
@@ -1025,12 +1029,9 @@ EdgeMap<label> addRegionPatches
                 << " " << mesh.boundaryMesh()[patchI].name()
                 << endl;
 
-            interfaceToPatch.insert(iter.key(), patchI);
-
-            minAddedPatchI = min(minAddedPatchI, patchI);
+            interfaceToPatch.insert(e, patchI);
         }
     }
-    //Info<< "minAddedPatchI:" << minAddedPatchI << endl;
     return interfaceToPatch;
 }
 
@@ -1348,24 +1349,26 @@ int main(int argc, char *argv[])
 
     // Sizes of interface between regions. From pair of regions to number of
     // faces.
-    EdgeMap<label> interfaceSizes
+    edgeList interfaces;
+    EdgeMap<label> interfaceSizes;
+    getInterfaceSizes
     (
-        getInterfaceSizes
-        (
-            mesh,
-            cellRegion,
-            true       // sum in parallel?
-        )
+        mesh,
+        cellRegion,
+        true,      // sum in parallel?
+
+        interfaces,
+        interfaceSizes
     );
 
     Info<< "Region\tRegion\tFaces" << nl
         << "------\t------\t-----" << endl;
 
-    forAllConstIter(EdgeMap<label>, interfaceSizes, iter)
+    forAll(interfaces, interI)
     {
-        const edge& e = iter.key();
+        const edge& e = interfaces[interI];
 
-        Info<< e[0] << '\t' << e[1] << '\t' << iter() << nl;
+        Info<< e[0] << '\t' << e[1] << '\t' << interfaceSizes[e] << nl;
     }
     Info<< endl;
 
@@ -1511,6 +1514,7 @@ int main(int argc, char *argv[])
             (
                 mesh,
                 cellRegion,
+                interfaces,
                 interfaceSizes,
                 regionNames
             )
diff --git a/src/OSspecific/Unix/Unix.C b/src/OSspecific/Unix/Unix.C
index e4ffe76857..c3479810d4 100644
--- a/src/OSspecific/Unix/Unix.C
+++ b/src/OSspecific/Unix/Unix.C
@@ -91,7 +91,9 @@ Foam::string Foam::getEnv(const word& envName)
     }
     else
     {
-        return string::null;
+        // Return null-constructed string rather than string::null
+        // to avoid cyclic dependencies in the construction of globals
+        return string();
     }
 }
 
@@ -277,7 +279,9 @@ Foam::fileName Foam::findEtcFile(const fileName& name, bool mandatory)
         ::exit(1);
     }
 
-    return fileName::null;
+    // Return null-constructed fileName rather than fileName::null
+    // to avoid cyclic dependencies in the construction of globals
+    return fileName();
 }
 
 
diff --git a/src/OpenFOAM/global/global.Cver b/src/OpenFOAM/global/global.Cver
index 677c30ecc7..319cbfd1ad 100644
--- a/src/OpenFOAM/global/global.Cver
+++ b/src/OpenFOAM/global/global.Cver
@@ -38,11 +38,6 @@ Description
 const char* const Foam::FOAMversion = "VERSION_STRING";
 const char* const Foam::FOAMbuild = "BUILD_STRING";
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Static initializers for string::null, word::null and fileName::null
-
-#include "stringsGlobals.C"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Setup an error handler for the global new operator
 
diff --git a/src/OpenFOAM/include/OSspecific.H b/src/OpenFOAM/include/OSspecific.H
index 1cca34885d..4594f6eff8 100644
--- a/src/OpenFOAM/include/OSspecific.H
+++ b/src/OpenFOAM/include/OSspecific.H
@@ -62,7 +62,7 @@ pid_t pgid();
 bool env(const word&);
 
 //- Return environment variable of given name
-//  Return string::null if the environment is undefined
+//  Return string() if the environment is undefined
 string getEnv(const word& name);
 
 //- Set an environment variable
@@ -101,7 +101,7 @@ bool chDir(const fileName& dir);
 //  -# shipped settings:
 //    - $WM_PROJECT_DIR/etc/
 //
-//  @return the full path name or fileName::null if the name cannot be found
+//  @return the full path name or fileName() if the name cannot be found
 //  Optionally abort if the file cannot be found
 fileName findEtcFile(const fileName& name, bool mandatory=false);
 
diff --git a/src/OpenFOAM/primitives/strings/fileName/fileName.C b/src/OpenFOAM/primitives/strings/fileName/fileName.C
index 6a7e2056c6..dbf0a664e8 100644
--- a/src/OpenFOAM/primitives/strings/fileName/fileName.C
+++ b/src/OpenFOAM/primitives/strings/fileName/fileName.C
@@ -33,6 +33,7 @@ License
 
 const char* const Foam::fileName::typeName = "fileName";
 int Foam::fileName::debug(debug::debugSwitch(fileName::typeName, 0));
+const Foam::fileName Foam::fileName::null;
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/primitives/strings/string/string.C b/src/OpenFOAM/primitives/strings/string/string.C
index 84450a2242..e80b1e8221 100644
--- a/src/OpenFOAM/primitives/strings/string/string.C
+++ b/src/OpenFOAM/primitives/strings/string/string.C
@@ -29,6 +29,9 @@ License
 
 /* * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * */
 
+const char* const Foam::string::typeName = "string";
+int Foam::string::debug(debug::debugSwitch(string::typeName, 0));
+const Foam::string Foam::string::null;
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/primitives/strings/string/string.H b/src/OpenFOAM/primitives/strings/string/string.H
index cbf2d3bac8..714c72a335 100644
--- a/src/OpenFOAM/primitives/strings/string/string.H
+++ b/src/OpenFOAM/primitives/strings/string/string.H
@@ -80,6 +80,8 @@ public:
 
     // Static data members
 
+        static const char* const typeName;
+        static int debug;
         static const string null;
 
 
@@ -116,8 +118,6 @@ public:
 
     // Member Functions
 
-    // Access
-
         //- Count and return the number of a given character in the string
         size_type count(const char) const;
 
@@ -130,9 +130,6 @@ public:
         template<class String>
         static inline bool meta(const string&, const char quote='\\');
 
-
-    // Edit
-
         //- Strip invalid characters from the given string
         template<class String>
         static inline bool stripInvalid(string&);
@@ -145,7 +142,6 @@ public:
         template<class String>
         static inline string quotemeta(const string&, const char quote='\\');
 
-
         //- Replace first occurence of sub-string oldStr with newStr
         //  starting at start
         string& replace
diff --git a/src/OpenFOAM/primitives/strings/stringsGlobals.C b/src/OpenFOAM/primitives/strings/stringsGlobals.C
index 9d0d201a5e..7f4431da7c 100644
--- a/src/OpenFOAM/primitives/strings/stringsGlobals.C
+++ b/src/OpenFOAM/primitives/strings/stringsGlobals.C
@@ -38,8 +38,4 @@ Description
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const Foam::string Foam::string::null;
-const Foam::word Foam::word::null;
-const Foam::fileName Foam::fileName::null;
-
 // ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/strings/word/word.C b/src/OpenFOAM/primitives/strings/word/word.C
index c890046657..d80966ebf0 100644
--- a/src/OpenFOAM/primitives/strings/word/word.C
+++ b/src/OpenFOAM/primitives/strings/word/word.C
@@ -31,5 +31,6 @@ License
 
 const char* const Foam::word::typeName = "word";
 int Foam::word::debug(Foam::debug::debugSwitch(word::typeName, 0));
+const Foam::word Foam::word::null;
 
 // ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index 6277740b9a..2566a66e21 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -211,7 +211,6 @@ Foam::refinementSurfaces::refinementSurfaces
 
             minLevel_[globalRegionI] = iter();
             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
-            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 
             // Check validity
             if
@@ -231,6 +230,13 @@ Foam::refinementSurfaces::refinementSurfaces
                     << exit(FatalError);
             }
         }
+        forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
+        {
+            label globalRegionI = regionOffset_[surfI] + iter.key();
+
+            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
+        }
+
 
         //// Optional patch names and patch types
         //forAllConstIter(HashTable<word>, regionPatchName[surfI], iter)
@@ -387,7 +393,6 @@ Foam::refinementSurfaces::refinementSurfaces
 
             minLevel_[globalRegionI] = iter();
             maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
-            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
 
             // Check validity
             if
@@ -407,6 +412,12 @@ Foam::refinementSurfaces::refinementSurfaces
                     << exit(FatalError);
             }
         }
+        forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
+        {
+            label globalRegionI = regionOffset_[surfI] + iter.key();
+
+            perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
+        }
     }
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
index 11843958d5..a6d7b44bb5 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
@@ -26,7 +26,6 @@ License
 
 #include "SpalartAllmaras.H"
 #include "addToRunTimeSelectionTable.H"
-#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,7 +53,7 @@ tmp<volScalarField> SpalartAllmaras::fv1() const
 tmp<volScalarField> SpalartAllmaras::fv2() const
 {
     volScalarField chi = nuTilda_/nu();
-    return 1.0/pow3(scalar(1) + chi/Cv2_);
+    return 1/pow3(scalar(1) + chi/Cv2_);
 }
 
 
@@ -71,70 +70,68 @@ tmp<volScalarField> SpalartAllmaras::fv3() const
 }
 
 
-tmp<volScalarField> SpalartAllmaras::calcS(const volTensorField& gradU)
+tmp<volScalarField> SpalartAllmaras::S(const volTensorField& gradU) const
 {
-    return ::sqrt(2.0)*mag(skew(gradU));
+    return sqrt(2.0)*mag(skew(gradU));
 }
 
 
-tmp<volScalarField> SpalartAllmaras::calcSTilda(const volTensorField& gradU)
+tmp<volScalarField> SpalartAllmaras::STilda
+(
+    const volScalarField& S,
+    const volScalarField& dTilda
+) const
 {
-    return fv3()*calcS(gradU) + fv2()*nuTilda_/sqr(kappa_*dTilda_);
+    return fv3()*S + fv2()*nuTilda_/sqr(kappa_*dTilda);
 }
 
 
 tmp<volScalarField> SpalartAllmaras::r
 (
     const volScalarField& visc,
-    const volScalarField& S
+    const volScalarField& S,
+    const volScalarField& dTilda
 ) const
 {
-    tmp<volScalarField> tr
+    return min
     (
-        new volScalarField
-        (
-            min
-            (
-                visc
-               /(
-                    max
-                    (
-                        S,
-                        dimensionedScalar("SMALL", S.dimensions(), SMALL)
-                    )
-                   *sqr(kappa_*dTilda_)
-                  + dimensionedScalar
-                    (
-                        "ROOTVSMALL",
-                        dimensionSet(0, 2 , -1, 0, 0),
-                        ROOTVSMALL
-                    )
-                ),
-                scalar(10.0)
-            )
-        )
+        visc
+       /(
+           max
+           (
+               S,
+               dimensionedScalar("SMALL", S.dimensions(), SMALL)
+           )
+          *sqr(kappa_*dTilda)
+         + dimensionedScalar
+           (
+               "ROOTVSMALL",
+               dimensionSet(0, 2 , -1, 0, 0),
+               ROOTVSMALL
+           )
+        ),
+        scalar(10)
     );
-
-    return tr;
 }
 
 
-tmp<volScalarField> SpalartAllmaras::fw(const volScalarField& S) const
+tmp<volScalarField> SpalartAllmaras::fw
+(
+    const volScalarField& S,
+    const volScalarField& dTilda
+) const
 {
-    volScalarField r = this->r(nuTilda_, S);
+    volScalarField r = this->r(nuTilda_, S, dTilda);
 
     volScalarField g = r + Cw2_*(pow6(r) - r);
 
-    return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
+    return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
 }
 
 
-void SpalartAllmaras::dTildaUpdate(const volScalarField&)
+tmp<volScalarField> SpalartAllmaras::dTilda(const volScalarField&) const
 {
-    if (mesh_.changing())
-    {
-        dTilda_ = min(CDES_*delta(), wallDist(mesh_).y());
-    }
+    return min(CDES_*delta(), y_);
 }
 
 
@@ -242,6 +239,8 @@ SpalartAllmaras::SpalartAllmaras
         )
     ),
 
+    y_(mesh_),
+
     nuTilda_
     (
         IOobject
@@ -255,8 +254,6 @@ SpalartAllmaras::SpalartAllmaras
         mesh_
     ),
 
-    dTilda_(min(CDES_*delta(), wallDist(mesh_).y())),
-
     nuSgs_
     (
         IOobject
@@ -277,11 +274,15 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
 {
     LESModel::correct(gradU);
 
-    const volScalarField STilda = calcSTilda(gradU);
+    if (mesh_.changing())
+    {
+        y_.correct();
+        y_.boundaryField() = max(y_.boundaryField(), VSMALL);
+    }
 
-    const volScalarField S = calcS(gradU);
-
-    dTildaUpdate(S);
+    const volScalarField S = this->S(gradU);
+    const volScalarField dTilda = this->dTilda(S);
+    const volScalarField STilda = this->STilda(S, dTilda);
 
     fvScalarMatrix nuTildaEqn
     (
@@ -296,7 +297,7 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
       - alphaNut_*Cb2_*magSqr(fvc::grad(nuTilda_))
      ==
         Cb1_*STilda*nuTilda_
-      - fvm::Sp(Cw1_*fw(STilda)*nuTilda_/sqr(dTilda_), nuTilda_)
+      - fvm::Sp(Cw1_*fw(STilda, dTilda)*nuTilda_/sqr(dTilda), nuTilda_)
     );
 
     nuTildaEqn.relax();
@@ -310,9 +311,15 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
 }
 
 
+tmp<volScalarField> SpalartAllmaras::k() const
+{
+    return sqr(nuSgs()/ck_/dTilda(S(fvc::grad(U()))));
+}
+
+
 tmp<volScalarField> SpalartAllmaras::epsilon() const
 {
-    return 2.0*nuEff()*magSqr(symm(fvc::grad(U())));
+    return 2*nuEff()*magSqr(symm(fvc::grad(U())));
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
index aa721b1b5c..b813075b1d 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
@@ -38,6 +38,7 @@ SourceFiles
 
 #include "LESModel.H"
 #include "volFields.H"
+#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,30 +87,39 @@ protected:
 
         // Fields
 
+            wallDist y_;
             volScalarField nuTilda_;
-            volScalarField dTilda_;
             volScalarField nuSgs_;
 
 
     // Protected member functions
 
-        // Helper functions
+        virtual tmp<volScalarField> fv1() const;
+        virtual tmp<volScalarField> fv2() const;
+        virtual tmp<volScalarField> fv3() const;
+        virtual tmp<volScalarField> S(const volTensorField& gradU) const;
 
-            virtual tmp<volScalarField> fv1() const;
-            virtual tmp<volScalarField> fv2() const;
-            virtual tmp<volScalarField> fv3() const;
-            //-
-            virtual tmp<volScalarField> calcS(const volTensorField& gradU);
-            virtual tmp<volScalarField> calcSTilda(const volTensorField& gradU);
-            virtual tmp<volScalarField> r
-            (
-                const volScalarField& visc,
-                const volScalarField& S
-            ) const;
-            virtual tmp<volScalarField> fw(const volScalarField& S) const;
+        virtual tmp<volScalarField> STilda
+        (
+            const volScalarField& S,
+            const volScalarField& dTilda
+        ) const;
 
-        //- Length scale calculation
-        virtual void dTildaUpdate(const volScalarField& S);
+        virtual tmp<volScalarField> r
+        (
+            const volScalarField& visc,
+            const volScalarField& S,
+            const volScalarField& dTilda
+        ) const;
+
+        virtual tmp<volScalarField> fw
+        (
+            const volScalarField& S,
+            const volScalarField& dTilda
+        ) const;
+
+        //- Length scale
+        virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
 
 
 public:
@@ -138,10 +148,7 @@ public:
     // Member Functions
 
         //- Return SGS kinetic energy
-        virtual tmp<volScalarField> k() const
-        {
-            return sqr(nuSgs()/ck_/dTilda_);
-        }
+        virtual tmp<volScalarField> k() const;
 
         //- Return sub-grid disipation rate
         virtual tmp<volScalarField> epsilon() const;
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.C b/src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.C
index 6f5435a268..edba5b1d3c 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.C
@@ -26,7 +26,6 @@ License
 
 #include "SpalartAllmarasDDES.H"
 #include "addToRunTimeSelectionTable.H"
-#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,51 +49,42 @@ tmp<volScalarField> SpalartAllmarasDDES::rd
     const volScalarField& S
 ) const
 {
-    volScalarField d = wallDist(mesh_).y();
-
-    tmp<volScalarField> trd
+    return min
     (
-        new volScalarField
-        (
-            min
+        visc
+        /(
+            max
             (
-                visc
-               /(
-                    max
-                    (
-                        S,
-                        dimensionedScalar("SMALL", S.dimensions(), SMALL)
-                    )*sqr(kappa_*d)
-                  + dimensionedScalar
-                    (
-                        "ROOTVSMALL",
-                        dimensionSet(0, 2 , -1, 0, 0),
-                        ROOTVSMALL
-                    )
-                ), scalar(10.0)
+                S,
+                dimensionedScalar("SMALL", S.dimensions(), SMALL)
+            )*sqr(kappa_*y_)
+          + dimensionedScalar
+            (
+                "ROOTVSMALL",
+                dimensionSet(0, 2 , -1, 0, 0),
+                ROOTVSMALL
             )
-        )
+        ),
+        scalar(10)
     );
-
-    return trd;
 }
 
 
-tmp<volScalarField> SpalartAllmarasDDES::fd(const volScalarField& S)
+tmp<volScalarField> SpalartAllmarasDDES::fd(const volScalarField& S) const
 {
-    return 1.0 - tanh(pow3(8.0*rd(nuSgs_ + transport().nu(), S)));
+    return 1 - tanh(pow3(8*rd(nuEff(), S)));
 }
 
 
-void SpalartAllmarasDDES::dTildaUpdate(const volScalarField& S)
+tmp<volScalarField> SpalartAllmarasDDES::dTilda(const volScalarField& S) const
 {
-    dTilda_ =
-        wallDist(mesh_).y()
-      - fd(S)*max
-        (
-            dimensionedScalar("zero", dimLength, 0.0),
-            wallDist(mesh_).y() - CDES_*delta()
-        );
+    return max
+    (
+        y_
+      - fd(S)
+       *max(y_ - CDES_*delta(), dimensionedScalar("zero", dimLength, 0)),
+        dimensionedScalar("small", dimLength, SMALL)
+    );
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.H b/src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.H
index ec7eeda68c..52d7c99c6a 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.H
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmarasDDES/SpalartAllmarasDDES.H
@@ -62,6 +62,14 @@ class SpalartAllmarasDDES
 {
     // Private member functions
 
+        tmp<volScalarField> fd(const volScalarField& S) const;
+
+        tmp<volScalarField> rd
+        (
+            const volScalarField& visc,
+            const volScalarField& S
+        ) const;
+
         // Disallow default bitwise copy construct and assignment
         SpalartAllmarasDDES(const SpalartAllmarasDDES&);
         SpalartAllmarasDDES& operator=(const SpalartAllmarasDDES&);
@@ -71,15 +79,8 @@ protected:
 
     // Protected member functions
 
-        tmp<volScalarField> fd(const volScalarField& S);
-        tmp<volScalarField> rd
-        (
-            const volScalarField& visc,
-            const volScalarField& S
-        ) const;
-
-        //- Length scale calculation
-        virtual void dTildaUpdate(const volScalarField& S);
+        //- Length scale
+        virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
 
 
 public:
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C
index 6912872c9b..4fc7c02057 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C
@@ -26,7 +26,6 @@ License
 
 #include "SpalartAllmarasIDDES.H"
 #include "addToRunTimeSelectionTable.H"
-#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,22 +45,29 @@ addToRunTimeSelectionTable(LESModel, SpalartAllmarasIDDES, dictionary);
 
 tmp<volScalarField> SpalartAllmarasIDDES::alpha() const
 {
-    return
-        0.25
-      - wallDist(mesh_).y()
-       /dimensionedScalar("hMax", dimLength, max(cmptMax(delta())));
+    return max
+    (
+        0.25 - y_/dimensionedScalar("hMax", dimLength, max(cmptMax(delta()))),
+        scalar(-5)
+    );
 }
 
 
-tmp<volScalarField> SpalartAllmarasIDDES::ft(const volScalarField& S) const
+tmp<volScalarField> SpalartAllmarasIDDES::ft
+(
+    const volScalarField& S
+) const
 {
-    return tanh(pow3(sqr(ct_)*r(nuSgs_, S)));
+    return tanh(pow3(sqr(ct_)*rd(nuSgs_, S)));
 }
 
 
-tmp<volScalarField> SpalartAllmarasIDDES::fl(const volScalarField& S) const
+tmp<volScalarField> SpalartAllmarasIDDES::fl
+(
+    const volScalarField& S
+) const
 {
-    return tanh(pow(sqr(cl_)*r(transport().nu(), S), 10));
+    return tanh(pow(sqr(cl_)*rd(nu(), S), 10));
 }
 
 
@@ -71,33 +77,24 @@ tmp<volScalarField> SpalartAllmarasIDDES::rd
     const volScalarField& S
 ) const
 {
-    volScalarField d = wallDist(mesh_).y();
-
-    tmp<volScalarField> trd
+    return min
     (
-        new volScalarField
-        (
-            min
-            (
-                visc
-               /(
-                   max
-                   (
-                       S,
-                       dimensionedScalar("SMALL", S.dimensions(), SMALL)
-                   )*sqr(kappa_*d)
-                 + dimensionedScalar
-                   (
-                       "ROOTVSMALL",
-                       dimensionSet(0, 2 , -1, 0, 0),
-                       ROOTVSMALL
-                   )
-               ), scalar(10.0)
-            )
-        )
+        visc
+       /(
+           max
+           (
+               S,
+               dimensionedScalar("SMALL", S.dimensions(), SMALL)
+           )*sqr(kappa_*y_)
+         + dimensionedScalar
+           (
+               "ROOTVSMALL",
+               dimensionSet(0, 2 , -1, 0, 0),
+               ROOTVSMALL
+           )
+       ),
+       scalar(10)
     );
-
-    return trd;
 }
 
 
@@ -105,43 +102,38 @@ tmp<volScalarField> SpalartAllmarasIDDES::rd
 
 tmp<volScalarField> SpalartAllmarasIDDES::fd(const volScalarField& S) const
 {
-    return 1.0 - tanh(pow3(8.0*rd(nuSgs_+transport().nu(), S)));
+    return 1 - tanh(pow3(8*rd(nuEff(), S)));
 }
 
 
-void SpalartAllmarasIDDES::dTildaUpdate(const volScalarField& S)
+tmp<volScalarField> SpalartAllmarasIDDES::dTilda(const volScalarField& S) const
 {
     volScalarField alpha = this->alpha();
-
     volScalarField expTerm = exp(sqr(alpha));
 
     volScalarField fHill =
-        2.0*(pos(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
+        2*(pos(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
 
+    volScalarField fStep = min(2*pow(expTerm, -9.0), scalar(1));
+    volScalarField fHyb = max(1 - fd(S), fStep);
+    volScalarField fAmp = 1 - max(ft(S), fl(S));
+    volScalarField fRestore = max(fHill - 1, scalar(0))*fAmp;
 
-    volScalarField fStep = min(2.0*pow(expTerm, -9.0), scalar(1));
-    volScalarField fHyb = max(1.0 - fd(S), fStep);
-
-    volScalarField fAmp = 1.0 - max(ft(S), fl(S));
-
-    volScalarField fRestore = max(fHill - 1.0, scalar(0))*fAmp;
-
-    // volScalarField ft2 = IGNORING ft2 terms
-
+    // IGNORING ft2 terms
     volScalarField Psi = sqrt
     (
         min
         (
             scalar(100),
-            (1.0 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*fv2())/max(SMALL, fv1())
+            (1 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*fv2())/max(SMALL, fv1())
         )
     );
 
-    dTilda_ = max
+    return max
     (
-        dimensionedScalar("zero", dimLength, 0.0),
-        fHyb*(1.0 + fRestore*Psi)*wallDist(mesh_).y()
-      + (1.0 - fHyb)*CDES_*Psi*delta()
+        dimensionedScalar("SMALL", dimLength, SMALL),
+        fHyb*(1 + fRestore*Psi)*y_
+      + (1 - fHyb)*CDES_*Psi*delta()
     );
 }
 
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.H b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.H
index 815d2c28a6..8c7ca63ceb 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.H
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.H
@@ -69,12 +69,16 @@ class SpalartAllmarasIDDES
         tmp<volScalarField> alpha() const;
         tmp<volScalarField> ft(const volScalarField& S) const;
         tmp<volScalarField> fl(const volScalarField& S) const;
+
         tmp<volScalarField> rd
         (
             const volScalarField& visc,
             const volScalarField& S
         ) const;
 
+        //- Delay function
+        tmp<volScalarField> fd(const volScalarField& S) const;
+
         // Disallow default bitwise copy construct and assignment
         SpalartAllmarasIDDES(const SpalartAllmarasIDDES&);
         SpalartAllmarasIDDES& operator=(const SpalartAllmarasIDDES&);
@@ -84,11 +88,9 @@ protected:
 
     // Protected member functions
 
-        //- Delay function
-        tmp<volScalarField> fd(const volScalarField& S) const;
+        //- Length scale
+        virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
 
-        //- Length scale calculation
-        virtual void dTildaUpdate(const volScalarField& S);
 
 public: