diff --git a/applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.C b/applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.C
index ee4544ee34529b4b397beed5c8a012ef67437fd7..712ed6564e972dfc28a875731472ea045763d527 100644
--- a/applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.C
+++ b/applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.C
@@ -51,6 +51,7 @@ Description
 #include "uindirectPrimitivePatch.H"
 #include "globalMeshData.H"
 #include "globalIndex.H"
+#include "timeSelector.H"
 
 using namespace Foam;
 
@@ -63,6 +64,8 @@ int main(int argc, char *argv[])
     (
         "extract surface from a polyMesh"
     );
+    timeSelector::addOptions();
+
     argList::validArgs.append("output file");
     #include "addRegionOption.H"
     argList::addBoolOption
@@ -86,7 +89,14 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
 
-    const fileName outFileName(args[1]);
+    const fileName userOutFileName(args[1]);
+
+    if (!userOutFileName.hasExt())
+    {
+        FatalErrorInFunction
+            << "Missing extension on output name " << userOutFileName
+            << exit(FatalError);
+    }
 
     Info<< "Extracting surface from boundaryMesh ..."
         << endl << endl;
@@ -106,275 +116,321 @@ int main(int argc, char *argv[])
         Info<< "Excluding all processor patches." << nl << endl;
     }
 
+
     Info<< "Reading mesh from time " << runTime.value() << endl;
 
     #include "createNamedPolyMesh.H"
 
+    // User specified times
+    instantList timeDirs = timeSelector::select0(runTime, args);
 
