diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLine.C b/src/postProcessing/functionObjects/field/streamLine/streamLine.C
index 30d8052513..ab815a8df8 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLine.C
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLine.C
@@ -34,6 +34,7 @@ License
 #include "globalIndex.H"
 #include "mapDistribute.H"
 #include "interpolationCellPoint.H"
+#include "PatchTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -41,6 +42,57 @@ defineTypeNameAndDebug(Foam::streamLine, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+Foam::autoPtr<Foam::indirectPrimitivePatch>
+Foam::streamLine::wallPatch() const
+{
+    const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    label nFaces = 0;
+
+    forAll(patches, patchI)
+    {
+        //if (!polyPatch::constraintType(patches[patchI].type()))
+        if (isA<wallPolyPatch>(patches[patchI]))
+        {
+            nFaces += patches[patchI].size();
+        }
+    }
+
+    labelList addressing(nFaces);
+
+    nFaces = 0;
+
+    forAll(patches, patchI)
+    {
+        //if (!polyPatch::constraintType(patches[patchI].type()))
+        if (isA<wallPolyPatch>(patches[patchI]))
+        {
+            const polyPatch& pp = patches[patchI];
+
+            forAll(pp, i)
+            {
+                addressing[nFaces++] = pp.start()+i;
+            }
+        }
+    }
+
+    return autoPtr<indirectPrimitivePatch>
+    (
+        new indirectPrimitivePatch
+        (
+            IndirectList<face>
+            (
+                mesh.faces(),
+                addressing
+            ),
+            mesh.points()
+        )
+    );
+}
+
+
 void Foam::streamLine::track()
 {
     const Time& runTime = obr_.time();
@@ -72,7 +124,7 @@ void Foam::streamLine::track()
 
     label nSeeds = returnReduce(particles.size(), sumOp<label>());
 
-    Info<< "Seeded " << nSeeds << " particles." << endl;
+    Info<< type() << " : seeded " << nSeeds << " particles." << endl;
 
     // Read or lookup fields
     PtrList<volScalarField> vsFlds;
@@ -163,7 +215,6 @@ void Foam::streamLine::track()
                 vsInterp.set
                 (
                     nScalar++,
-                    //new interpolationCellPoint<scalar>(f)
                     interpolation<scalar>::New
                     (
                         interpolationScheme_,
@@ -186,7 +237,6 @@ void Foam::streamLine::track()
                 vvInterp.set
                 (
                     nVector++,
-                    //new interpolationCellPoint<vector>(f)
                     interpolation<vector>::New
                     (
                         interpolationScheme_,
@@ -241,6 +291,120 @@ void Foam::streamLine::track()
         allVectors_[i].setCapacity(nSeeds);
     }
 
+
+    //- For surface following: per face the normal to constrain the velocity
+    //  with.
+    pointField patchNormals;
+
+
+    if (followSurface_)
+    {
+        // Determine geometric data on the non-constraint patches
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // - point normals
+        // - neighbouring cells
+        autoPtr<indirectPrimitivePatch> boundaryPatch(wallPatch());
+
+        // Calculate parallel consistent normal (on boundaryPatch points only)
+        pointField boundaryNormals
+        (
+            PatchTools::pointNormals
+            (
+                mesh,
+                boundaryPatch(),
+                boundaryPatch().addressing()
+            )
+        );
+        // Make sure all points know about point normal. There might be points
+        // which are not on boundaryPatch() but coupled to points that are.
+        pointField pointNormals(mesh.nPoints(), vector::zero);
+        UIndirectList<point>
+        (
+            pointNormals,
+            boundaryPatch().meshPoints()
+        ) = boundaryNormals;
+        boundaryNormals.clear();
+
+        mesh.globalData().syncPointData
+        (
+            pointNormals,
+            maxMagSqrEqOp<point>(),
+            mapDistribute::transform()
+        );
+
+        // Calculate per-face a single constraint normal:
+        // - faces on patch: face-normal
+        // - faces on cell with point on patch: point-normal
+        // - rest: vector::zero
+        //
+        // Note: the disadvantage is that we can only handle faces with one
+        //       point on the patch - otherwise we might pick up the wrong
+        //       velocity.
+
+
+
+        patchNormals.setSize(mesh.nFaces(), vector::zero);
+        label nPointNormals = 0;
+
+        // 1. faces on patch
+        UIndirectList<point>(patchNormals, boundaryPatch().addressing()) =
+             boundaryPatch().faceNormals();
+
+        // 2. faces with a point on the patch
+        forAll(mesh.faces(), faceI)
+        {
+            if (patchNormals[faceI] == vector::zero)
+            {
+                const face& f = mesh.faces()[faceI];
+                forAll(f, fp)
+                {
+                    label pointI = f[fp];
+                    if (pointNormals[pointI] != vector::zero)
+                    {
+                        patchNormals[faceI] = pointNormals[pointI];
+                        nPointNormals++;
+                    }
+                }
+            }
+        }
+
+        // 3. cells on faces with a point on the patch
+        forAll(mesh.points(), pointI)
+        {
+            if (pointNormals[pointI] != vector::zero)
+            {
+                const labelList& pCells = mesh.pointCells(pointI);
+                forAll(pCells, i)
+                {
+                    const cell& cFaces = mesh.cells()[pCells[i]];
+                    forAll(cFaces, j)
+                    {
+                        label faceI = cFaces[j];
+
+                        if (patchNormals[faceI] == vector::zero)
+                        {
+                            patchNormals[faceI] = pointNormals[pointI];
+                            nPointNormals++;
+                        }
+                    }
+                }
+            }
+        }
+
+        Info<< type() << " : calculated constrained-tracking normals" << nl
+            << "Number of faces in mesh                     : "
+            << returnReduce(mesh.nFaces(), sumOp<label>()) << nl
+            << "    of which on walls                       : "
+            <<  returnReduce
+                (
+                    boundaryPatch().addressing().size(),
+                    sumOp<label>()
+                ) << nl
+            << "    of which on cell with point on boundary : "
+            << returnReduce(nPointNormals, sumOp<label>()) << nl
+            << endl;
+    }
+
     // additional particle info
     streamLineParticle::trackingData td
     (
@@ -248,8 +412,16 @@ void Foam::streamLine::track()
         vsInterp,
         vvInterp,
         UIndex,         // index of U in vvInterp
-        trackForward_,  // track in +u direction
-        nSubCycle_,     // step through cells in steps?
+        trackForward_,  // track in +u direction?
+        nSubCycle_,     // automatic track control:step through cells in steps?
+        trackLength_,   // fixed track length
+
+        //- For surface-constrained streamlines
+        followSurface_, // stay on surface?
+        patchNormals,   // per face normal to constrain tracking to
+        //pointNormals,   // per point          ,,
+        trackHeight_,   // track height
+
         allTracks_,
         allScalars_,
         allVectors_
@@ -279,7 +451,8 @@ Foam::streamLine::streamLine
     name_(name),
     obr_(obr),
     loadFromFiles_(loadFromFiles),
-    active_(true)
+    active_(true),
+    nSubCycle_(0)
 {
     // Only active if a fvMesh is available
     if (isA<fvMesh>(obr_))
@@ -334,6 +507,16 @@ void Foam::streamLine::read(const dictionary& dict)
                 dict.lookup("U") >> UName_;
             }
         }
+
+        if (findIndex(fields_, UName_) == -1)
+        {
+            FatalIOErrorIn("streamLine::read(const dictionary&)", dict)
+                << "Velocity field for tracking " << UName_
+                << " should be present in the list of fields " << fields_
+                << exit(FatalIOError);
+        }
+
+
         dict.lookup("trackForward") >> trackForward_;
         dict.lookup("lifeTime") >> lifeTime_;
         if (lifeTime_ < 1)
@@ -343,10 +526,48 @@ void Foam::streamLine::read(const dictionary& dict)
                 << exit(FatalError);
         }
 
-        dict.lookup("nSubCycle") >> nSubCycle_;
-        if (nSubCycle_ < 1)
+
+        bool subCycling = dict.found("nSubCycle");
+        bool fixedLength = dict.found("trackLength");
+
+        if (subCycling && fixedLength)
         {
-            nSubCycle_ = 1;
+            FatalIOErrorIn("streamLine::read(const dictionary&)", dict)
+                << "Cannot both specify automatic time stepping (through '"
+                << "nSubCycle' specification) and fixed track length (through '"
+                << "trackLength')"
+                << exit(FatalIOError);
+        }
+
+
+        nSubCycle_ = 1;
+        if (dict.readIfPresent("nSubCycle", nSubCycle_))
+        {
+            trackLength_ = VGREAT;
+            if (nSubCycle_ < 1)
+            {
+                nSubCycle_ = 1;
+            }
+            Info<< type() << " : automatic track length specified through"
+                << " number of sub cycles : " << nSubCycle_ << nl << endl;
+        }
+        else
+        {
+            dict.lookup("trackLength") >> trackLength_;
+
+            Info<< type() << " : fixed track length specified : "
+                << trackLength_ << nl << endl;
+        }
+
+
+        followSurface_ = false;
+        trackHeight_ = VGREAT;
+        if (dict.readIfPresent("followSurface", followSurface_))
+        {
+            dict.lookup("trackHeight") >> trackHeight_;
+
+            Info<< type() << " : surface-constrained streamline at "
+                << trackHeight_ << " distance above the surface" << nl << endl;
         }
 
         interpolationScheme_ = dict.lookupOrDefault
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLine.H b/src/postProcessing/functionObjects/field/streamLine/streamLine.H
index 9fc981034c..77984ed2c8 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLine.H
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLine.H
@@ -43,6 +43,7 @@ SourceFiles
 #include "vectorList.H"
 #include "polyMesh.H"
 #include "writer.H"
+#include "indirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -95,9 +96,17 @@ class streamLine
         //- Maximum lifetime (= number of cells) of particle
         label lifeTime_;
 
+        //- Surface-constrained?
+        bool followSurface_;
+
+        scalar trackHeight_;
+
         //- Number of subcycling steps
         label nSubCycle_;
 
+        //- Track length
+        scalar trackLength_;
+
         //- Optional specified name of particles
         word cloudName_;
 
@@ -141,6 +150,8 @@ class streamLine
             List<DynamicList<vectorList> > allVectors_;
 
 
+        //- Construct patch out of all wall patch faces
+        autoPtr<indirectPrimitivePatch> wallPatch() const;
 
         //- Do all seeding and tracking
         void track();
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
index c7808b9900..279e17ef80 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
@@ -43,7 +43,7 @@ Foam::scalar Foam::streamLineParticle::calcSubCycleDeltaT
     const vector& U
 ) const
 {
-    streamLineParticle testParticle(*this);
+    particle testParticle(*this);
 
     bool oldKeepParticle = td.keepParticle;
     bool oldSwitchProcessor = td.switchProcessor;
@@ -59,7 +59,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
 (
     const trackingData& td,
     const point& position,
-    const label cellI
+    const label cellI,
+    const label faceI
 )
 {
     if (cellI == -1)
@@ -76,7 +77,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
             td.vsInterp_[scalarI].interpolate
             (
                 position,
-                cellI
+                cellI,
+                faceI
             )
         );
     }
@@ -89,7 +91,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
             td.vvInterp_[vectorI].interpolate
             (
                 position,
-                cellI
+                cellI,
+                faceI
             )
         );
     }
