diff --git a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C
index c68438c0e4..7a8e25a5b6 100644
--- a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C
+++ b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C
@@ -478,6 +478,7 @@ int main(int argc, char *argv[])
                     if
                     (
                         name == "positions"
+                     || name == "coordinates"
                      || name == "origProcId"
                      || name == "origId"
                     )
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
index 8a2c513fa0..004af363fb 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
@@ -745,12 +745,12 @@ int main(int argc, char *argv[])
                         false
                     );
 
-                    IOobject* positionsPtr = sprayObjs.lookup
-                    (
-                        word("positions")
-                    );
+                    // Note: looking up "positions" for backwards compatibility
+                    IOobject* positionsPtr =
+                        sprayObjs.lookup(word("positions"));
+                    IOobject* coordsPtr = sprayObjs.lookup(word("coordinates"));
 
-                    if (positionsPtr)
+                    if (positionsPtr || coordsPtr)
                     {
                         // Read lagrangian particles
                         // ~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
index 716f8522ab..131598746a 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
@@ -586,8 +586,10 @@ int main(int argc, char *argv[])
 
                             IOobject* positionsPtr =
                                 sprayObjs.lookup(word("positions"));
+                            IOobject* coordsPtr =
+                                sprayObjs.lookup(word("coordinates"));
 
-                            if (positionsPtr)
+                            if (coordsPtr || positionsPtr)
                             {
                                 cloudObjects.insert(cloudDirs[i], sprayObjs);
                             }
diff --git a/applications/utilities/parallelProcessing/redistributePar/parLagrangianRedistributor.C b/applications/utilities/parallelProcessing/redistributePar/parLagrangianRedistributor.C
index 9052b8df0b..73cfda6582 100644
--- a/applications/utilities/parallelProcessing/redistributePar/parLagrangianRedistributor.C
+++ b/applications/utilities/parallelProcessing/redistributePar/parLagrangianRedistributor.C
@@ -96,9 +96,14 @@ void Foam::parLagrangianRedistributor::findClouds
             cloud::prefix/localCloudDirs[i]
         );
 
-        if (sprayObjs.lookup(word("positions")))
+        if
+        (
+            sprayObjs.lookup(word("coordinates"))
+         || sprayObjs.lookup(word("positions"))
+        )
         {
-            // One of the objects is positions so must be valid cloud
+            // One of the objects is coordinates/positions so must be valid
+            // cloud
 
             label cloudI = findIndex(cloudNames, localCloudDirs[i]);
 
@@ -107,7 +112,7 @@ void Foam::parLagrangianRedistributor::findClouds
             forAllConstIter(IOobjectList, sprayObjs, iter)
             {
                 const word& name = iter.key();
-                if (name != "positions")
+                if (name != "coordinates" && name != "positions")
                 {
                     objectNames[cloudI][objectI++] = name;
                 }
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/findCloudFields.H b/applications/utilities/postProcessing/dataConversion/foamToEnsight/findCloudFields.H
index 89e1084abc..650258ef1f 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/findCloudFields.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/findCloudFields.H
@@ -35,8 +35,8 @@ if (timeDirs.size() && !noLagrangian)
                 cloudPrefix/cloudName
             );
 
-            // Clouds always have "positions"
-            if (cloudObjs.found("positions"))
+            // Clouds always have "positions" (v1706 and lower) or "coordinates"
+            if (cloudObjs.found("positions") || cloudObjs.found("coordinates"))
             {
                 // Save the cloud fields on a per cloud basis
                 auto fieldsPerCloud = cloudFields(cloudName);
@@ -52,9 +52,10 @@ if (timeDirs.size() && !noLagrangian)
         }
     }
 
-    // Prune out "positions" again since it gets treated specially
+    // Prune out geometry again since it gets treated specially
     forAllIters(cloudFields, cloudIter)
     {
+        cloudIter().erase("coordinates");
         cloudIter().erase("positions");
     }
 
@@ -74,4 +75,5 @@ if (cloudNames.size())
     cloudNames.writeList(Info) << endl;
 }
 
+
 // ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/ensightOutputSerialCloud.C b/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/ensightOutputSerialCloud.C