-    // Create local surface from:
-    // - explicitly named patches only (-patches (at your option)
-    // - all patches (default in sequential mode)
-    // - all non-processor patches (default in parallel mode)
-    // - all non-processor patches (sequential mode, -excludeProcPatches
-    //   (at your option)
-
-    // Construct table of patches to include.
-    const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
-
-    labelHashSet includePatches(bMesh.size());
-
-    if (args.optionFound("patches"))
+    forAll(timeDirs, timeIndex)
     {
-        includePatches = bMesh.patchSet
-        (
-            wordReList(args.optionLookup("patches")())
-        );
-    }
-    else
-    {
-        forAll(bMesh, patchi)
-        {
-            const polyPatch& patch = bMesh[patchi];
+        runTime.setTime(timeDirs[timeIndex], timeIndex);
+        Info<< "Time [" << timeIndex << "] = " << runTime.timeName();
 
-            if (includeProcPatches || !isA<processorPolyPatch>(patch))
-            {
-                includePatches.insert(patchi);
-            }
+        fileName outFileName;
+        if (timeDirs.size() == 1)
+        {
+            outFileName = userOutFileName;
+            Info<< nl;
         }
-    }
-
-
-
-    const faceZoneMesh& fzm = mesh.faceZones();
-    labelHashSet includeFaceZones(fzm.size());
-
-    if (args.optionFound("faceZones"))
-    {
-        wordReList zoneNames(args.optionLookup("faceZones")());
-        const wordList allZoneNames(fzm.names());
-        forAll(zoneNames, i)
+        else
         {
-            const wordRe& zoneName = zoneNames[i];
-
-            labelList zoneIDs = findStrings(zoneName, allZoneNames);
-
-            forAll(zoneIDs, j)
+            polyMesh::readUpdateState meshState = mesh.readUpdate();
+            if (timeIndex && meshState == polyMesh::UNCHANGED)
             {
-                includeFaceZones.insert(zoneIDs[j]);
+                Info<<"  ... no mesh change." << nl;
+                continue;
             }
-
-            if (zoneIDs.empty())
+            else
             {
-                WarningInFunction
-                    << "Cannot find any faceZone name matching "
-                    << zoneName << endl;
+                Info<< nl;
             }
 
+            // The filename based on the original, but with additional
+            // time information. The extension was previously checked that
+            // it exists
+            std::string::size_type dot = userOutFileName.rfind('.');
+
+            outFileName =
+                userOutFileName.substr(0, dot) + "_"
+              + Foam::name(runTime.value()) + "."
+              + userOutFileName.ext();
         }
-        Info<< "Additionally triangulating faceZones "
-            << UIndirectList<word>(allZoneNames, includeFaceZones.sortedToc())
-            << endl;
-    }
 
+        // Create local surface from:
+        // - explicitly named patches only (-patches (at your option)
+        // - all patches (default in sequential mode)
+        // - all non-processor patches (default in parallel mode)
+        // - all non-processor patches (sequential mode, -excludeProcPatches
+        //   (at your option)
 
+        // Construct table of patches to include.
+        const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
 
-    // From (name of) patch to compact 'zone' index
-    HashTable<label> compactZoneID(1000);
-    // Mesh face and compact zone indx
-    DynamicList<label> faceLabels;
-    DynamicList<label> compactZones;
+        labelHashSet includePatches(bMesh.size());
 
-    {
-        // Collect sizes. Hash on names to handle local-only patches (e.g.
-        //  processor patches)
-        HashTable<label> patchSize(1000);
-        label nFaces = 0;
-        forAllConstIter(labelHashSet, includePatches, iter)
+        if (args.optionFound("patches"))
         {
-            const polyPatch& pp = bMesh[iter.key()];
-            patchSize.insert(pp.name(), pp.size());
-            nFaces += pp.size();
+            includePatches = bMesh.patchSet
+            (
+                wordReList(args.optionLookup("patches")())
+            );
         }
-
-        HashTable<label> zoneSize(1000);
-        forAllConstIter(labelHashSet, includeFaceZones, iter)
+        else
         {
-            const faceZone& pp = fzm[iter.key()];
-            zoneSize.insert(pp.name(), pp.size());
-            nFaces += pp.size();
-        }
+            forAll(bMesh, patchi)
+            {
+                const polyPatch& patch = bMesh[patchi];
 
+                if (includeProcPatches || !isA<processorPolyPatch>(patch))
+                {
+                    includePatches.insert(patchi);
+                }
+            }
+        }
 
-        Pstream::mapCombineGather(patchSize, plusEqOp<label>());
-        Pstream::mapCombineGather(zoneSize, plusEqOp<label>());
 
+        const faceZoneMesh& fzm = mesh.faceZones();
+        labelHashSet includeFaceZones(fzm.size());
 
-        // Allocate compact numbering for all patches/faceZones
-        forAllConstIter(HashTable<label>, patchSize, iter)
+        if (args.optionFound("faceZones"))
         {
-            label sz = compactZoneID.size();
-            compactZoneID.insert(iter.key(), sz);
-        }
+            wordReList zoneNames(args.optionLookup("faceZones")());
+            const wordList allZoneNames(fzm.names());
+            forAll(zoneNames, i)
+            {
+                const wordRe& zoneName = zoneNames[i];
 
-        forAllConstIter(HashTable<label>, zoneSize, iter)
-        {
-            label sz = compactZoneID.size();
-            //Info<< "For faceZone " << iter.key() << " allocating zoneID "
-            //    << sz << endl;
-            compactZoneID.insert(iter.key(), sz);
+                labelList zoneIDs = findStrings(zoneName, allZoneNames);
+
+                forAll(zoneIDs, j)
+                {
+                    includeFaceZones.insert(zoneIDs[j]);
+                }
+
+                if (zoneIDs.empty())
+                {
+                    WarningInFunction
+                        << "Cannot find any faceZone name matching "
+                        << zoneName << endl;
+                }
+
+            }
+            Info<< "Additionally triangulating faceZones "
+                << UIndirectList<word>
+                  (
+                      allZoneNames,
+                      includeFaceZones.sortedToc()
+                  )
+                << endl;
         }
 
 
-        Pstream::mapCombineScatter(compactZoneID);
 
+        // From (name of) patch to compact 'zone' index
+        HashTable<label> compactZoneID(1000);
+        // Mesh face and compact zone indx
+        DynamicList<label> faceLabels;
+        DynamicList<label> compactZones;
 
-        // Rework HashTable into labelList just for speed of conversion
-        labelList patchToCompactZone(bMesh.size(), -1);
-        labelList faceZoneToCompactZone(bMesh.size(), -1);
-        forAllConstIter(HashTable<label>, compactZoneID, iter)
         {
-            label patchi = bMesh.findPatchID(iter.key());
-            if (patchi != -1)
+            // Collect sizes. Hash on names to handle local-only patches (e.g.
+            //  processor patches)
+            HashTable<label> patchSize(1000);
+            label nFaces = 0;
+            forAllConstIter(labelHashSet, includePatches, iter)
             {
-                patchToCompactZone[patchi] = iter();
+                const polyPatch& pp = bMesh[iter.key()];
+                patchSize.insert(pp.name(), pp.size());
+                nFaces += pp.size();
             }
-            else
+
+            HashTable<label> zoneSize(1000);
+            forAllConstIter(labelHashSet, includeFaceZones, iter)
             {
-                label zoneI = fzm.findZoneID(iter.key());
-                faceZoneToCompactZone[zoneI] = iter();
+                const faceZone& pp = fzm[iter.key()];
+                zoneSize.insert(pp.name(), pp.size());
+                nFaces += pp.size();
             }
-        }
 
 
-        faceLabels.setCapacity(nFaces);
-        compactZones.setCapacity(nFaces);
+            Pstream::mapCombineGather(patchSize, plusEqOp<label>());
+            Pstream::mapCombineGather(zoneSize, plusEqOp<label>());
 
-        // Collect faces on patches
-        forAllConstIter(labelHashSet, includePatches, iter)
-        {
-            const polyPatch& pp = bMesh[iter.key()];
-            forAll(pp, i)
+
+            // Allocate compact numbering for all patches/faceZones
+            forAllConstIter(HashTable<label>, patchSize, iter)
             {
-                faceLabels.append(pp.start()+i);
-                compactZones.append(patchToCompactZone[pp.index()]);
+                label sz = compactZoneID.size();
+                compactZoneID.insert(iter.key(), sz);
             }
-        }
-        // Collect faces on faceZones
-        forAllConstIter(labelHashSet, includeFaceZones, iter)
-        {
-            const faceZone& pp = fzm[iter.key()];
-            forAll(pp, i)
+
+            forAllConstIter(HashTable<label>, zoneSize, iter)
             {
-                faceLabels.append(pp[i]);
-                compactZones.append(faceZoneToCompactZone[pp.index()]);
+                label sz = compactZoneID.size();
+                //Info<< "For faceZone " << iter.key() << " allocating zoneID "
+                //    << sz << endl;
+                compactZoneID.insert(iter.key(), sz);
             }
-        }
-    }
 
 
-    // Addressing engine for all faces
-    uindirectPrimitivePatch allBoundary
-    (
-        UIndirectList<face>(mesh.faces(), faceLabels),
-        mesh.points()
-    );
+            Pstream::mapCombineScatter(compactZoneID);
 
 
-    // Find correspondence to master points
-    labelList pointToGlobal;
-    labelList uniqueMeshPoints;
-    autoPtr<globalIndex> globalNumbers = mesh.globalData().mergePoints
-    (
-        allBoundary.meshPoints(),
-        allBoundary.meshPointMap(),
-        pointToGlobal,
-        uniqueMeshPoints
-    );
+            // Rework HashTable into labelList just for speed of conversion
+            labelList patchToCompactZone(bMesh.size(), -1);
+            labelList faceZoneToCompactZone(bMesh.size(), -1);
+            forAllConstIter(HashTable<label>, compactZoneID, iter)
+            {
+                label patchi = bMesh.findPatchID(iter.key());
+                if (patchi != -1)
+                {
+                    patchToCompactZone[patchi] = iter();
+                }
+                else
+                {
+                    label zoneI = fzm.findZoneID(iter.key());
+                    faceZoneToCompactZone[zoneI] = iter();
+                }
+            }
 