@@ -201,13 +204,24 @@ bool Foam::streamLineParticle::move
 
             // Store current position and sampled velocity.
             sampledPositions_.append(position());
-            vector U = interpolateFields(td, position(), cell());
+            vector U = interpolateFields(td, position(), cell(), face());
 
             if (!td.trackForward_)
             {
                 U = -U;
             }
 
+            if (td.followSurface_)
+            {
+                //- Simple: remove wall-normal component. Should keep particle
+                //  at same height. Does not work well on warped faces.
+                const vector& n = td.patchNormals_[tetFaceI_];
+                if (n != vector::zero)
+                {
+                    U -= (U&n)*n;
+                }
+            }
+
             scalar magU = mag(U);
 
             if (magU < SMALL)
@@ -219,7 +233,13 @@ bool Foam::streamLineParticle::move
 
             U /= magU;
 
-            if (subIter == 0 && td.nSubCycle_ > 1)
+            if (td.trackLength_ < GREAT)
+            {
+                dt = td.trackLength_;
+                //Pout<< "    subiteration " << subIter
+                //    << " : fixed length: updated dt:" << dt << endl;
+            }
+            else if (subIter == 0 && td.nSubCycle_ > 1)
             {
                 // Adapt dt to cross cell in a few steps
                 dt = calcSubCycleDeltaT(td, dt, U);
@@ -264,7 +284,8 @@ bool Foam::streamLineParticle::move
             if (debug)
             {
                 Pout<< "streamLineParticle : Removing stagnant particle:"
-                    << p << " sampled positions:" << sampledPositions_.size()
+                    << p.position()
+                    << " sampled positions:" << sampledPositions_.size()
                     << endl;
             }
             td.keepParticle = false;
@@ -273,12 +294,13 @@ bool Foam::streamLineParticle::move
         {
             // Normal exit. Store last position and fields
             sampledPositions_.append(position());
-            interpolateFields(td, position(), cell());
+            interpolateFields(td, position(), cell(), face());
 
             if (debug)
             {
                 Pout<< "streamLineParticle : Removing particle:"
-                    << p << " sampled positions:" << sampledPositions_.size()
+                    << p.position()
+                    << " sampled positions:" << sampledPositions_.size()
                     << endl;
             }
         }
@@ -374,8 +396,42 @@ void Foam::streamLineParticle::hitWallPatch
     const tetIndices&
 )
 {
-    // Remove particle
-    td.keepParticle = false;
+    if (td.followSurface_)
+    {
+        // Push slightly in. Since we're pushing to tet centre we should still
+        // be in the same tet.
+        tetPointRef tet = currentTet();
+
+        const point oldPosition(position());
+
+        position() += trackingCorrectionTol*(tet.centre() - position());
+
+        if (debug)
+        {
+            label patchID = patch(face());
+
+            Pout<< "hitWallPatch : on wall "
+                << mesh().boundaryMesh()[patchID].name()
+                << " moved from " << oldPosition
+                << " to " << position()
+                << " due to centre " << tet.centre() << endl;
+        }
+
+        // Check that still in tet (should always be the case
+        if (debug && !tet.inside(position()))
+        {
+            FatalErrorIn("streamLineParticle::hitWallPatch()")
+                << "Pushing particle " << *this
+                << " away from wall " << wpp.name()
+                << " moved it outside of tet."
+                << exit(FatalError);
+        }
+    }
+    else
+    {
+        // Remove particle
+        td.keepParticle = false;
+    }
 }
 
 
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
index b1ee9acf88..c5feff4285 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
@@ -73,6 +73,13 @@ public:
         const label UIndex_;
         const bool trackForward_;
         const label nSubCycle_;