index 5cc4a45f5c..db87438e46 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/ensightOutputSerialCloud.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/ensightOutputSerialCloud.C
@@ -65,7 +65,7 @@ void Foam::ensightSerialCloud::writePositions
 
             forAllConstIter(Cloud<passiveParticle>, cloudPtr(), elmnt)
             {
-                const vector& p = elmnt().position();
+                const vector p(elmnt().position());
 
                 os.write(p.x());
                 os.write(p.y());
@@ -79,7 +79,7 @@ void Foam::ensightSerialCloud::writePositions
             label parcelId = 0;
             forAllConstIter(Cloud<passiveParticle>, cloudPtr(), elmnt)
             {
-                const vector& p = elmnt().position();
+                const vector p(elmnt().position());
 
                 os.write(++parcelId, 8); // unusual width
                 os.write(p.x());
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/findFields.H b/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/findFields.H
index 66e0616aba..9ac1c56808 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/findFields.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/findFields.H
@@ -60,16 +60,16 @@ if (timeDirs.size())
             cloudPrefix/cloudName
         );
 
-        bool hasPositions = false;
+        bool hasCoordinates = false;
         forAllConstIter(IOobjectList, objs, fieldIter)
         {
             const IOobject obj = *fieldIter();
             const word& fieldName = obj.name();
             const word& fieldType = obj.headerClassName();
 
-            if (fieldName == "positions")
+            if (fieldName == "positions" || fieldName == "coordinates")
             {
-                hasPositions = true;
+                hasCoordinates = true;
             }
             else if (cloudFieldTypes.found(fieldType))
             {
@@ -79,7 +79,7 @@ if (timeDirs.size())
         }
 
         // drop this cloud if it has no positions or is otherwise empty
-        if (!hasPositions || cloudIter().empty())
+        if (!hasCoordinates || cloudIter().empty())
         {
             Info<< "removing cloud " << cloudName << endl;
             cloudFields.erase(cloudIter);
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/foamToEnsightParts.C b/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/foamToEnsightParts.C
index 3aa9ea5855..bac5a91a44 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/foamToEnsightParts.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/foamToEnsightParts.C
@@ -381,8 +381,13 @@ int main(int argc, char *argv[])
                 cloudPrefix/cloudName
             );
 
-            // Check that the positions field is present for this time
-            if (!cloudObjs.found("positions"))
+            // Check that the positions/coordinates field is present for this
+            // time
+            if
+            (
+                !cloudObjs.found("positions")
+             || !cloudObjs.found("coordinates")
+            )
             {
                 continue;
             }
diff --git a/applications/utilities/postProcessing/dataConversion/foamToGMV/gmvOutput.H b/applications/utilities/postProcessing/dataConversion/foamToGMV/gmvOutput.H
index 03a76e5856..a51474d4be 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToGMV/gmvOutput.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToGMV/gmvOutput.H
@@ -62,7 +62,7 @@ for (label i=0; i < nTypes; i++)
 
     if (fieldTypes[i] == cloud::prefix)
     {
-        IOobject lagrangianHeader
+        IOobject positionsHeader
         (
             "positions",
             runTime.timeName(),
@@ -71,9 +71,22 @@ for (label i=0; i < nTypes; i++)
             IOobject::NO_READ
         );
 
+        IOobject coordinatesHeader
+        (
+            "coordinates",
+            runTime.timeName(),
+            cloud::prefix,
+            mesh,
+            IOobject::NO_READ
+        );
+
         if
         (
-            lagrangianHeader.typeHeaderOk<IOPosition<Cloud<passiveParticle>>>
+            positionsHeader.typeHeaderOk<IOPosition<Cloud<passiveParticle>>>
+            (
+                false
+            )
+         || coordinatesHeader.typeHeaderOk<IOPosition<Cloud<passiveParticle>>>
             (
                 false
             )
diff --git a/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C b/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C
index 763fbd4f9d..98f14db918 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C
@@ -1168,8 +1168,9 @@ int main(int argc, char *argv[])
             );
 
             IOobject* positionsPtr = sprayObjs.lookup("positions");
+            IOobject* coordinatesPtr = sprayObjs.lookup("coordinates");
 
-            if (positionsPtr)
+            if (positionsPtr || coordinatesPtr)
             {
                 mkDir(fvPath/cloud::prefix/cloudDirs[cloudI]);
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/findClouds.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/findClouds.H
index 5637f22498..115c4b3d3b 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/findClouds.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/findClouds.H
@@ -23,10 +23,8 @@ if (timeDirs.size() && !noLagrangian)
             fileName::DIRECTORY
         );
 
-        forAll(cloudDirs, cloudI)
+        for (const word& cloudName : cloudDirs)
         {
-            const word& cloudName = cloudDirs[cloudI];
-
             IOobjectList cloudObjs
             (
                 mesh,
@@ -34,8 +32,8 @@ if (timeDirs.size() && !noLagrangian)
                 cloudPrefix/cloudName
             );
 
-            // clouds always require "positions"
-            if (cloudObjs.found("positions"))
+            // Clouds always require "positions"/"coordinates"
+            if (cloudObjs.found("positions") || cloudObjs.found("coordinates"))
             {
                 if (allCloudDirs.insert(cloudName))
                 {
@@ -60,9 +58,9 @@ if (cloudNames.size())
 {
     // complete the echo information
     Info<< "(";
-    forAll(cloudNames, cloudNo)
+    for (const word& cloudName : cloudNames)
     {
-        Info<< ' ' << cloudNames[cloudNo];
+        Info<< ' ' << cloudName;
     }
     Info<< " ) " << endl;
 }
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
index 9575b36373..e84e12f118 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
@@ -1344,7 +1344,7 @@ int main(int argc, char *argv[])
                 cloud::prefix/cloudName
             );
 
-            if (sprayObjs.found("positions"))
+            if (sprayObjs.found("positions") || sprayObjs.found("coordinates"))
             {
                 wordList labelNames(sprayObjs.names(labelIOField::typeName));
                 Info<< "        labels      :";
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMeshLagrangian.C b/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMeshLagrangian.C
index c08e17559c..f7be26d4be 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMeshLagrangian.C
+++ b/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMeshLagrangian.C
@@ -64,7 +64,8 @@ vtkSmartPointer<vtkPolyData> Foam::vtkPVFoam::lagrangianVTKMesh
     );
 
     IOobject* positionsPtr = sprayObjs.lookup(word("positions"));
-    if (positionsPtr)
+    IOobject* coordinatesPtr = sprayObjs.lookup(word("coordinates"));
+    if (positionsPtr || coordinatesPtr)
     {
         Cloud<passiveParticle> parcels(mesh, cloudName, false);
 
diff --git a/applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_set_filenames.H b/applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_set_filenames.H
index ea12f0aef7..a5d8986a75 100644
--- a/applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_set_filenames.H
+++ b/applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_set_filenames.H
@@ -160,7 +160,7 @@ int USERD_set_filenames
         }
     }
 
-    IOobject sprayHeader
+    IOobject positionsHeader
     (
         "positions",
         runTime.timeName(),
@@ -171,7 +171,22 @@ int USERD_set_filenames
         false
     );
 
-    if (sprayHeader.typeHeaderOk<Cloud<passiveParticle>>(false))
+    IOobject coordinatesHeader
+    (
+        "coordinates",
+        runTime.timeName(),
+        cloud::prefix,
+        runTime,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE,
+        false
+    );
+
+    if
+    (
+        positionsHeader.typeHeaderOk<Cloud<passiveParticle>>(false)
+     || coordinatesHeader.typeHeaderOk<Cloud<passiveParticle>>(false)
+    )
     {
         Info<< "[Found lagrangian]" << endl;
 
@@ -181,10 +196,8 @@ int USERD_set_filenames
 
         IOobjectList objects(*meshPtr, runTime.timeName(), cloud::prefix);
 
-        lagrangianScalarNames =
-            objects.names(sprayScalarFieldName);
-        lagrangianVectorNames =
-            objects.names(sprayVectorFieldName);
+        lagrangianScalarNames = objects.names(sprayScalarFieldName);
+        lagrangianVectorNames = objects.names(sprayVectorFieldName);
 
         isSpray[fieldNames.size()] = true;
 
diff --git a/applications/utilities/preProcessing/mapFields/mapLagrangian.C b/applications/utilities/preProcessing/mapFields/mapLagrangian.C
index d450795f47..ebc89520af 100644
--- a/applications/utilities/preProcessing/mapFields/mapLagrangian.C
+++ b/applications/utilities/preProcessing/mapFields/mapLagrangian.C
@@ -120,8 +120,9 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
         );
 
         IOobject* positionsPtr = objects.lookup("positions");
+        IOobject* coordinatesPtr = objects.lookup("coordinates");
 
-        if (positionsPtr)
+        if (positionsPtr || coordinatesPtr)
         {
             Info<< nl << "    processing cloud " << cloudDirs[cloudI] << endl;
 
diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C b/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C
index 2fe1df15ec..56cf1a83ef 100644
--- a/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C
@@ -111,7 +111,10 @@ void mapLagrangian(const meshToMesh& interp)
         bool foundPositions =
             returnReduce(objects.found("positions"), orOp<bool>());;
 
-        if (foundPositions)
+        bool foundCoordinates =
+            returnReduce(objects.found("coordinates"), orOp<bool>());;
+
+        if (foundPositions || foundCoordinates)
         {
             Info<< nl << "    processing cloud " << cloudDirs[cloudI] << endl;
 
diff --git a/etc/controlDict b/etc/controlDict
index 0377378758..77a630c3fb 100644
--- a/etc/controlDict
+++ b/etc/controlDict
@@ -44,6 +44,9 @@ InfoSwitches
     writeDictionaries 0;
     writeOptionalEntries 0;
 
+    // Write lagrangian "positions" file in v1706 format (at earlier)
+    writeLagrangianPositions 0;
+
     // Report list of slaves/pids used (parallel)
     writeSlaves     1;
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
index 7fa09a835d..f2389857d9 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
@@ -63,7 +63,7 @@ Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
 
         const point& tetBasePt = pPts[f[faceBasePtI]];
 
-        for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
+        for (label tetPtI = 1; tetPtI < f.size() - 1; ++tetPtI)
         {
             label facePtI = (tetPtI + faceBasePtI) % f.size();
             label otherFacePtI = f.fcIndex(facePtI);
@@ -158,7 +158,7 @@ Foam::label Foam::polyMeshTetDecomposition::findBasePoint
 
         const point& tetBasePt = pPts[f[faceBasePtI]];
 
-        for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
+        for (label tetPtI = 1; tetPtI < f.size() - 1; ++tetPtI)
         {
             label facePtI = (tetPtI + faceBasePtI) % f.size();
             label otherFacePtI = f.fcIndex(facePtI);
@@ -219,17 +219,16 @@ Foam::labelList Foam::polyMeshTetDecomposition::findFaceBasePts
 
     label nInternalFaces = mesh.nInternalFaces();
 
-    for (label fI = 0; fI < nInternalFaces; fI++)
+    for (label fI = 0; fI < nInternalFaces; ++fI)
     {
         tetBasePtIs[fI] = findSharedBasePoint(mesh, fI, tol, report);
     }
 
     pointField neighbourCellCentres(mesh.nFaces() - nInternalFaces);
 
-    for(label facei = nInternalFaces; facei < mesh.nFaces(); facei++)
+    for (label facei = nInternalFaces; facei < mesh.nFaces(); ++facei)
     {
-        neighbourCellCentres[facei - nInternalFaces] =
-            pC[pOwner[facei]];
+        neighbourCellCentres[facei - nInternalFaces] = pC[pOwner[facei]];
     }
 
     syncTools::swapBoundaryFacePositions(mesh, neighbourCellCentres);
@@ -250,8 +249,7 @@ Foam::labelList Foam::polyMeshTetDecomposition::findFaceBasePts
         fI++, bFI++
     )
     {
-        label patchi =
-            mesh.boundaryMesh().patchID()[bFI];
+        label patchi = mesh.boundaryMesh().patchID()[bFI];
 
         if (patches[patchi].coupled())
         {
@@ -381,7 +379,7 @@ bool Foam::polyMeshTetDecomposition::checkFaceTets
     // Calculate coupled cell centre
     pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
 
-    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
+    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
     {
         neiCc[facei - mesh.nInternalFaces()] = cc[own[facei]];
     }
@@ -526,21 +524,21 @@ bool Foam::polyMeshTetDecomposition::checkFaceTets
 Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
 (
     const polyMesh& mesh,
-    label fI,
-    label cI
+    label facei,
+    label celli
 )
 {
     const faceList& pFaces = mesh.faces();
 
-    const face& f = pFaces[fI];
+    const face& f = pFaces[facei];
 
     label nTets = f.size() - 2;
 
     List<tetIndices> faceTets(nTets);
 
-    for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI ++)
+    for (label tetPti = 1; tetPti < f.size() - 1; ++tetPti)
     {
-        faceTets[tetPtI - 1] = tetIndices(cI, fI, tetPtI);
+        faceTets[tetPti - 1] = tetIndices(celli, facei, tetPti);
     }
 
     return faceTets;
@@ -550,28 +548,26 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
 Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::cellTetIndices
 (
     const polyMesh& mesh,
-    label cI
+    label celli
 )
 {
     const faceList& pFaces = mesh.faces();
     const cellList& pCells = mesh.cells();
 
-    const cell& thisCell = pCells[cI];
+    const cell& thisCell = pCells[celli];
 
     label nTets = 0;
 
-    forAll(thisCell, cFI)
+    for (const label facei : thisCell)
     {
-        nTets += pFaces[thisCell[cFI]].size() - 2;
+        nTets += pFaces[facei].size() - 2;
     }
 
     DynamicList<tetIndices> cellTets(nTets);
 
-    forAll(thisCell, cFI)
+    for (const label facei : thisCell)
     {
-        label fI = thisCell[cFI];
-
-        cellTets.append(faceTetIndices(mesh, fI, cI));
+        cellTets.append(faceTetIndices(mesh, facei, celli));
     }
 
     return cellTets;
@@ -581,27 +577,26 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::cellTetIndices
 Foam::tetIndices Foam::polyMeshTetDecomposition::findTet
 (
     const polyMesh& mesh,
-    label cI,
+    label celli,
     const point& pt
 )
 {
     const faceList& pFaces = mesh.faces();
     const cellList& pCells = mesh.cells();
 
-    const cell& thisCell = pCells[cI];
+    const cell& thisCell = pCells[celli];
 
     tetIndices tetContainingPt;
 
 
-    forAll(thisCell, cFI)
+    for (const label facei : thisCell)
     {
-        label fI = thisCell[cFI];
-        const face& f = pFaces[fI];
+        const face& f = pFaces[facei];
 
-        for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
+        for (label tetPti = 1; tetPti < f.size() - 1; ++tetPti)
         {
             // Get tetIndices of face triangle
-            tetIndices faceTetIs(cI, fI, tetPtI);
+            tetIndices faceTetIs(celli, facei, tetPti);
 
             // Check if inside
             if (faceTetIs.tet(mesh).inside(pt))
diff --git a/src/functionObjects/field/nearWallFields/findCellParticle.C b/src/functionObjects/field/nearWallFields/findCellParticle.C
index 7c0bd7f7c9..ade118370f 100644
--- a/src/functionObjects/field/nearWallFields/findCellParticle.C
+++ b/src/functionObjects/field/nearWallFields/findCellParticle.C
@@ -65,10 +65,11 @@ Foam::findCellParticle::findCellParticle
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    particle(mesh, is, readFields)
+    particle(mesh, is, readFields, newFormat)
 {
     if (readFields)
     {
diff --git a/src/functionObjects/field/nearWallFields/findCellParticle.H b/src/functionObjects/field/nearWallFields/findCellParticle.H
index 077335b2fe..6b7343e500 100644
--- a/src/functionObjects/field/nearWallFields/findCellParticle.H
+++ b/src/functionObjects/field/nearWallFields/findCellParticle.H
@@ -146,7 +146,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct and return a clone
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
index 02d8722560..6996e1500c 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.C
@@ -480,10 +480,11 @@ Foam::wallBoundedParticle::wallBoundedParticle
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    particle(mesh, is, readFields)
+    particle(mesh, is, readFields, newFormat)
 {
     if (readFields)
     {
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
index 654ee8ae57..c648fc5c3c 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticle.H
@@ -249,7 +249,8 @@ public:
         (
             const polyMesh& c,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct copy
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
index ef914b5b32..250e4a21b1 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.C
@@ -161,10 +161,11 @@ Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    wallBoundedParticle(mesh, is, readFields)
+    wallBoundedParticle(mesh, is, readFields, newFormat)
 {
     if (readFields)
     {
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
index 62d6a83940..09b0b21229 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedStreamLineParticle.H
@@ -181,7 +181,8 @@ public:
         (
             const polyMesh& c,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct copy
diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
index 3e1cad4338..482da35719 100644
--- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
+++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
@@ -186,7 +186,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct and return a clone
diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelIO.C b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelIO.C
index dcc6e5ecc6..944b40f798 100644
--- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelIO.C
+++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelIO.C
@@ -44,10 +44,11 @@ Foam::DSMCParcel<ParcelType>::DSMCParcel
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    ParcelType(mesh, is, readFields),
+    ParcelType(mesh, is, readFields, newFormat),
     U_(Zero),
     Ei_(0.0),
     typeId_(-1)
diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 05d1b58aed..4eba7f4cf7 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -40,12 +40,12 @@ void Foam::Cloud<ParticleType>::checkPatches() const
 {
     const polyBoundaryMesh& pbm = polyMesh_.boundaryMesh();
     bool ok = true;
-    forAll(pbm, patchi)
+    for (const polyPatch& pp : pbm)
     {
-        if (isA<cyclicAMIPolyPatch>(pbm[patchi]))
+        if (isA<cyclicAMIPolyPatch>(pp))
         {
             const cyclicAMIPolyPatch& cami =
-                refCast<const cyclicAMIPolyPatch>(pbm[patchi]);
+                refCast<const cyclicAMIPolyPatch>(pp);
 
             if (cami.owner())
             {
@@ -76,6 +76,7 @@ Foam::Cloud<ParticleType>::Cloud
 :
     cloud(pMesh, cloudName),
     IDLList<ParticleType>(),
+    geometryType_(IOPosition<Cloud<ParticleType>>::geometryType::COORDINATES),
     polyMesh_(pMesh),
     labels_(),
     globalPositionsPtr_()
@@ -113,7 +114,7 @@ void Foam::Cloud<ParticleType>::deleteParticle(ParticleType& p)
 template<class ParticleType>
 void Foam::Cloud<ParticleType>::deleteLostParticles()
 {
-    forAllIter(typename Cloud<ParticleType>, *this, pIter)
+    forAllIters(*this, pIter)
     {
         ParticleType& p = pIter();
 
@@ -170,7 +171,7 @@ void Foam::Cloud<ParticleType>::move
     }
 
     // Initialise the stepFraction moved for the particles
-    forAllIter(typename Cloud<ParticleType>, *this, pIter)
+    forAllIters(*this, pIter)
     {
         pIter().stepFraction() = 0;
     }
@@ -205,7 +206,7 @@ void Foam::Cloud<ParticleType>::move
         }
 
         // Loop over all particles
-        forAllIter(typename Cloud<ParticleType>, *this, pIter)
+        forAllIters(*this, pIter)
         {
             ParticleType& p = pIter();
 
@@ -293,9 +294,9 @@ void Foam::Cloud<ParticleType>::move
 
         bool transfered = false;
 
-        forAll(allNTrans, i)
+        for (const label n : allNTrans)
         {
-            if (allNTrans[i])
+            if (n)
             {
                 transfered = true;
                 break;
@@ -309,10 +310,8 @@ void Foam::Cloud<ParticleType>::move
         }
 
         // Retrieve from receive buffers
-        forAll(neighbourProcs, i)
+        for (const label neighbProci : neighbourProcs)
         {
-            label neighbProci = neighbourProcs[i];
-
             label nRec = allNTrans[neighbProci];
 
             if (nRec)
@@ -329,7 +328,7 @@ void Foam::Cloud<ParticleType>::move
 
                 label pI = 0;
 
-                forAllIter(typename Cloud<ParticleType>, newParticles, newpIter)
+                forAllIters(newParticles, newpIter)
                 {
                     ParticleType& newp = newpIter();
 
@@ -368,10 +367,10 @@ void Foam::Cloud<ParticleType>::autoMap(const mapPolyMesh& mapper)
     const vectorField& positions = globalPositionsPtr_();
 
     label i = 0;
-    forAllIter(typename Cloud<ParticleType>, *this, iter)
+    forAllIters(*this, iter)
     {
         iter().autoMap(positions[i], mapper);
-        ++ i;
+        ++i;
     }
 }
 
@@ -384,11 +383,12 @@ void Foam::Cloud<ParticleType>::writePositions() const
         this->db().time().path()/this->name() + "_positions.obj"
     );
 
-    forAllConstIter(typename Cloud<ParticleType>, *this, pIter)
+    forAllConstIters(*this, pIter)
     {
         const ParticleType& p = pIter();
-        pObj<< "v " << p.position().x() << " " << p.position().y() << " "
-            << p.position().z() << nl;
+        const point position(p.position());
+        pObj<< "v " << position.x() << " " << position.y() << " "
+            << position.z() << nl;
     }
 
     pObj.flush();
@@ -408,10 +408,10 @@ void Foam::Cloud<ParticleType>::storeGlobalPositions() const
     vectorField& positions = globalPositionsPtr_();
 
     label i = 0;
-    forAllConstIter(typename Cloud<ParticleType>, *this, iter)
+    forAllConstIters(*this, iter)
     {
         positions[i] = iter().position();
-        ++ i;
+        ++i;
     }
 }
 
diff --git a/src/lagrangian/basic/Cloud/Cloud.H b/src/lagrangian/basic/Cloud/Cloud.H
index ec92e831fe..a49247d0d6 100644
--- a/src/lagrangian/basic/Cloud/Cloud.H
+++ b/src/lagrangian/basic/Cloud/Cloud.H
@@ -75,6 +75,10 @@ class Cloud
 {
     // Private data
 
+        //- Geometry type
+        typename IOPosition<Cloud<ParticleType>>::geometryType geometryType_;
+
+        //- Reference to the mesh database
         const polyMesh& polyMesh_;
 
         //- Temporary storage for addressing. Used in findTris.
diff --git a/src/lagrangian/basic/Cloud/CloudIO.C b/src/lagrangian/basic/Cloud/CloudIO.C
index c530f1f940..8678c5bde6 100644
--- a/src/lagrangian/basic/Cloud/CloudIO.C
+++ b/src/lagrangian/basic/Cloud/CloudIO.C
@@ -54,6 +54,16 @@ void Foam::Cloud<ParticleType>::readCloudUniformProperties()
     {
         const IOdictionary uniformPropsDict(dictObj);
 
+        // Fall back to positions mode if the entry is not present for
+        // backwards compatibility
+        geometryType_ =
+            IOPosition<Cloud<ParticleType>>::geometryTypeNames_.lookupOrDefault
+            (
+                "geometry",
+                uniformPropsDict,
+                IOPosition<Cloud<ParticleType>>::geometryType::POSITIONS
+            );
+
         const word procName("processor" + Foam::name(Pstream::myProcNo()));
         if (uniformPropsDict.found(procName))
         {
@@ -91,6 +101,12 @@ void Foam::Cloud<ParticleType>::writeCloudUniformProperties() const
     Pstream::listCombineGather(np, maxEqOp<label>());
     Pstream::listCombineScatter(np);
 
+    uniformPropsDict.add
+    (
+        "geometry",
+        IOPosition<Cloud<ParticleType>>::geometryTypeNames_[geometryType_]
+    );
+
     forAll(np, i)
     {
         word procName("processor" + Foam::name(i));
@@ -113,7 +129,7 @@ void Foam::Cloud<ParticleType>::initCloud(const bool checkClass)
 {
     readCloudUniformProperties();
 
-    IOPosition<Cloud<ParticleType>> ioP(*this);
+    IOPosition<Cloud<ParticleType>> ioP(*this, geometryType_);
 
     bool valid = ioP.headerOk();
     Istream& is = ioP.readStream(checkClass ? typeName : "", valid);
@@ -130,6 +146,9 @@ void Foam::Cloud<ParticleType>::initCloud(const bool checkClass)
             << "Assuming the initial cloud contains 0 particles." << endl;
     }
 
+    // Always operate in co-ordinates mode after reading
+    geometryType_ = IOPosition<Cloud<ParticleType>>::geometryType::COORDINATES;
+
     // Ask for the tetBasePtIs to trigger all processors to build
     // them, otherwise, if some processors have no particles then
     // there is a comms mismatch.
@@ -148,6 +167,7 @@ Foam::Cloud<ParticleType>::Cloud
 )
 :
     cloud(pMesh, cloudName),
+    geometryType_(IOPosition<Cloud<ParticleType>>::geometryType::COORDINATES),
     polyMesh_(pMesh),
     labels_(),
     cellWallFacesPtr_()
diff --git a/src/lagrangian/basic/IOPosition/IOPosition.C b/src/lagrangian/basic/IOPosition/IOPosition.C
index 08d6e8e437..69eeca2474 100644
--- a/src/lagrangian/basic/IOPosition/IOPosition.C
+++ b/src/lagrangian/basic/IOPosition/IOPosition.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,22 +25,38 @@ License
 
 #include "IOPosition.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class CloudType>
+const Foam::Enum<typename Foam::IOPosition<CloudType>::geometryType>
+Foam::IOPosition<CloudType>::geometryTypeNames_
+{
+    { geometryType::POSITIONS, "positions" },
+    { geometryType::COORDINATES, "coordinates" }
+};
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
-Foam::IOPosition<CloudType>::IOPosition(const CloudType& c)
+Foam::IOPosition<CloudType>::IOPosition
+(
+    const CloudType& c,
+    const geometryType& geomType
+)
 :
     regIOobject
     (
         IOobject
         (
-            "positions",
+            geometryTypeNames_[geomType],
             c.time().timeName(),
             c,
             IOobject::MUST_READ,
             IOobject::NO_WRITE
         )
     ),
+    geometryType_(geomType),
     cloud_(c)
 {}
 
@@ -59,10 +75,26 @@ bool Foam::IOPosition<CloudType>::writeData(Ostream& os) const
 {
     os  << cloud_.size() << nl << token::BEGIN_LIST << nl;
 
-    forAllConstIter(typename CloudType, cloud_, iter)
+    switch (geometryType_)
     {
-        iter().writePosition(os);
-        os  << nl;
+        case geometryType::POSITIONS:
+        {
+            forAllConstIters(cloud_, iter)
+            {
+                iter().writePosition(os);
+                os  << nl;
+            }
+            break;
+        }
+        case geometryType::COORDINATES:
+        {
+            forAllConstIters(cloud_, iter)
+            {
+                iter().writeCoordinates(os);
+                os  << nl;
+            }
+            break;
+        }
     }
 
     os  << token::END_LIST << endl;
@@ -78,33 +110,39 @@ void Foam::IOPosition<CloudType>::readData(Istream& is, CloudType& c)
 
     token firstToken(is);
 
+    bool newFormat = geometryType_ == geometryType::COORDINATES;
+
     if (firstToken.isLabel())
     {
         label s = firstToken.labelToken();
 
         // Read beginning of contents
-        is.readBeginList
-        (
-            "IOPosition<CloudType>::readData(Istream&, CloudType&)"
-        );
+        is.readBeginList(FUNCTION_NAME);
 
         for (label i=0; i<s; i++)
         {
             // Read position only
-            c.append(new typename CloudType::particleType(mesh, is, false));
+            c.append
+            (
+                new typename CloudType::particleType
+                (
+                    mesh,
+                    is,
+                    false,
+                    newFormat
+                )
+            );
         }
 
         // Read end of contents
-        is.readEndList("IOPosition<CloudType>::readData(Istream&, CloudType&)");
+        is.readEndList(FUNCTION_NAME);
     }
     else if (firstToken.isPunctuation())
     {
         if (firstToken.pToken() != token::BEGIN_LIST)
         {
-            FatalIOErrorInFunction
-            (
-                is
-            )   << "incorrect first token, '(', found "
+            FatalIOErrorInFunction(is)
+                << "incorrect first token, '(', found "
                 << firstToken.info() << exit(FatalIOError);
         }
 
@@ -120,16 +158,17 @@ void Foam::IOPosition<CloudType>::readData(Istream& is, CloudType& c)
             is.putBack(lastToken);
 
             // Read position only
-            c.append(new typename CloudType::particleType(mesh, is, false));
+            c.append
+            (
+                new typename CloudType::particleType(mesh, is, false, newFormat)
+            );
             is  >> lastToken;
         }
     }
     else
     {
-        FatalIOErrorInFunction
-        (
-            is
-        )   << "incorrect first token, expected <int> or '(', found "
+        FatalIOErrorInFunction(is)
+            << "incorrect first token, expected <int> or '(', found "
             << firstToken.info() << exit(FatalIOError);
     }
 
diff --git a/src/lagrangian/basic/IOPosition/IOPosition.H b/src/lagrangian/basic/IOPosition/IOPosition.H
index fff078e18a..4f1ef401f4 100644
--- a/src/lagrangian/basic/IOPosition/IOPosition.H
+++ b/src/lagrangian/basic/IOPosition/IOPosition.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,6 +36,7 @@ SourceFiles
 #define IOPosition_H
 
 #include "regIOobject.H"
+#include "Enum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,9 +52,23 @@ class IOPosition
 :
     public regIOobject
 {
+public:
+
+    enum class geometryType
+    {
+        POSITIONS,
+        COORDINATES
+    };
+
+    static const Enum<geometryType> geometryTypeNames_;
+
+
+private:
 
     // Private data
 
+        geometryType geometryType_;
+
         //- Reference to the cloud
         const CloudType& cloud_;
 
@@ -72,7 +87,11 @@ public:
     // Constructors
 
         //- Construct from cloud
-        IOPosition(const CloudType&);
+        IOPosition
+        (
+            const CloudType& c,
+            const geometryType& geomType = geometryType::COORDINATES
+        );
 
 
     // Member functions
diff --git a/src/lagrangian/basic/indexedParticle/indexedParticle.H b/src/lagrangian/basic/indexedParticle/indexedParticle.H
index c619cc6694..c8e3d7634b 100644
--- a/src/lagrangian/basic/indexedParticle/indexedParticle.H
+++ b/src/lagrangian/basic/indexedParticle/indexedParticle.H
@@ -81,10 +81,11 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         )
         :
-            particle(mesh, is, readFields)
+            particle(mesh, is, readFields, newFormat)
         {}
 
         //- Construct as a copy
diff --git a/src/lagrangian/basic/injectedParticle/injectedParticle.H b/src/lagrangian/basic/injectedParticle/injectedParticle.H
index 24d29a7b5f..1d8bfbf3fb 100644
--- a/src/lagrangian/basic/injectedParticle/injectedParticle.H
+++ b/src/lagrangian/basic/injectedParticle/injectedParticle.H
@@ -133,7 +133,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct as a copy
diff --git a/src/lagrangian/basic/injectedParticle/injectedParticleIO.C b/src/lagrangian/basic/injectedParticle/injectedParticleIO.C
index 39313222c9..711f393cd0 100644
--- a/src/lagrangian/basic/injectedParticle/injectedParticleIO.C
+++ b/src/lagrangian/basic/injectedParticle/injectedParticleIO.C
@@ -48,10 +48,11 @@ Foam::injectedParticle::injectedParticle
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    particle(mesh, is, readFields),
+    particle(mesh, is, readFields, newFormat),
     tag_(-1),
     soi_(0.0),
     d_(0.0),
diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index c9aab06e06..5561ac9f15 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.C
@@ -37,6 +37,10 @@ Foam::label Foam::particle::particleCount_ = 0;
 namespace Foam
 {
     defineTypeNameAndDebug(particle, 0);
+    bool particle::writeLagrangianPositions
+    (
+        debug::infoSwitch("writeLagrangianPositions", 0)
+    );
 }
 
 
@@ -540,14 +544,14 @@ void Foam::particle::locate
         if (nWarnings < maxNWarnings)
         {
             WarningInFunction << boundaryMsg << endl;
-            ++ nWarnings;
+            ++nWarnings;
         }
         if (nWarnings == maxNWarnings)
         {
             WarningInFunction
                 << "Suppressing any further warnings about particles being "
                 << "located outside of the mesh." << endl;
-            ++ nWarnings;
+            ++nWarnings;
         }
     }
 }
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index bd69473b45..9ddea7bbd7 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -87,6 +87,9 @@ class particle
 {
     // Private member data
 
+        //- Write particle positions file (v1706 format and earlier)
+        static bool writeLagrangianPositions;
+
         //- Size in bytes of the position data
         static const std::size_t sizeofPosition;
 
@@ -368,7 +371,13 @@ public:
         );
 
         //- Construct from Istream
-        particle(const polyMesh& mesh, Istream&, bool readFields = true);
+        particle
+        (
+            const polyMesh& mesh,
+            Istream&,
+            bool readFields = true,
+            bool newFormat = true
+        );
 
         //- Construct as a copy
         particle(const particle& p);
@@ -558,7 +567,7 @@ public:
             trackingData& td
         );
 
-        //- Convenience function. Cobines trackToFace and hitFace
+        //- Convenience function. Combines trackToFace and hitFace
         template<class TrackCloudType>
         void trackToAndHitFace
         (
@@ -650,6 +659,9 @@ public:
         template<class CloudType>
         static void writeObjects(const CloudType& c, objectRegistry& obr);
 
+        //- Write the particle barycentric coordinates and cell info
+        void writeCoordinates(Ostream&) const;
+
         //- Write the particle position and cell
         void writePosition(Ostream&) const;
 
diff --git a/src/lagrangian/basic/particle/particleIO.C b/src/lagrangian/basic/particle/particleIO.C
index 4277457cab..15ab90ff9e 100644
--- a/src/lagrangian/basic/particle/particleIO.C
+++ b/src/lagrangian/basic/particle/particleIO.C
@@ -44,7 +44,13 @@ const std::size_t Foam::particle::sizeofFields
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields)
+Foam::particle::particle
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields,
+    bool newFormat
+)
 :
     mesh_(mesh),
     coordinates_(),
@@ -56,25 +62,88 @@ Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields)
     origProc_(Pstream::myProcNo()),
     origId_(-1)
 {
-    if (is.format() == IOstream::ASCII)
+    if (newFormat)
     {
-        is  >> coordinates_ >> celli_ >> tetFacei_ >> tetPti_;
-
-        if (readFields)
+        if (is.format() == IOstream::ASCII)
         {
-            is  >> facei_ >> stepFraction_ >> origProc_ >> origId_;
+            is  >> coordinates_ >> celli_ >> tetFacei_ >> tetPti_;
+            if (readFields)
+            {
+                is  >> facei_ >> stepFraction_ >> origProc_ >> origId_;
+            }
+        }
+        else
+        {
+            if (readFields)
+            {
+                is.read(reinterpret_cast<char*>(&coordinates_), sizeofFields);
+            }
+            else
+            {
+                is.read(reinterpret_cast<char*>(&coordinates_), sizeofPosition);
+            }
         }
     }
     else
     {
-        if (readFields)
+        struct oldParticle
         {
-            is.read(reinterpret_cast<char*>(&coordinates_), sizeofFields);
+            vector position;
+            label celli;
+            label facei;
+            scalar stepFraction;
+            label tetFacei;
+            label tetPti;
+            label origProc;
+            label origId;
+        } p;
+
+        if (is.format() == IOstream::ASCII)
+        {
+            is >> p.position >> p.celli;
+
+            if (readFields)
+            {
+                is  >> p.facei
+                    >> p.stepFraction
+                    >> p.tetFacei
+                    >> p.tetPti
+                    >> p.origProc
+                    >> p.origId;
+            }
         }
         else
         {
-            is.read(reinterpret_cast<char*>(&coordinates_), sizeofPosition);
+            const size_t s =
+                offsetof(oldParticle, facei) - offsetof(oldParticle, position);
+
+            is.read(reinterpret_cast<char*>(&p.position), s);
+
+            if (readFields)
+            {
+                const size_t s =
+                    sizeof(oldParticle) - offsetof(oldParticle, facei);
+
+                is.read(reinterpret_cast<char*>(&p.facei), s);
+            }
         }
+
+        if (readFields)
+        {
+            // Note: other position-based properties are set using locate(...)
+            stepFraction_ = p.stepFraction;
+            origProc_ = p.origProc;
+            origId_ = p.origId;
+        }
+
+        locate
+        (
+            p.position,
+            nullptr,
+            p.celli,
+            false,
+            "Particle initialised with a location outside of the mesh."
+        );
     }
 
     // Check state of Istream
@@ -82,7 +151,7 @@ Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields)
 }
 
 
-void Foam::particle::writePosition(Ostream& os) const
+void Foam::particle::writeCoordinates(Ostream& os) const
 {
     if (os.format() == IOstream::ASCII)
     {
@@ -101,6 +170,41 @@ void Foam::particle::writePosition(Ostream& os) const
 }
 
 
+void Foam::particle::writePosition(Ostream& os) const
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << position() << token::SPACE << celli_;
+    }
+    else
+    {
+        struct oldParticle
+        {
+            vector position;
+            label celli;
+            label facei;
+            scalar stepFraction;
+            label tetFacei;
+            label tetPti;
+            label origProc;
+            label origId;
+        } p;
+
+        const size_t s =
+            offsetof(oldParticle, facei) - offsetof(oldParticle, position);
+
+        p.position = position();
+        p.celli = celli_;
+
+        vector pos(position());
+        os.write(reinterpret_cast<const char*>(&p.position), s);
+    }
+
+    // Check state of Ostream
+    os.check(FUNCTION_NAME);
+}
+
+
 Foam::Ostream& Foam::operator<<(Ostream& os, const particle& p)
 {
     if (os.format() == IOstream::ASCII)
diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C
index f78706aa62..0748f4df3b 100644
--- a/src/lagrangian/basic/particle/particleTemplates.C
+++ b/src/lagrangian/basic/particle/particleTemplates.C
@@ -56,7 +56,7 @@ void Foam::particle::readFields(TrackCloudType& c)
     c.checkFieldIOobject(c, origId);
 
     label i = 0;
-    forAllIter(typename TrackCloudType, c, iter)
+    forAllIters(c, iter)
     {
         particle& p = iter();
 
@@ -72,8 +72,21 @@ void Foam::particle::writeFields(const TrackCloudType& c)
 {
     label np = c.size();
 
-    IOPosition<TrackCloudType> ioP(c);
-    ioP.write(np > 0);
+    {
+        IOPosition<TrackCloudType> ioP(c);
+        ioP.write(np > 0);
+    }
+
+    // Optionally write positions file in v1706 format and earlier
+    if (writeLagrangianPositions)
+    {
+        IOPosition<TrackCloudType> ioP
+        (
+            c,
+            IOPosition<TrackCloudType>::geometryType::POSITIONS
+        );
+        ioP.write(np > 0);
+    }
 
     IOField<label> origProc
     (
@@ -87,11 +100,11 @@ void Foam::particle::writeFields(const TrackCloudType& c)
     );
 
     label i = 0;
-    forAllConstIter(typename TrackCloudType, c, iter)
+    forAllConstIters(c, iter)
     {
         origProc[i] = iter().origProc_;
         origId[i] = iter().origId_;
-        i++;
+        ++i;
     }
 
     origProc.write(np > 0);
@@ -106,13 +119,15 @@ void Foam::particle::writeObjects(const CloudType& c, objectRegistry& obr)
 
     IOField<label>& origProc(cloud::createIOField<label>("origProc", np, obr));
     IOField<label>& origId(cloud::createIOField<label>("origId", np, obr));
+    IOField<point>& position(cloud::createIOField<point>("position", np, obr));
 
     label i = 0;
-    forAllConstIter(typename CloudType, c, iter)
+    forAllConstIters(c, iter)
     {
         origProc[i] = iter().origProc_;
         origId[i] = iter().origId_;
-        i++;
+        position[i] = iter().position();
+        ++i;
     }
 }
 
@@ -140,7 +155,7 @@ void Foam::particle::hitFace
     }
     else if (onBoundaryFace())
     {
-        if(!p.hitPatch(cloud, ttd))
+        if (!p.hitPatch(cloud, ttd))
         {
             const polyPatch& patch = mesh_.boundaryMesh()[p.patch()];
 
diff --git a/src/lagrangian/basic/passiveParticle/passiveParticle.H b/src/lagrangian/basic/passiveParticle/passiveParticle.H
index 35ef33c0ae..0207cdb14a 100644
--- a/src/lagrangian/basic/passiveParticle/passiveParticle.H
+++ b/src/lagrangian/basic/passiveParticle/passiveParticle.H
@@ -88,10 +88,11 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         )
         :
-            particle(mesh, is, readFields)
+            particle(mesh, is, readFields, newFormat)
         {}
 
         //- Construct as copy
diff --git a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H
index e6c70e58c2..168d62fa0e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.H
@@ -220,7 +220,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct as a copy
diff --git a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelIO.C
index a5004adac7..dcec4a3bcd 100644
--- a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelIO.C
@@ -48,10 +48,11 @@ Foam::CollidingParcel<ParcelType>::CollidingParcel
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    ParcelType(mesh, is, readFields),
+    ParcelType(mesh, is, readFields, newFormat),
     f_(Zero),
     angularMomentum_(Zero),
     torque_(Zero),
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index e32748ae92..b215cd658b 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -399,7 +399,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct as a copy
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
index 37fdfa8c00..046015f3a0 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C
@@ -53,10 +53,11 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    ParcelType(mesh, is, readFields),
+    ParcelType(mesh, is, readFields, newFormat),
     active_(false),
     typeId_(0),
     nParticle_(0.0),
diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H
index 29ea117cec..8f5bbe459f 100644
--- a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcel.H
@@ -232,7 +232,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct as a copy
diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelIO.C
index 238cf3ebf3..97a6f446c8 100644
--- a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelIO.C
@@ -51,10 +51,11 @@ Foam::MPPICParcel<ParcelType>::MPPICParcel
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    ParcelType(mesh, is, readFields),
+    ParcelType(mesh, is, readFields, newFormat),
     UCorrect_(Zero)
 {
     if (readFields)
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index 9fef053dd9..e6a35d645e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -334,7 +334,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct as a copy
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
index 6f3633c031..20e1e4adfe 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
@@ -50,10 +50,11 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    ParcelType(mesh, is, readFields),
+    ParcelType(mesh, is, readFields, newFormat),
     YGas_(0),
     YLiquid_(0),
     YSolid_(0),
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index 812b2c7ef5..d12b029a89 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -281,7 +281,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct as a copy
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
index 9541a85df8..6859f25545 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C
@@ -50,10 +50,11 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    ParcelType(mesh, is, readFields),
+    ParcelType(mesh, is, readFields, newFormat),
     mass0_(0.0),
     Y_(0)
 {
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index c45c3377e2..c76cd6893f 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -342,7 +342,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct as a copy
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
index 0b658e033a..d8261c8127 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C
@@ -51,10 +51,11 @@ Foam::ThermoParcel<ParcelType>::ThermoParcel
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    ParcelType(mesh, is, readFields),
+    ParcelType(mesh, is, readFields, newFormat),
     T_(0.0),
     Cp_(0.0)
 {
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
index f7af9ec6f3..a33f20ed2d 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
@@ -276,7 +276,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct and return a clone
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
index e1c44b0587..5dbc89f7ab 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C
@@ -41,10 +41,11 @@ Foam::molecule::molecule
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    particle(mesh, is, readFields),
+    particle(mesh, is, readFields, newFormat),
     Q_(Zero),
     v_(Zero),
     a_(Zero),
diff --git a/src/lagrangian/solidParticle/solidParticle.C b/src/lagrangian/solidParticle/solidParticle.C
index 94cbc02b79..568c1b4c7d 100644
--- a/src/lagrangian/solidParticle/solidParticle.C
+++ b/src/lagrangian/solidParticle/solidParticle.C
@@ -44,8 +44,6 @@ bool Foam::solidParticle::move
     td.switchProcessor = false;
     td.keepParticle = true;
 
-    const polyBoundaryMesh& pbMesh = mesh().boundaryMesh();
-
     while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
     {
         if (debug)
diff --git a/src/lagrangian/solidParticle/solidParticle.H b/src/lagrangian/solidParticle/solidParticle.H
index 636ff14760..232f7db690 100644
--- a/src/lagrangian/solidParticle/solidParticle.H
+++ b/src/lagrangian/solidParticle/solidParticle.H
@@ -146,7 +146,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct and return a clone
diff --git a/src/lagrangian/solidParticle/solidParticleIO.C b/src/lagrangian/solidParticle/solidParticleIO.C
index 537bf9d7dc..285681f666 100644
--- a/src/lagrangian/solidParticle/solidParticleIO.C
+++ b/src/lagrangian/solidParticle/solidParticleIO.C
@@ -40,10 +40,11 @@ Foam::solidParticle::solidParticle
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    particle(mesh, is, readFields)
+    particle(mesh, is, readFields, newFormat)
 {
     if (readFields)
     {
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
index 0fb14b2cf4..db358e4eb0 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H
@@ -277,7 +277,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct as a copy
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelIO.C b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelIO.C
index 3e1cc6c372..7fca953769 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelIO.C
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelIO.C
@@ -51,10 +51,11 @@ Foam::SprayParcel<ParcelType>::SprayParcel
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    ParcelType(mesh, is, readFields),
+    ParcelType(mesh, is, readFields, newFormat),
     d0_(0.0),
     position0_(Zero),
     sigma_(0.0),
diff --git a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C
index 7bc85fc137..4bd52700e1 100644
--- a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C
+++ b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.C
@@ -78,10 +78,11 @@ Foam::trackedParticle::trackedParticle
 (
     const polyMesh& mesh,
     Istream& is,
-    bool readFields
+    bool readFields,
+    bool newFormat
 )
 :
-    particle(mesh, is, readFields)
+    particle(mesh, is, readFields, newFormat)
 {
     if (readFields)
     {
diff --git a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H
index 1718db8c53..d5266f4a66 100644
--- a/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H
+++ b/src/mesh/snappyHexMesh/trackedParticle/trackedParticle.H
@@ -153,7 +153,8 @@ public:
         (
             const polyMesh& mesh,
             Istream& is,
-            bool readFields = true
+            bool readFields = true,
+            bool newFormat = true
         );
 
         //- Construct and return a clone