-    // Gather all unique points on master
-    List<pointField> gatheredPoints(Pstream::nProcs());
-    gatheredPoints[Pstream::myProcNo()] = pointField
-    (
-        mesh.points(),
-        uniqueMeshPoints
-    );
-    Pstream::gatherList(gatheredPoints);
 
-    // Gather all faces
-    List<faceList> gatheredFaces(Pstream::nProcs());
-    gatheredFaces[Pstream::myProcNo()] = allBoundary.localFaces();
-    forAll(gatheredFaces[Pstream::myProcNo()], i)
-    {
-        inplaceRenumber(pointToGlobal, gatheredFaces[Pstream::myProcNo()][i]);
-    }
-    Pstream::gatherList(gatheredFaces);
+            faceLabels.setCapacity(nFaces);
+            compactZones.setCapacity(nFaces);
 
-    // Gather all ZoneIDs
-    List<labelList> gatheredZones(Pstream::nProcs());
-    gatheredZones[Pstream::myProcNo()] = compactZones.xfer();
-    Pstream::gatherList(gatheredZones);
+            // Collect faces on patches
+            forAllConstIter(labelHashSet, includePatches, iter)
+            {
+                const polyPatch& pp = bMesh[iter.key()];
+                forAll(pp, i)
+                {
+                    faceLabels.append(pp.start()+i);
+                    compactZones.append(patchToCompactZone[pp.index()]);
+                }
+            }
+            // Collect faces on faceZones
+            forAllConstIter(labelHashSet, includeFaceZones, iter)
+            {
+                const faceZone& pp = fzm[iter.key()];
+                forAll(pp, i)
+                {
+                    faceLabels.append(pp[i]);
+                    compactZones.append(faceZoneToCompactZone[pp.index()]);
+                }
+            }
+        }
 