+        const scalar trackLength_;
+
+        //- For surface-constrained tracking
+        const bool followSurface_;
+        const pointField& patchNormals_;
+        const scalar trackHeight_;
+
 
         DynamicList<vectorList>& allPositions_;
         List<DynamicList<scalarList> >& allScalars_;
@@ -89,6 +96,12 @@ public:
                 const label UIndex,
                 const bool trackForward,
                 const label nSubCycle,
+                const scalar trackLength,
+                const bool followSurface,
+
+                const pointField& patchNormals,
+                const scalar trackHeight,
+
                 DynamicList<List<point> >& allPositions,
                 List<DynamicList<scalarList> >& allScalars,
                 List<DynamicList<vectorList> >& allVectors
@@ -100,6 +113,12 @@ public:
                 UIndex_(UIndex),
                 trackForward_(trackForward),
                 nSubCycle_(nSubCycle),
+                trackLength_(trackLength),
+                followSurface_(followSurface),
+
+                patchNormals_(patchNormals),
+                trackHeight_(trackHeight),
+
                 allPositions_(allPositions),
                 allScalars_(allScalars),
                 allVectors_(allVectors)
@@ -135,12 +154,20 @@ private:
             const vector& U
         ) const;
 
+        void constrainVelocity
+        (
+            trackingData& td,
+            const scalar dt,
+            vector& U
+        );
+
         //- Interpolate all quantities; return interpolated velocity.
         vector interpolateFields
         (
             const trackingData&,
             const point&,
-            const label cellI
+            const label cellI,
+            const label faceI
         );