diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.C b/applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.C
index 944ff9cb48726bf01fd57bc17cc9c1106989898d..ffea6709521b973551b9b52ddea76dbebb4d622e 100644
--- a/applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.C
+++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.C
@@ -158,7 +158,6 @@ void Foam::createShellMesh::calcPointRegions
 }
 
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::createShellMesh::createShellMesh
@@ -184,7 +183,6 @@ Foam::createShellMesh::createShellMesh
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
 void Foam::createShellMesh::setRefinement
 (
     const pointField& thickness,
diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C b/applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C
index 6a3b7dd7eadd0a926c8dabc1b9823be7d101caeb..4cf3c67ff998f2f1f6b9f8146ef03dc0c59f6242 100644
--- a/applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C
+++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C
@@ -130,6 +130,7 @@ Usage
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "cyclicPolyPatch.H"
+#include "syncTools.H"
 
 using namespace Foam;
 
@@ -595,6 +596,237 @@ void createDummyFvMeshFiles(const polyMesh& mesh, const word& regionName)
 }
 
 
+//XXXXXXXXX
+label findUncoveredPatchFace
+(
+    const fvMesh& mesh,
+    const UIndirectList<label>& extrudeMeshFaces,// mesh faces that are extruded
+    const label meshEdgeI                       // mesh edge
+)
+{
+    // Make set of extruded faces.
+    labelHashSet extrudeFaceSet(extrudeMeshFaces.size());
+    forAll(extrudeMeshFaces, i)
+    {
+        extrudeFaceSet.insert(extrudeMeshFaces[i]);
+    }
+
+    label patchI = -1;
+
+    const labelList& eFaces = mesh.edgeFaces()[meshEdgeI];
+    forAll(eFaces, i)
+    {
+        label faceI = eFaces[i];
+        if (!mesh.isInternalFace(faceI) && !extrudeFaceSet.found(faceI))
+        {
+            patchI = mesh.boundaryMesh().whichPatch(faceI);
+            break;
+        }
+    }
+    return patchI;
+}
+// Count the number of faces in patches that need to be created
+void countExtrudePatches
+(
+    const fvMesh& mesh,
+    const primitiveFacePatch& extrudePatch,
+    const label nZones,
+    const labelList& zoneID,
+    const labelList& extrudeMeshFaces,
+    const labelList& extrudeMeshEdges,
+
+    labelList& zoneSidePatch,
+    labelList& zoneZonePatch
+)
+{
+    const labelListList& edgeFaces = extrudePatch.edgeFaces();
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+        if (eFaces.size() == 2)
+        {
+            label zone0 = zoneID[eFaces[0]];
+            label zone1 = zoneID[eFaces[1]];
+
+            if (zone0 != zone1)
+            {
+                label minZone = min(zone0,zone1);
+                label maxZone = max(zone0,zone1);
+                zoneZonePatch[minZone*nZones+maxZone]++;
+            }
+        }
+        else
+        {
+            // Check whether we are on a mesh edge with external patches. If
+            // so choose any uncovered one. If none found put face in
+            // undetermined zone 'side' patch
+
+            label patchI = findUncoveredPatchFace
+            (
+                mesh,
+                UIndirectList<label>(extrudeMeshFaces, eFaces),
+                extrudeMeshEdges[edgeI]
+            );
+
+            if (patchI == -1)
+            {
+                // Determine the min zone of all connected zones.
+                label minZone = zoneID[eFaces[0]];
+                for (label i = 1; i < eFaces.size(); i++)
+                {
+                    minZone = min(minZone, zoneID[eFaces[i]]);
+                }
+                zoneSidePatch[minZone]++;
+            }
+        }
+    }
+    Pstream::listCombineGather(zoneSidePatch, plusEqOp<label>());
+    Pstream::listCombineScatter(zoneSidePatch);
+    Pstream::listCombineGather(zoneZonePatch, plusEqOp<label>());
+    Pstream::listCombineScatter(zoneZonePatch);
+}
+bool lessThan(const point& x, const point& y)
+{
+    for (direction dir = 0; dir < point::nComponents; dir++)
+    {
+        if (x[dir] < y[dir]) return true;
+        if (x[dir] > y[dir]) return false;
+    }
+    return false;
+}
+
+class minEqVectorListOp
+{
+    public:
+    void operator()(List<vector>& x, const List<vector>& y) const
+    {
+        if (y.size())
+        {
+            if (x.size())
+            {
+                forAll(x, i)
+                {
+                    if (lessThan(y[i], x[i]))
+                    {
+                        x[i] = y[i];
+                    }
+                }
+            }
+            else
+            {
+                x = y;
+            }
+        }
+    }
+};
+// Constrain&sync normals on points that are on coupled patches.
+void constrainCoupledNormals
+(
+    const fvMesh& mesh,
+    const primitiveFacePatch& extrudePatch,
+    const labelList& regionToPoint,
+
+    vectorField& regionNormals
+)
+{
+    // Invert regionToPoint to create pointToRegions.
+    labelListList pointToRegions
+    (
+        invertOneToMany
+        (
+            extrudePatch.nPoints(),
+            regionToPoint
+        )
+    );
+    // Sort acc. to region so (hopefully) coupled points will do the same.
+    forAll(pointToRegions, pointI)
+    {
+        sort(pointToRegions[pointI]);
+    }
+
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Constrain displacement on cyclic patches
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Note: bit contentious to always do this on cyclic - should user use
+    // different patch type, e.g. 'cyclicSlip' instead?
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (isA<cyclicPolyPatch>(pp))
+        {
+            forAll(pp.meshPoints(), pointI)
+            {
+                Map<label>::const_iterator fnd =
+                    extrudePatch.meshPointMap().find
+                    (
+                        pp.meshPoints()[pointI]
+                    );
+                if (fnd != extrudePatch.meshPointMap().end())
+                {
+                    // fnd() is a point on this cyclic.
+                    const vector& cycNormal = pp.pointNormals()[pointI];
+                    const labelList& pRegions = pointToRegions[fnd()];
+                    forAll(pRegions, i)
+                    {
+                        // Remove cyclic normal component.
+                        vector& regionNormal = regionNormals[pRegions[i]];
+                        regionNormal -= (regionNormal&cycNormal)*cycNormal;
+                    }
+                }
+            }
+        }
+    }
+
+
+    // Synchronise regionNormals
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Re-work regionNormals into multiple normals per point
+    List<List<vector> > pointNormals(mesh.nPoints());
+    forAll(pointToRegions, pointI)
+    {
+        const labelList& pRegions = pointToRegions[pointI];
+
+        label meshPointI = extrudePatch.meshPoints()[pointI];
+        List<vector>& pNormals = pointNormals[meshPointI];
+        pNormals.setSize(pRegions.size());
+        forAll(pRegions, i)
+        {
+            pNormals[i] = regionNormals[pRegions[i]];
+        }
+    }
+
+    // Synchronise
+    syncTools::syncPointList
+    (
+        mesh,
+        pointNormals,
+        minEqVectorListOp(),
+        List<vector>(),         // nullValue
+        false                   // applySeparation
+    );
+
+    // Re-work back into regionNormals
+    forAll(pointToRegions, pointI)
+    {
+        const labelList& pRegions = pointToRegions[pointI];
+
+        label meshPointI = extrudePatch.meshPoints()[pointI];
+        const List<vector>& pNormals = pointNormals[meshPointI];
+        forAll(pRegions, i)
+        {
+            regionNormals[pRegions[i]] = pNormals[i];
+        }
+    }
+}
+//XXXXXXXXX
+
+
 tmp<pointField> calcOffset
 (
     const primitiveFacePatch& extrudePatch,
@@ -739,6 +971,16 @@ int main(int argc, char *argv[])
         << endl;
 
 
+    // Determine corresponding mesh edges
+    const labelList extrudeMeshEdges
+    (
+        extrudePatch.meshEdges
+        (
+            mesh.edges(),
+            mesh.pointEdges()
+        )
+    );
+
 
 
     // Check whether the zone is internal or external faces to determine
@@ -772,7 +1014,7 @@ int main(int argc, char *argv[])
             Info<< "FaceZone " << fz.name() << " has boundary faces" << endl;
         }
     }