-    // On master combine all points, faces, zones
-    if (Pstream::master())
-    {
-        pointField allPoints = ListListOps::combine<pointField>
+
+        // Addressing engine for all faces
+        uindirectPrimitivePatch allBoundary
         (
-            gatheredPoints,
-            accessOp<pointField>()
+            UIndirectList<face>(mesh.faces(), faceLabels),
+            mesh.points()
         );
-        gatheredPoints.clear();
 
-        faceList allFaces = ListListOps::combine<faceList>
+
+        // Find correspondence to master points
+        labelList pointToGlobal;
+        labelList uniqueMeshPoints;
+        autoPtr<globalIndex> globalNumbers = mesh.globalData().mergePoints
         (
-            gatheredFaces,
-            accessOp<faceList>()
+            allBoundary.meshPoints(),
+            allBoundary.meshPointMap(),
+            pointToGlobal,
+            uniqueMeshPoints
         );
-        gatheredFaces.clear();
 
-        labelList allZones = ListListOps::combine<labelList>
+        // Gather all unique points on master
+        List<pointField> gatheredPoints(Pstream::nProcs());
+        gatheredPoints[Pstream::myProcNo()] = pointField
         (
-            gatheredZones,
-            accessOp<labelList>()
+            mesh.points(),
+            uniqueMeshPoints
         );
-        gatheredZones.clear();
-
+        Pstream::gatherList(gatheredPoints);
 
-        // Zones
-        surfZoneIdentifierList surfZones(compactZoneID.size());
-        forAllConstIter(HashTable<label>, compactZoneID, iter)
+        // Gather all faces
+        List<faceList> gatheredFaces(Pstream::nProcs());
+        gatheredFaces[Pstream::myProcNo()] = allBoundary.localFaces();
+        forAll(gatheredFaces[Pstream::myProcNo()], i)
         {
-            surfZones[iter()] = surfZoneIdentifier(iter.key(), iter());
-            Info<< "surfZone " << iter()  <<  " : " << surfZones[iter()].name()
-                << endl;
+            inplaceRenumber
+           (
+                pointToGlobal,
+                gatheredFaces[Pstream::myProcNo()][i]
+           );
         }
+        Pstream::gatherList(gatheredFaces);
 
-        UnsortedMeshedSurface<face> unsortedFace
-        (
-            xferMove(allPoints),
-            xferMove(allFaces),
-            xferMove(allZones),
-            xferMove(surfZones)
-        );
+        // Gather all ZoneIDs
+        List<labelList> gatheredZones(Pstream::nProcs());
+        gatheredZones[Pstream::myProcNo()] = compactZones.xfer();
+        Pstream::gatherList(gatheredZones);
 
+        // On master combine all points, faces, zones
+        if (Pstream::master())
+        {
+            pointField allPoints = ListListOps::combine<pointField>
+            (
+                gatheredPoints,
+                accessOp<pointField>()
+            );
+            gatheredPoints.clear();
+
+            faceList allFaces = ListListOps::combine<faceList>
+            (
+                gatheredFaces,
+                accessOp<faceList>()
+            );
+            gatheredFaces.clear();
+
+            labelList allZones = ListListOps::combine<labelList>
+            (
+                gatheredZones,
+                accessOp<labelList>()
+            );
+            gatheredZones.clear();
+
+
+            // Zones
+            surfZoneIdentifierList surfZones(compactZoneID.size());
+            forAllConstIter(HashTable<label>, compactZoneID, iter)
+            {
+                surfZones[iter()] = surfZoneIdentifier(iter.key(), iter());
+                Info<< "surfZone " << iter()
+                    <<  " : "      << surfZones[iter()].name()
+                    << endl;
+            }
 
-        MeshedSurface<face> sortedFace(unsortedFace);
+            UnsortedMeshedSurface<face> unsortedFace
+            (
+                xferMove(allPoints),
+                xferMove(allFaces),
+                xferMove(allZones),
+                xferMove(surfZones)
+            );
 
-        fileName globalCasePath
-        (
-            outFileName.isAbsolute()
-          ? outFileName
-          : (
-                runTime.processorCase()
-              ? runTime.rootPath()/runTime.globalCaseName()/outFileName
-              : runTime.path()/outFileName
-            )
-        );
 
-        Info<< "Writing merged surface to " << globalCasePath << endl;
+            MeshedSurface<face> sortedFace(unsortedFace);
 
-        sortedFace.write(globalCasePath);
-    }
+            fileName globalCasePath
+            (
+                outFileName.isAbsolute()
+              ? outFileName
+              : (
+                    runTime.processorCase()
+                  ? runTime.rootPath()/runTime.globalCaseName()/outFileName
+                  : runTime.path()/outFileName
+                )
+            );
+
+            Info<< "Writing merged surface to " << globalCasePath << endl;
 
+            sortedFace.write(globalCasePath);
+        }
+    }
     Info<< "End\n" << endl;
 
     return 0;