diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index 6746b0b4a1..86eb165e04 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -96,4 +96,6 @@ flux/flux.C
 ddt2/ddt2.C
 zeroGradient/zeroGradient.C
 
+stabilityBlendingFactor/stabilityBlendingFactor.C
+
 LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects
diff --git a/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.C b/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.C
new file mode 100644
index 0000000000..a67339b00f
--- /dev/null
+++ b/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.C
@@ -0,0 +1,647 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "stabilityBlendingFactor.H"
+#include "addToRunTimeSelectionTable.H"
+#include "zeroGradientFvPatchFields.H"
+#include "fvcGrad.H"
+#include "surfaceInterpolate.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(stabilityBlendingFactor, 0);
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        stabilityBlendingFactor,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::stabilityBlendingFactor::writeFileHeader
+(
+    Ostream& os
+) const
+{
+    writeHeader(os, "Stabilization blending factor");
+    writeCommented(os, "Time");
+    writeTabbed(os, "Scheme1");
+    writeTabbed(os, "Scheme2");
+    writeTabbed(os, "Blended");
+    os  << endl;
+}
+
+bool Foam::functionObjects::stabilityBlendingFactor::calc()
+{
+    init(true);
+    return true;
+}
+
+
+bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
+{
+
+    const IOField<scalar>* residualPtr =
+        mesh_.lookupObjectPtr<IOField<scalar>>(residualName_);
+
+    if (residuals_)
+    {
+        if (!residualPtr)
+        {
+             WarningInFunction
+                << residualName_ << " is not available in the database."
+                << "It won't be considered for the blended field. " << nl
+                << "Add the corresponding FO (residuals). If the FO is already "
+                << "set, you need to wait for the first iteration."
+                << endl;
+
+        }
+        else
+        {
+            scalar meanRes = gAverage(mag(*residualPtr)) + VSMALL;
+
+            if (log)
+            {
+                Log << "    Average(mag(residuals)) :  " << meanRes << endl;
+            }
+
+            oldError_ = error_;
+            oldErrorIntegral_ = errorIntegral_;
+            error_ = mag(meanRes - mag(*residualPtr));
+            errorIntegral_ = oldErrorIntegral_ + 0.5*(error_ + oldError_);
+            const scalarField errorDifferential = error_ - oldError_;
+
+            const scalarField factorList
+            (
+                + P_*error_
+                + I_*errorIntegral_
+                + D_*errorDifferential
+            );
+
+            const scalarField indicatorResidual
+            (
+                max
+                (
+                    min
+                    (
+                        mag(factorList - meanRes)/(maxResidual_*meanRes),
+                        1.0
+                    ),
+                    0.0
+                )
+            );
+
+            forAll (indicator_, i)
+            {
+                indicator_[i] = indicatorResidual[i];
+            }
+        }
+    }
+
+    const volScalarField* nonOrthPtr =
+        mesh_.lookupObjectPtr<volScalarField>(nonOrthogonalityName_);
+
+    if (nonOrthogonality_)
+    {
+        indicator_ =
+        max
+        (
+            indicator_,
+            min
+            (
+                max
+                (
+                    0.0,
+                    (*nonOrthPtr - maxNonOrthogonality_)
+                   /(minNonOrthogonality_ - maxNonOrthogonality_)
+                ),
+                1.0
+            )
+        );
+
+        if (log)
+        {
+            Log << "    Max non-orthogonality :  " << max(*nonOrthPtr).value()
+                << endl;
+        }
+    }
+
+    const volScalarField* skewnessPtr =
+        mesh_.lookupObjectPtr<volScalarField>(skewnessName_);
+
+    if (skewness_)
+    {
+        indicator_ =
+        max
+        (
+            indicator_,
+            min
+            (
+                max
+                (
+                    0.0,
+                    (*skewnessPtr - maxSkewness_)
+                  / (minSkewness_ - maxSkewness_)
+                ),
+                1.0
+            )
+        );
+
+        if (log)
+        {
+            Log << "    Max skewness :  " << max(*skewnessPtr).value()
+                << endl;
+        }
+    }
+
+    const volScalarField* faceWeightsPtr =
+        mesh_.lookupObjectPtr<volScalarField>(faceWeightName_);
+
+    if (faceWeight_)
+    {
+        indicator_ =
+            max
+            (
+                indicator_,
+                min
+                (
+                    max
+                    (
+                        0.0,
+                        (minFaceWeight_ - *faceWeightsPtr)
+                      / (minFaceWeight_ - maxFaceWeight_)
+                    ),
+                    1.0
+                )
+            );
+
+        if (log)
+        {
+            Log << "    Min face weight:  " << min(*faceWeightsPtr).value()
+                << endl;
+        }
+    }
+
+
+    if (gradCc_)
+    {
+        tmp<volScalarField> magGradCCPtr
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "magGradCC",
+                    time_.timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh_,
+                dimensionedScalar("0", dimless, Zero),
+                zeroGradientFvPatchScalarField::typeName
+            )
+        );
+
+        for (direction i=0; i<vector::nComponents; i++)
+        {
+            // Create field with zero grad
+            volScalarField cci
+            (
+                IOobject
+                (
+                    "cc" + word(vector::componentNames[i]),
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh_,
+                dimensionedScalar("zero", dimLength, 0),
+                zeroGradientFvPatchScalarField::typeName
+            );
+            cci = mesh_.C().component(i);
+            cci.correctBoundaryConditions();
+            magGradCCPtr.ref() +=  mag(fvc::grad(cci)).ref();
+        }
+
+        if (log)
+        {
+            Log << "    Max magGradCc :  " << max(magGradCCPtr.ref()).value()
+                << endl;
+        }
+
+        indicator_ =
+            max
+            (
+                indicator_,
+                min
+                (
+                    max
+                    (
+                        0.0,
+                        (magGradCCPtr.ref() - maxGradCc_)
+                      / (minGradCc_ - maxGradCc_)
+                    ),
+                    1.0
+                )
+            );
+    }
+
+
+    const volVectorField* UNamePtr =
+        mesh_.lookupObjectPtr<volVectorField>(UName_);
+
+    if (Co_)
+    {
+        tmp<volScalarField> CoPtr
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "Co",
+                    time_.timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh_,
+                dimensionedScalar("0", dimless, Zero),
+                zeroGradientFvPatchScalarField::typeName
+            )
+        );
+
+        volScalarField& Co = CoPtr.ref();
+
+        Co.primitiveFieldRef() =
+            mesh_.time().deltaT()*mag(*UNamePtr)/pow(mesh_.V(), 1.0/3.0);
+
+        indicator_ =
+            max
+            (
+                indicator_,
+                min(max(0.0, (Co - Co1_)/(Co2_ - Co1_)), 1.0)
+            );
+
+        if (log)
+        {
+            Log << "    Max Co :  " << max(Co).value()
+                << endl;
+        }
+    }
+
+    indicator_.correctBoundaryConditions();
+    indicator_.min(1.0);
+    indicator_.max(0.0);
+
+    // Update the blended surface field
+    surfaceScalarField* surBlendedPtr =
+    (
+        mesh_.lookupObjectRefPtr<surfaceScalarField>(resultName_)
+    );
+
+    *surBlendedPtr = fvc::interpolate(indicator_);
+
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::stabilityBlendingFactor::stabilityBlendingFactor
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fieldExpression(name, runTime, dict),
+    writeFile(obr_, name, typeName, dict),
+    indicator_
+    (
+        IOobject
+        (
+            "blendedIndicator",
+            time_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("0", dimless, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    nonOrthogonality_(dict.lookupOrDefault<Switch>("switchNonOrtho", false)),
+    gradCc_(dict.lookupOrDefault<Switch>("switchGradCc", false)),
+    residuals_(dict.lookupOrDefault<Switch>("switchResiduals", false)),
+    faceWeight_(dict.lookupOrDefault<Switch>("switchFaceWeight", false)),
+    skewness_(dict.lookupOrDefault<Switch>("switchSkewness", false)),
+    Co_(dict.lookupOrDefault<Switch>("switchCo", false)),
+
+    maxNonOrthogonality_
+    (
+        dict.lookupOrDefault<scalar>("maxNonOrthogonality", 20.0)
+    ),
+    minNonOrthogonality_
+    (
+        dict.lookupOrDefault<scalar>("minNonOrthogonality", 60.0)
+    ),
+    maxGradCc_(dict.lookupOrDefault<scalar>("maxGradCc", 3.0)),
+    minGradCc_(dict.lookupOrDefault<scalar>("minGradCc", 4.0)),
+    maxResidual_(dict.lookupOrDefault<scalar>("maxResidual", 10.0)),
+    minFaceWeight_(dict.lookupOrDefault<scalar>("minFaceWeight", 0.3)),
+    maxFaceWeight_(dict.lookupOrDefault<scalar>("maxFaceWeight", 0.2)),
+    maxSkewness_(dict.lookupOrDefault<scalar>("maxSkewness", 2.0)),
+    minSkewness_(dict.lookupOrDefault<scalar>("minSkewness", 3.0)),
+    Co1_(dict.lookupOrDefault<scalar>("Co1", 1.0)),
+    Co2_(dict.lookupOrDefault<scalar>("Co2", 10.0)),
+
+    nonOrthogonalityName_
+    (
+        dict.lookupOrDefault<word>("nonOrthogonality", "nonOrthoAngle")
+    ),
+    faceWeightName_
+    (
+        dict.lookupOrDefault<word>("faceWeight", "faceWeight")
+    ),
+    skewnessName_
+    (
+        dict.lookupOrDefault<word>("skewness", "skewness")
+    ),
+    residualName_
+    (
+        dict.lookupOrDefault<word>("residual", "initialResidual:p")
+    ),
+    UName_
+    (
+         dict.lookupOrDefault<word>("U", "U")
+    ),
+
+    tolerance_(0.001),
+    error_(mesh_.nCells(), 0.0),
+    errorIntegral_(mesh_.nCells(), 0.0),
+    oldError_(mesh_.nCells(), 0.0),
+    oldErrorIntegral_(mesh_.nCells(), 0.0),
+    P_(dict.lookupOrDefault<scalar>("P", 3)),
+    I_(dict.lookupOrDefault<scalar>("I", 0.0)),
+    D_(dict.lookupOrDefault<scalar>("D", 0.25))
+{
+    read(dict);
+    setResultName(typeName, "");
+
+    tmp<surfaceScalarField> faceBlendedPtr
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                resultName_,
+                time_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("0", dimless, 0.0)
+        )
+    );
+    store(resultName_, faceBlendedPtr);
+
+    const volScalarField* nonOrthPtr =
+        mesh_.lookupObjectPtr<volScalarField>(nonOrthogonalityName_);
+
+    if (nonOrthogonality_)
+    {
+        if (!nonOrthPtr)
+        {
+            IOobject fieldHeader
+            (
+                nonOrthogonalityName_,
+                mesh_.time().constant(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            );
+
+            if (fieldHeader.typeHeaderOk<volScalarField>(true, true, false))
+            {
+                volScalarField* vfPtr(new volScalarField(fieldHeader, mesh_));
+                mesh_.objectRegistry::store(vfPtr);
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "Field : " << nonOrthogonalityName_ << " not found."
+                    << exit(FatalError);
+            }
+        }
+    }
+
+
+    const volScalarField* faceWeightsPtr =
+        mesh_.lookupObjectPtr<volScalarField>(faceWeightName_);
+
+    if (faceWeight_)
+    {
+        if (!faceWeightsPtr)
+        {
+            IOobject fieldHeader
+            (
+                faceWeightName_,
+                mesh_.time().constant(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            );
+
+            if (fieldHeader.typeHeaderOk<volScalarField>(true, true, false))
+            {
+                volScalarField* vfPtr(new volScalarField(fieldHeader, mesh_));
+                mesh_.objectRegistry::store(vfPtr);
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "Field : " << faceWeightName_ << " not found."
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    const volScalarField* skewnessPtr =
+        mesh_.lookupObjectPtr<volScalarField>(skewnessName_);
+
+    if (skewness_)
+    {
+        if (!skewnessPtr)
+        {
+            IOobject fieldHeader
+            (
+                skewnessName_,
+                mesh_.time().constant(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            );
+
+            if (fieldHeader.typeHeaderOk<volScalarField>(true, true, false))
+            {
+                volScalarField* vfPtr(new volScalarField(fieldHeader, mesh_));
+                mesh_.objectRegistry::store(vfPtr);
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "Field : " << skewnessName_ << " not found."
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    if (log)
+    {
+        indicator_.writeOpt() = IOobject::AUTO_WRITE;
+    }
+
+    init(true);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::stabilityBlendingFactor::~stabilityBlendingFactor()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::stabilityBlendingFactor::read
+(
+    const dictionary& dict
+)
+{
+    if (fieldExpression::read(dict) && writeFile::read(dict))
+    {
+        dict.lookup("switchNonOrtho") >> nonOrthogonality_;
+        dict.lookup("switchGradCc") >> gradCc_;
+        dict.lookup("switchResiduals") >> residuals_;
+        dict.lookup("switchFaceWeight") >> faceWeight_;
+        dict.lookup("switchSkewness") >> skewness_;
+        dict.lookup("switchCo") >> Co_;
+
+        dict.readIfPresent("maxNonOrthogonality", maxNonOrthogonality_);
+        dict.readIfPresent("maxGradCc", maxGradCc_);
+        dict.readIfPresent("maxResidual", maxResidual_);
+        dict.readIfPresent("maxSkewness", maxSkewness_);
+        dict.readIfPresent("maxFaceWeight", maxFaceWeight_);
+        dict.readIfPresent("Co2", Co2_);
+
+        dict.readIfPresent("minFaceWeight", minFaceWeight_);
+        dict.readIfPresent("minNonOrthogonality", minNonOrthogonality_);
+        dict.readIfPresent("minGradCc", minGradCc_);
+        dict.readIfPresent("minSkewness", minSkewness_);
+        dict.readIfPresent("Co1", Co1_);
+
+
+        dict.readIfPresent("P", P_);
+        dict.readIfPresent("I", I_);
+        dict.readIfPresent("D", D_);
+
+        tolerance_ = 0.001;
+        if
+        (
+            dict.readIfPresent("tolerance", tolerance_)
+         && (tolerance_ < 0 || tolerance_ > 1)
+        )
+        {
+            FatalErrorInFunction
+                << "tolerance must be in the range 0 to 1.  Supplied value: "
+                << tolerance_ << exit(FatalError);
+        }
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+bool Foam::functionObjects::stabilityBlendingFactor::write()
+{
+    // Generate scheme statistics
+    label nCellsScheme1 = 0;
+    label nCellsScheme2 = 0;
+    label nCellsBlended = 0;
+    forAll(indicator_, celli)
+    {
+        scalar i = indicator_[celli];
+
+        if (i < tolerance_)
+        {
+            nCellsScheme2++;
+        }
+        else if (i > (1 - tolerance_))
+        {
+            nCellsScheme1++;
+        }
+        else
+        {
+            nCellsBlended++;
+        }
+    }
+
+    reduce(nCellsScheme1, sumOp<label>());
+    reduce(nCellsScheme2, sumOp<label>());
+    reduce(nCellsBlended, sumOp<label>());
+
+    Log << "    scheme 1 cells :  " << nCellsScheme1 << nl
+        << "    scheme 2 cells :  " << nCellsScheme2 << nl
+        << "    blended cells  :  " << nCellsBlended << nl
+        << endl;
+
+    writeTime(file());
+
+    file()
+        << token::TAB << nCellsScheme1
+        << token::TAB << nCellsScheme2
+        << token::TAB << nCellsBlended
+        << endl;
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.H b/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.H
new file mode 100644
index 0000000000..49bd1dac0e
--- /dev/null
+++ b/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.H
@@ -0,0 +1,467 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::stabilityBlendingFactor
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    Calculates and outputs the stabilityBlendingFactor to be used by the
+    local blended convection scheme. The output is a surface field weight
+    between 0-1
+
+    The weight of a blended scheme is given by a function of the blending
+    factor, f:
+
+    \f[
+        weight = f scheme1 + (1 - f) scheme2
+    \f]
+
+    The factor is calculated based on five criteria:
+
+    1) mesh non-orthogonality field
+    2) magnitude of cell centres gradient
+    3) convergence rate of residuals
+    4) faceWeight
+    5) skewness
+    6) Co number
+
+    The user can enable them individually.
+
+    For option 1, the following relation is used:
+        \f[
+            fNon =
+                min
+                (
+                    max
+                    (
+                        0.0,
+                        (non-orthogonality - maxNonOrthogonality)
+                       /(minNonOrthogonality - maxNonOrthogonality)
+                    ),
+                    1.0
+                )
+        \f[
+
+    For option 2, the following relation is used:
+        \f[
+            fMagGradCc =
+                min
+                (
+                    max
+                    (
+                        0.0,
+                        (magGradCC - maxGradCc)
+                      / (minGradCc - maxGradCc)
+                    ),
+                    1.0
+                )
+        \f[
+    Note that  magGradCC is equal to 3 for ortoghonal meshes
+
+    For option 3 a PID control is used in order to control residual
+    unbouded fluctuations for individual cells.
+
+        \f[
+            factor =
+                P*residual
+              + I*residualIntegral
+              + D*residualDifferential
+        \f[
+
+        where P, I and D are user inputs
+
+    The following relation is used:
+    \f[
+        fRes = (factor - meanRes)/(maxRes*meanRes);
+    \f[
+
+        where:
+            meanRes = average(residual)
+            maxRes is an user input
+
+
+    fRes will blend more towards one as the cell residual is larger then
+    the domain mean residuals.
+
+
+    For option 4 , the following relation is used:
+
+        \f[
+            ffaceWeight = min
+            (
+                max
+                (
+                    0.0,
+                    (minFaceWeight - faceWeights)
+                  / (minFaceWeight - maxFaceWeight)
+                ),
+                1.0
+            )
+        \f[
+    Note that faceWeights for a orthogonal mesh is 0.5.
+
+
+    For option 5 , the following relation is used:
+
+        \f[
+            fskewness =
+            min
+            (
+                max
+                (
+                    0.0,
+                    (skewness    - maxSkewness)
+                  / (minSkewness - maxSkewness)
+                ),
+                1.0
+            )
+        \f[
+
+
+    For option 6 , the following relation is used:
+
+    \f[
+        fCoWeight = min(max((Co - Co1)/(Co2 - Co1), 0), 1)
+    \f]
+    where
+
+    \vartable
+        Co1 | Courant number below which scheme2 is used
+        Co2 | Courant number above which scheme1 is used
+    \endvartable
+
+    The final factor is determined by:
+
+    \f[
+        f = max(fNon, fMagGradCc, fRes, ffaceWeight, fskewness, fCoWeight)
+    \f]
+
+    An indicator (volume) field, named blendedIndicator is generated if the log
+    flag is on:
+    - 1 represent scheme1 as active,
+    - 0 represent scheme2 as active.
+
+    Additional reporting is written to the standard output, providing
+    statistics as to the number of cells used by each scheme.
+
+Usage
+    Example of function object specification to calculate the blending factor:
+    \verbatim
+    stabilityBlendingFactor1
+    {
+        type                stabilityBlendingFactor;
+        libs                ("libfieldFunctionObjects.so");
+        log                 true;
+        switchNonOrtho      yes;
+        switchGradCc        no;
+        switchResiduals     yes;
+        switchFaceWeight    no;
+        switchSkewness      no;
+        switchCo            no;
+
+        maxNonOrthogonality 20;
+        minNonOrthogonality 60;
+
+        maxGradCc           3;
+        minGradCc           4;
+
+        maxFaceWeight       0.3;
+        minFaceWeight       0.2;
+
+        maxSkewness         2;
+        minSkewness         3;
+
+        maxResidual         10;
+
+        result              UBlendingFactor;
+        residual            initialResidual:p;
+        P                   1.5;
+        I                   0;
+        D                   0.5;
+
+        Co1                 1;
+        Co2                 10;
+        ...
+        ...
+        field               U;
+    }
+
+    Example of function object specification to calculate the residuals used
+    by stabilityBlendingFactor. The following writes 'initialResidual:p'
+    field
+
+        residuals
+        {
+            type            residuals;
+            libs            ("libutilityFunctionObjects.so");
+            writeFields     true;
+            fields          (p);
+        }
+
+    \endverbatim
+
+    Where the entries comprise:
+    \table
+        Property     | Description             | Required    | Default value
+        type         | Type name: stabilityBlendingFactor | yes       |
+
+
+        log          | Log to standard output  | no          | yes
+        switchNonOrtho | non-orthogonal method | no | false
+        switchGradCc | cell centre gradient method   | no | false
+        switchResiduals | residual evolution method   | no | false
+        switchFaceWeight | face weight method | no    | false
+        switchSkewness | skewness method | no    | false
+        switchCo       | Co blended |  no        | false
+
+        maxNonOrthogonality| maximum non-orthogonal for scheme2 | no | 20
+        minNonOrthogonality| minimum non-orthogonal for scheme1 | no | 60
+
+        maxGradCc| maximum gradient for scheme2 | no | 2
+        minGradCc| minimum gradient for scheme1 | no | 4
+
+        maxResidual| maximum residual-mean ratio for scheme1 | no | 10
+        P       |   proportional factor for PID     | no | 3
+        I       |   integral factor for PID         | no | 0
+        D       |   differential factor for PID     | no | 0.25
+
+        maxFaceWeight | maximum face weight for scheme1 | no | 0.2
+        minFaceWeight | minimum face weight for scheme2 | no | 0.3
+
+        maxSkewness | maximum skewness for scheme2 | no | 2
+        minSkewness | minimum skewness for scheme1 | no | 3
+
+
+        faceWeight | Name of the faceWeight field | no | faceWeight
+        skewness | Name of the skewness field | no | skewness
+        nonOrthogonality | Name of the non-orthogonal field | no |
+            nonOrthoAngle
+        residual    | Name of the residual field | no | initialResidual:p
+        U           | Name of the flux field for Co blended | no | U
+
+
+        Co1 | Courant number below which scheme2 is used | no | 1
+        Co2 | Courant number above which scheme1 is used | no | 10
+
+
+        tolerance    | Tolerance for number of blended cells | no | 0.001
+        field        | Name of field to evaluate | yes       |
+        result    | Name of surface field to be used in the localBlended scheme
+        | yes
+
+    \endtable
+
+See also
+    Foam::functionObjects::fieldExpression
+    Foam::functionObjects::fvMeshFunctionObject
+    Foam::functionObjects::writeFile
+
+SourceFiles
+    stabilityBlendingFactor.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_stabilityBlendingFactor_H
+#define functionObjects_stabilityBlendingFactor_H
+
+#include "fieldExpression.H"
+#include "writeFile.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class stabilityBlendingFactor Declaration
+\*---------------------------------------------------------------------------*/
+
+class stabilityBlendingFactor
+:
+    public fieldExpression,
+    public writeFile
+{
+    // Private member data
+
+        //- Cell based blended indicator
+        volScalarField indicator_;
+
+        // Switches
+
+            //- Switch for non-orthogonality
+            Switch nonOrthogonality_;
+
+            //- Switch for grad of cell centres
+            Switch gradCc_;
+
+            //- Switch for residuals
+            Switch residuals_;
+
+            //- Switch for face weight
+            Switch faceWeight_;
+
+            //- Switch for skewness
+            Switch skewness_;
+
+            //- Switch for Co
+            Switch Co_;
+
+
+        // Lower and upper limits
+
+            //- Maximum non-orthogonality for fully scheme 2
+            scalar maxNonOrthogonality_;
+
+            //- Minimum non-orthogonality for fully scheme 1
+            scalar minNonOrthogonality_;
+
+            //- Maximum  gradcc for fully scheme 2
+            scalar maxGradCc_;
+
+            //- Minimum  gradcc for fully scheme 1
+            scalar minGradCc_;
+
+            //- Maximum ratio to average residual for scheme 2
+            scalar maxResidual_;
+
+            //- Minimum face weight for fully scheme 2
+            scalar minFaceWeight_;
+
+            //- Maximum face weight for fully scheme 1
+            scalar maxFaceWeight_;
+
+            //- Maximum skewness for fully scheme 2
+            scalar maxSkewness_;
+
+            //- Minimum skewness for fully scheme 1
+            scalar minSkewness_;
+
+            //- Maximum Co for fully scheme 2
+            scalar Co1_;
+
+            //- Minimum Co for fully scheme 1
+            scalar Co2_;
+
+
+        // File names
+
+            //- Name of the non-orthogonalit field
+            word nonOrthogonalityName_;
+
+            //- Name of the face weight field
+            word faceWeightName_;
+
+            //- Name of the skewnes field
+            word skewnessName_;
+
+            //- Name of the residual field
+            word residualName_;
+
+            //- Name of the U used for Co based blended
+            word UName_;
+
+
+        //- Tolerance used when calculating the number of blended cells
+        scalar tolerance_;
+
+
+        //- Error fields
+        scalarField error_;
+        scalarField errorIntegral_;
+        scalarField oldError_;
+        scalarField oldErrorIntegral_;
+
+        //- Proportional gain
+        scalar P_;
+
+        //- Integral gain
+        scalar I_;
+
+        //- Derivative gain
+        scalar D_;
+
+
+    // Private Member Functions
+
+        //- Init fields
+        bool init(bool first);
+
+        //- Calculate the blending factor field and return true if successful
+        virtual bool calc();
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Write the file header
+        virtual void writeFileHeader(Ostream& os) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("stabilityBlendingFactor");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        stabilityBlendingFactor
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~stabilityBlendingFactor();
+
+
+    // Member Functions
+
+        //- Read the stabilityBlendingFactor data
+        virtual bool read(const dictionary&);
+
+        //- Write the stabilityBlendingFactor
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/motorBike/Allrun b/tutorials/incompressible/simpleFoam/motorBike/Allrun
index 268780fb13..e0e6038798 100755
--- a/tutorials/incompressible/simpleFoam/motorBike/Allrun
+++ b/tutorials/incompressible/simpleFoam/motorBike/Allrun
@@ -32,6 +32,8 @@ restore0Dir -processor
 
 runParallel $decompDict patchSummary
 runParallel $decompDict potentialFoam
+runParallel $decompDict checkMesh -writeFields '(nonOrthoAngle)' -constant
+
 runParallel $decompDict $(getApplication)
 
 runApplication reconstructParMesh -constant
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
index a626f4d688..30048b9ae6 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/controlDict
@@ -16,7 +16,7 @@ FoamFile
 
 application     simpleFoam;
 
-startFrom       latestTime;
+startFrom       startTime;
 
 startTime       0;
 
@@ -51,6 +51,7 @@ functions
     #include "cuttingPlane"
     #include "forceCoeffs"
     #include "ensightWrite"
+    #include "stabilizationSchemes"
 }
 
 
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes b/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
index 12fdd994df..88d86cc630 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
@@ -28,7 +28,7 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      bounded Gauss linearUpwindV grad(U);
+    div(phi,U)      bounded Gauss localBlended linearUpwindV grad(U) linear;
     div(phi,k)      bounded Gauss upwind;
     div(phi,omega)  bounded Gauss upwind;
     div((nuEff*dev2(T(grad(U))))) Gauss linear;
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/stabilizationSchemes b/tutorials/incompressible/simpleFoam/motorBike/system/stabilizationSchemes
new file mode 100644
index 0000000000..2de422fd76
--- /dev/null
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/stabilizationSchemes
@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+residuals
+{
+    type            residuals;
+    libs            ("libutilityFunctionObjects.so");
+    writeFields     true;
+    writeControl    outputTime;
+    fields          (p);
+}
+
+blendingFactor
+{
+    type                stabilityBlendingFactor;
+    libs                ("libfieldFunctionObjects.so");
+    log                 true;
+    switchNonOrtho      yes;
+    switchGradCc        no;
+    switchResiduals     yes;
+    switchSkewness      no;
+    switchFaceWeight    no;
+    switchCo            no;
+
+    maxNonOrthogonality 20;
+    minNonOrthogonality 60;
+
+    maxGradCc           3;
+    minGradCc           4;
+
+    maxResidual         100;
+
+    P                   5;
+    I                   0.01;
+    D                   0.5;
+
+    minFaceWeight       0.3;
+    maxFaceWeight       0.2;
+
+    Co1                 1;
+    Co2                 2;
+
+    maxSkewness         2;
+    minSkewness         3;
+
+
+    field               U;
+    result              UBlendingFactor;
+}
+
+// ************************************************************************* //