-
+    Info<< endl;
 
 
 
@@ -866,36 +1108,18 @@ int main(int argc, char *argv[])
     labelList zoneSidePatch(faceZones.size(), 0);
     labelList zoneZonePatch(faceZones.size()*faceZones.size(), 0);
 
-    forAll(edgeFaces, edgeI)
-    {
-        const labelList& eFaces = edgeFaces[edgeI];
-        if (eFaces.size() == 2)
-        {
-            label zone0 = zoneID[eFaces[0]];
-            label zone1 = zoneID[eFaces[1]];
+    countExtrudePatches
+    (
+        mesh,
+        extrudePatch,
+        faceZones.size(),
+        zoneID,
+        extrudeMeshFaces,
+        extrudeMeshEdges,
 
-            if (zone0 != zone1)
-            {
-                label minZone = min(zone0,zone1);
-                label maxZone = max(zone0,zone1);
-                zoneZonePatch[minZone*faceZones.size()+maxZone]++;
-            }
-        }
-        else
-        {
-            // Determine the min zone of all connected zones.
-            label minZone = zoneID[eFaces[0]];
-            for (label i = 1; i < eFaces.size(); i++)
-            {
-                minZone = min(minZone, zoneID[eFaces[i]]);
-            }
-            zoneSidePatch[minZone]++;
-        }
-    }
-    Pstream::listCombineGather(zoneSidePatch, plusEqOp<label>());
-    Pstream::listCombineScatter(zoneSidePatch);
-    Pstream::listCombineGather(zoneZonePatch, plusEqOp<label>());
-    Pstream::listCombineScatter(zoneZonePatch);
+        zoneSidePatch,
+        zoneZonePatch
+    );
 
     // Now check which patches to add.
     Info<< "Adding patches for edges on zones:" << nl << nl
@@ -980,6 +1204,7 @@ int main(int argc, char *argv[])
     // Is edge an non-manifold edge
     PackedBoolList nonManifoldEdge(extrudePatch.nEdges());
 
+    // Note: logic has to be same as in countExtrudePatches.
     forAll(edgeFaces, edgeI)
     {
         const labelList& eFaces = edgeFaces[edgeI];
@@ -1004,65 +1229,31 @@ int main(int argc, char *argv[])
         }
         else
         {
-            ePatches.setSize(eFaces.size());
-            forAll(eFaces, i)
-            {
-                ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]];
-            }
-            nonManifoldEdge[edgeI] = 1;
-        }
-    }
-
-    // Override constraint types
-    {
-        const edgeList& extrudeEdges = extrudePatch.edges();
-        const labelList& extrudeMeshPts = extrudePatch.meshPoints();
-
-        // Map from mesh edge to local patch edge index
-        EdgeMap<label> extrudeMeshEdges(extrudePatch.nEdges());
+            label patchI = findUncoveredPatchFace
+            (
+                mesh,
+                UIndirectList<label>(extrudeMeshFaces, eFaces),
+                extrudeMeshEdges[edgeI]
+            );
 
-        forAll(extrudeEdges, edgeI)
-        {
-            if (edgeFaces[edgeI].size() == 1)
+            if (patchI != -1)
             {
-                const edge& e = extrudeEdges[edgeI];
-                const edge meshE(extrudeMeshPts[e[0]], extrudeMeshPts[e[1]]);
-                extrudeMeshEdges.insert(meshE, edgeI);
+                ePatches.setSize(eFaces.size(), patchI);
             }
-        }
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            if (polyPatch::constraintType(pp.type()))
+            else
             {
-                const edgeList& edges = pp.edges();
-
-                forAll(edges, ppEdgeI)
+                ePatches.setSize(eFaces.size());
+                forAll(eFaces, i)
                 {
-                    const edge& e = edges[ppEdgeI];
-                    edge meshE(pp.meshPoints()[e[0]], pp.meshPoints()[e[1]]);
-                    EdgeMap<label>::const_iterator iter = extrudeMeshEdges.find
-                    (
-                        meshE
-                    );
-                    if (iter != extrudeMeshEdges.end())
-                    {
-                        label extrudeEdgeI = iter();
-                        extrudeEdgePatches[extrudeEdgeI] = labelList
-                        (
-                            edgeFaces[extrudeEdgeI].size(),
-                            patchI
-                        );
-                    }
+                    ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]];
                 }
             }
+            nonManifoldEdge[edgeI] = 1;
         }
     }
 
 
 
-
     // Assign point regions
     // ~~~~~~~~~~~~~~~~~~~~
 
@@ -1075,6 +1266,7 @@ int main(int argc, char *argv[])
     (
         extrudePatch,
         nonManifoldEdge,
+
         pointRegions,
         regionPoints
     );
@@ -1096,6 +1288,18 @@ int main(int argc, char *argv[])
     }
     regionNormals /= mag(regionNormals);
 
+
+    // Constrain&sync normals on points that are on coupled patches.
+    constrainCoupledNormals
+    (
+        mesh,
+        extrudePatch,
+        regionPoints,
+
+        regionNormals
+    );
+
+
     // For debugging: dump hedgehog plot of normals
     {
         OFstream str(runTime.path()/"regionNormals.obj");
@@ -1112,7 +1316,7 @@ int main(int argc, char *argv[])
 
                 meshTools::writeOBJ(str, pt);
                 vertI++;
-                meshTools::writeOBJ(str, pt+0.01*regionNormals[region]);
+                meshTools::writeOBJ(str, pt+thickness*regionNormals[region]);
                 vertI++;
                 str << "l " << vertI-1 << ' ' << vertI << nl;
             }
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 4d89dedde191f2b1bca08e07a3bf8b586910afb7..31c274ea38d9fd0372992665f3b7edd61d975890 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -346,7 +346,7 @@ bool Foam::KinematicParcel<ParcelType>::hitPatch
 {
     ParcelType& p = static_cast<ParcelType&>(*this);
 
-    // Invoke poost-processing mdoel
+    // Invoke post-processing model
     td.cloud().postProcessing().postPatch(p, patchI);
 
     // Invoke surface film model
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index 84d959976d769204edb8f466f3d14b5dda1382ab..a232a5f0609c754717d9d545817911d299b3054f 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -212,7 +212,12 @@ Foam::scalar Foam::InjectionModel<CloudType>::setNumberOfParticles
         }
         case pbNumber:
         {
-            nP = massTotal_/(rho*volumeTotal_*parcels);
+            nP = massTotal_/(rho*volumeTotal_);
+            break;
+        }
+        case pbFixed:
+        {
+            nP = nParticlesFixed_;
             break;
         }
         default:
@@ -285,6 +290,7 @@ Foam::InjectionModel<CloudType>::InjectionModel(CloudType& owner)
     nInjections_(0),
     parcelsAddedTotal_(0),
     parcelBasis_(pbNumber),
+    nParticlesFixed_(0.0),
     time0_(0.0),
     timeStep0_(0.0)
 {
@@ -310,6 +316,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
     nInjections_(0),
     parcelsAddedTotal_(0),
     parcelBasis_(pbNumber),
+    nParticlesFixed_(0.0),
     time0_(owner.db().time().value()),
     timeStep0_(0.0)
 {
@@ -320,6 +327,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
         << endl;
 
     const word parcelBasisType = coeffDict_.lookup("parcelBasisType");
+
     if (parcelBasisType == "mass")
     {
         parcelBasis_ = pbMass;
@@ -328,6 +336,16 @@ Foam::InjectionModel<CloudType>::InjectionModel
     {
         parcelBasis_ = pbNumber;
     }
+    else if (parcelBasisType == "fixed")
+    {
+        parcelBasis_ = pbFixed;
+
+        Info<< "    Choosing nParticles to be a fixed value, massTotal "
+            << "variable now does not determine anything."
+            << endl;
+
+        nParticlesFixed_ = readScalar(coeffDict_.lookup("nParticles"));
+    }
     else
     {
         FatalErrorIn
@@ -338,7 +356,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
                 "CloudType&, "
                 "const word&"
             ")"
-        )<< "parcelBasisType must be either 'number' or 'mass'" << nl
+        )<< "parcelBasisType must be either 'number', 'mass' or 'fixed'" << nl
          << exit(FatalError);
     }
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
index d7309f5a0a8081f67e85a20041e10688a6bed50c..c6fb6e424bc3270da4d41f4f546a493a6a79c1c3 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
@@ -73,7 +73,8 @@ public:
         enum parcelBasis
         {
             pbNumber,
-            pbMass
+            pbMass,
+            pbFixed
         };
 
 
@@ -137,6 +138,10 @@ protected:
             //- Parcel basis enumeration
             parcelBasis parcelBasis_;
 
+            //- nParticles to assign to parcels when the 'fixed' basis
+            //  is selected
+            scalar nParticlesFixed_;
+
             //- Continuous phase time at start of injection time step [s]
             scalar time0_;
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.H
index 7d8305088f3f92f01eafbe171a6e565605e845ff..0202c9eb60be1a3c10324b67fd057208cdc1aa42 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.H
@@ -76,9 +76,6 @@ class ManualInjection
         //- Parcel size PDF model
         const autoPtr<pdfs::pdf> parcelPDF_;
 
-        //- Number of particles represented by each parcel
-        scalar nParticlesPerParcel_;
-
 
 protected:
 
diff --git a/src/turbulenceModels/incompressible/RAS/Make/files b/src/turbulenceModels/incompressible/RAS/Make/files
index 7ddd1d4e3b24fa0d6156d5f6fed061beebf8ace1..f6107822dfbc1cee80f472537f5ef4448d187c6d 100644
--- a/src/turbulenceModels/incompressible/RAS/Make/files
+++ b/src/turbulenceModels/incompressible/RAS/Make/files
@@ -43,6 +43,7 @@ $(kqRWallFunctions)/kqRWallFunction/kqRWallFunctionFvPatchFields.C
 derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C
 derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
 derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
+derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C
 
 backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C
new file mode 100644
index 0000000000000000000000000000000000000000..8291b745f8ca9ece50b8fecc160aeed1476d269e
--- /dev/null
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C
@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 "fixedShearStressFvPatchVectorField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "RASModel.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+fixedShearStressFvPatchVectorField::
+fixedShearStressFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(p, iF),
+    tau0_(vector::zero)
+{}
+
+
+fixedShearStressFvPatchVectorField::
+fixedShearStressFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchVectorField(p, iF),
+    tau0_(dict.lookupOrDefault<vector>("tau", vector::zero))
+{
+    fvPatchField<vector>::operator=(patchInternalField());
+}
+
+
+fixedShearStressFvPatchVectorField::
+fixedShearStressFvPatchVectorField
+(
+    const fixedShearStressFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchVectorField(ptf, p, iF, mapper),
+    tau0_(ptf.tau0_)
+{}
+
+
+fixedShearStressFvPatchVectorField::
+fixedShearStressFvPatchVectorField
+(
+    const fixedShearStressFvPatchVectorField& ptf
+)
+:
+    fixedValueFvPatchVectorField(ptf),
+    tau0_(ptf.tau0_)
+{}
+
+
+fixedShearStressFvPatchVectorField::
+fixedShearStressFvPatchVectorField
+(
+    const fixedShearStressFvPatchVectorField& ptf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(ptf, iF),
+    tau0_(ptf.tau0_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void fixedShearStressFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+
+    const vectorField Ui = Uw.patchInternalField();
+
+    vector tauHat = tau0_/mag(tau0_);
+
+    const scalarField& ry = patch().deltaCoeffs();
+
+    scalarField nuEffw = rasModel.nuEff()().boundaryField()[patchI];
+
+    vectorField UwUpdated =
+        tauHat*(tauHat & (tau0_*(1.0/(ry*nuEffw)) + Ui));
+
+    operator==(UwUpdated);
+
+    if (debug)
+    {
+        vectorField nHat = this->patch().nf();
+        volSymmTensorField Reff = rasModel.devReff();
+        Info<< "tau : " << (nHat & Reff.boundaryField()[patchI])() << endl;
+    }
+
+    fixedValueFvPatchVectorField::updateCoeffs();
+}
+
+
+void fixedShearStressFvPatchVectorField::write(Ostream& os) const
+{
+    fixedValueFvPatchVectorField::write(os);
+    os.writeKeyword("tau") << tau0_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchVectorField,
+    fixedShearStressFvPatchVectorField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+// ************************************************************************* //
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.H
new file mode 100644
index 0000000000000000000000000000000000000000..a8878cc450e9c4ce8355a55c1495a10047321166
--- /dev/null
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.H
@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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::fixedShearStressFvPatchVectorField
+
+Description
+    Set a constant shear stress as tau0 = -nuEff dU/dn.
+
+SourceFiles
+    fixedShearStressFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedShearStressFvPatchVectorField_H
+#define fixedShearStressFvPatchVectorField_H
+
+#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+/*---------------------------------------------------------------------------*\
+             Class fixedShearStressFvPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class fixedShearStressFvPatchVectorField
+:
+    public fixedValueFvPatchVectorField
+{
+    // Private data
+
+        //- Constant shear stress
+        vector tau0_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fixedShearStress");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        fixedShearStressFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        fixedShearStressFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        fixedShearStressFvPatchVectorField
+        (
+            const fixedShearStressFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        fixedShearStressFvPatchVectorField
+        (
+            const fixedShearStressFvPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new fixedShearStressFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        fixedShearStressFvPatchVectorField
+        (
+            const fixedShearStressFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
+        (
+            const DimensionedField<vector, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new fixedShearStressFvPatchVectorField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //