diff --git a/applications/solvers/incompressible/pimpleDyMFoam/pimpleDyMFoam.C b/applications/solvers/incompressible/pimpleDyMFoam/pimpleDyMFoam.C
index 9dc779fa23..66692cc4a8 100644
--- a/applications/solvers/incompressible/pimpleDyMFoam/pimpleDyMFoam.C
+++ b/applications/solvers/incompressible/pimpleDyMFoam/pimpleDyMFoam.C
@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    turbDyMFoam
+    pimpleDyMFoam.C
 
 Description
     Transient solver for incompressible, flow of Newtonian fluids
diff --git a/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C b/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C
index 35eae96a58..ac32c1ccc9 100644
--- a/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C
+++ b/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C
@@ -48,6 +48,7 @@ Description
 #include "OFstream.H"
 #include "IFstream.H"
 #include "IOobjectList.H"
+#include "SortableList.H"
 
 using namespace Foam;
 
@@ -107,6 +108,7 @@ int main(int argc, char *argv[])
     {
         // Not in memory. Load it.
         pointSet set(*iter());
+        SortableList<label> pointLabels(set.toc());
 
         label zoneID = mesh.pointZones().findZoneID(set.name());
         if (zoneID == -1)
@@ -120,7 +122,7 @@ int main(int argc, char *argv[])
                 new pointZone
                 (
                     set.name(),             //name
-                    set.toc(),              //addressing
+                    pointLabels,            //addressing
                     sz,                     //index
                     mesh.pointZones()       //pointZoneMesh
                 )
@@ -131,7 +133,7 @@ int main(int argc, char *argv[])
         {
             Info<< "Overwriting contents of existing pointZone " << zoneID
                 << " with that of set " << set.name() << "." << endl;
-            mesh.pointZones()[zoneID] = set.toc();
+            mesh.pointZones()[zoneID] = pointLabels;
             mesh.pointZones().writeOpt() = IOobject::AUTO_WRITE;
         }
     }
@@ -150,6 +152,7 @@ int main(int argc, char *argv[])
     {
         // Not in memory. Load it.
         cellSet set(*iter());
+        SortableList<label> cellLabels(set.toc());
 
         label zoneID = mesh.cellZones().findZoneID(set.name());
         if (zoneID == -1)
@@ -163,7 +166,7 @@ int main(int argc, char *argv[])
                 new cellZone
                 (
                     set.name(),             //name
-                    set.toc(),              //addressing
+                    cellLabels,             //addressing
                     sz,                     //index
                     mesh.cellZones()        //pointZoneMesh
                 )
@@ -174,7 +177,7 @@ int main(int argc, char *argv[])
         {
             Info<< "Overwriting contents of existing cellZone " << zoneID
                 << " with that of set " << set.name() << "." << endl;
-            mesh.cellZones()[zoneID] = set.toc();
+            mesh.cellZones()[zoneID] = cellLabels;
             mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
         }
     }
@@ -193,6 +196,7 @@ int main(int argc, char *argv[])
     {
         // Not in memory. Load it.
         faceSet set(*iter());
+        SortableList<label> faceLabels(set.toc());
 
         DynamicList<label> addressing(set.size());
         DynamicList<bool> flipMap(set.size());
@@ -214,9 +218,9 @@ int main(int argc, char *argv[])
             // Load corresponding cells
             cellSet cells(mesh, setName);
 
-            forAllConstIter(faceSet, set, iter)
+            forAll(faceLabels, i)
             {
-                label faceI = iter.key();
+                label faceI = faceLabels[i];
 
                 bool flip = false;
 
@@ -273,9 +277,10 @@ int main(int argc, char *argv[])
         else
         {
             // No flip map.
-            forAllConstIter(faceSet, set, iter)
+            forAll(faceLabels, i)
             {
-                addressing.append(iter.key());
+                label faceI = faceLabels[i];
+                addressing.append(faceI);
                 flipMap.append(false);
             }
         }
diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
index 8458e7f570..9311c9b4e9 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
@@ -51,6 +51,7 @@ int main(int argc, char *argv[])
     argList::noParallel();
 #   include "addRegionOption.H"
     argList::validOptions.insert("fields", "\"(list of fields)\"");
+    argList::validOptions.insert("noLagrangian", "");
 
 #   include "setRootCase.H"
 #   include "createTime.H"
@@ -61,6 +62,8 @@ int main(int argc, char *argv[])
         IStringStream(args.options()["fields"])() >> selectedFields;
     }
 
+    bool noLagrangian = args.options().found("noLagrangian");
+
     // determine the processor count directly
     label nProcs = 0;
     while (isDir(args.path()/(word("processor") + name(nProcs))))
@@ -269,118 +272,121 @@ int main(int argc, char *argv[])
         // the first processor that has them. They are in pass2 only used
         // for name and type (scalar, vector etc).
 
-        HashTable<IOobjectList> cloudObjects;
-
-        forAll (databases, procI)
+        if (!noLagrangian)
         {
-            fileNameList cloudDirs
-            (
-                readDir
-                (
-                    databases[procI].timePath()/regionPrefix/cloud::prefix,
-                    fileName::DIRECTORY
-                )
-            );
+            HashTable<IOobjectList> cloudObjects;
 
-            forAll (cloudDirs, i)
+            forAll (databases, procI)
             {
-                // Check if we already have cloud objects for this cloudname.
-                HashTable<IOobjectList>::const_iterator iter =
-                    cloudObjects.find(cloudDirs[i]);
-
-                if (iter == cloudObjects.end())
-                {
-                    // Do local scan for valid cloud objects.
-                    IOobjectList sprayObjs
+                fileNameList cloudDirs
+                (
+                    readDir
                     (
-                        procMeshes.meshes()[procI],
-                        databases[procI].timeName(),
-                        cloud::prefix/cloudDirs[i]
-                    );
+                        databases[procI].timePath()/regionPrefix/cloud::prefix,
+                        fileName::DIRECTORY
+                    )
+                );
 
-                    IOobject* positionsPtr = sprayObjs.lookup("positions");
+                forAll (cloudDirs, i)
+                {
+                    // Check if we already have cloud objects for this cloudname
+                    HashTable<IOobjectList>::const_iterator iter =
+                        cloudObjects.find(cloudDirs[i]);
 
-                    if (positionsPtr)
+                    if (iter == cloudObjects.end())
                     {
-                        cloudObjects.insert(cloudDirs[i], sprayObjs);
+                        // Do local scan for valid cloud objects
+                        IOobjectList sprayObjs
+                        (
+                            procMeshes.meshes()[procI],
+                            databases[procI].timeName(),
+                            cloud::prefix/cloudDirs[i]
+                        );
+
+                        IOobject* positionsPtr = sprayObjs.lookup("positions");
+
+                        if (positionsPtr)
+                        {
+                            cloudObjects.insert(cloudDirs[i], sprayObjs);
+                        }
                     }
                 }
             }
-        }
 
 
-        if (cloudObjects.size())
-        {
-            // Pass2: reconstruct the cloud
-            forAllConstIter(HashTable<IOobjectList>, cloudObjects, iter)
+            if (cloudObjects.size())
             {
-                const word cloudName = string::validate<word>(iter.key());
+                // Pass2: reconstruct the cloud
+                forAllConstIter(HashTable<IOobjectList>, cloudObjects, iter)
+                {
+                    const word cloudName = string::validate<word>(iter.key());
 
-                // Objects (on arbitrary processor)
-                const IOobjectList& sprayObjs = iter();
+                    // Objects (on arbitrary processor)
+                    const IOobjectList& sprayObjs = iter();
 
-                Info<< "Reconstructing lagrangian fields for cloud "
-                    << cloudName << nl << endl;
+                    Info<< "Reconstructing lagrangian fields for cloud "
+                        << cloudName << nl << endl;
 
-                reconstructLagrangianPositions
-                (
-                    mesh,
-                    cloudName,
-                    procMeshes.meshes(),
-                    procMeshes.faceProcAddressing(),
-                    procMeshes.cellProcAddressing()
-                );
-                reconstructLagrangianFields<label>
-                (
-                    cloudName,
-                    mesh,
-                    procMeshes.meshes(),
-                    sprayObjs
-                );
-                reconstructLagrangianFields<scalar>
-                (
-                    cloudName,
-                    mesh,
-                    procMeshes.meshes(),
-                    sprayObjs
-                );
-                reconstructLagrangianFields<vector>
-                (
-                    cloudName,
-                    mesh,
-                    procMeshes.meshes(),
-                    sprayObjs
-                );
-                reconstructLagrangianFields<sphericalTensor>
-                (
-                    cloudName,
-                    mesh,
-                    procMeshes.meshes(),
-                    sprayObjs
-                );
-                reconstructLagrangianFields<symmTensor>
-                (
-                    cloudName,
-                    mesh,
-                    procMeshes.meshes(),
-                    sprayObjs
-                );
-                reconstructLagrangianFields<tensor>
-                (
-                    cloudName,
-                    mesh,
-                    procMeshes.meshes(),
-                    sprayObjs
-                );
+                    reconstructLagrangianPositions
+                    (
+                        mesh,
+                        cloudName,
+                        procMeshes.meshes(),
+                        procMeshes.faceProcAddressing(),
+                        procMeshes.cellProcAddressing()
+                    );
+                    reconstructLagrangianFields<label>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
+                    reconstructLagrangianFields<scalar>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
+                    reconstructLagrangianFields<vector>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
+                    reconstructLagrangianFields<sphericalTensor>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
+                    reconstructLagrangianFields<symmTensor>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
+                    reconstructLagrangianFields<tensor>
+                    (
+                        cloudName,
+                        mesh,
+                        procMeshes.meshes(),
+                        sprayObjs
+                    );
+                }
+            }
+            else
+            {
+                Info << "No lagrangian fields" << nl << endl;
             }
-        }
-        else
-        {
-            Info << "No lagrangian fields" << nl << endl;
         }
 
         // If there are any "uniform" directories copy them from
-        // the master processor.
+        // the master processor
 
         fileName uniformDir0 = databases[0].timePath()/"uniform";
         if (isDir(uniformDir0))
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
index 0d99f3e95b..47e919c9c3 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
@@ -50,6 +50,7 @@ vtkPV3FoamReader::vtkPV3FoamReader()
 
     output0_  = NULL;
 
+#ifdef VTKPV3FOAM_DUALPORT
     // Add second output for the Lagrangian
     this->SetNumberOfOutputPorts(2);
     vtkMultiBlockDataSet *lagrangian = vtkMultiBlockDataSet::New();
@@ -57,6 +58,7 @@ vtkPV3FoamReader::vtkPV3FoamReader()
 
     this->GetExecutive()->SetOutputData(1, lagrangian);
     lagrangian->Delete();
+#endif
 
     TimeStepRange[0] = 0;
     TimeStepRange[1] = 0;
@@ -319,15 +321,6 @@ int vtkPV3FoamReader::RequestData
         )
     );
 
-    vtkMultiBlockDataSet* lagrangianOutput = vtkMultiBlockDataSet::SafeDownCast
-    (
-        outputVector->GetInformationObject(1)->Get
-        (
-            vtkMultiBlockDataSet::DATA_OBJECT()
-        )
-    );
-
-
     if (Foam::vtkPV3Foam::debug)
     {
         cout<< "update output with "
@@ -376,7 +369,21 @@ int vtkPV3FoamReader::RequestData
 
 #else
 
-    foamData_->Update(output, lagrangianOutput);
+#ifdef VTKPV3FOAM_DUALPORT
+    foamData_->Update
+    (
+        output,
+        vtkMultiBlockDataSet::SafeDownCast
+        (
+            outputVector->GetInformationObject(1)->Get
+            (
+                vtkMultiBlockDataSet::DATA_OBJECT()
+            )
+        );
+    );
+#else
+    foamData_->Update(output, output);
+#endif
 
     if (ShowPatchNames)
     {
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
index 624358b4a0..118e99fc91 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
@@ -465,7 +465,6 @@ void Foam::vtkPV3Foam::Update
         cout<< "<beg> Foam::vtkPV3Foam::Update - output with "
             << output->GetNumberOfBlocks() << " and "
             << lagrangianOutput->GetNumberOfBlocks() << " blocks\n";
-
         output->Print(cout);
         lagrangianOutput->Print(cout);
         printMemory();
@@ -504,8 +503,10 @@ void Foam::vtkPV3Foam::Update
         reader_->UpdateProgress(0.7);
     }
 
+#ifdef VTKPV3FOAM_DUALPORT
     // restart port1 at block=0
     blockNo = 0;
+#endif
     convertMeshLagrangian(lagrangianOutput, blockNo);
 
     reader_->UpdateProgress(0.8);
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
index c6dd65af9e..4e203eacec 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.H
@@ -71,6 +71,8 @@ SourceFiles
 #include "PrimitivePatchInterpolation.H"
 #include "volPointInterpolation.H"
 
+#undef VTKPV3FOAM_DUALPORT
+
 // * * * * * * * * * * * * * Forward Declarations  * * * * * * * * * * * * * //
 
 class vtkDataArraySelection;
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
index a7cfc94e2c..4925a31284 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
@@ -326,10 +326,10 @@ void Foam::vtkPV3Foam::convertVolField
     const labelList& superCells = decompInfo.superCells();
 
     vtkFloatArray* celldata = vtkFloatArray::New();
-    celldata->SetNumberOfTuples( superCells.size() );
-    celldata->SetNumberOfComponents( nComp );
-    celldata->Allocate( nComp*superCells.size() );
-    celldata->SetName( tf.name().c_str() );
+    celldata->SetNumberOfTuples(superCells.size());
+    celldata->SetNumberOfComponents(nComp);
+    celldata->Allocate(nComp*superCells.size());
+    celldata->SetName(tf.name().c_str());
 
     if (debug)
     {
diff --git a/src/OpenFOAM/db/dictionary/dictionary.C b/src/OpenFOAM/db/dictionary/dictionary.C
index 8c351510e7..fd89f0f608 100644
--- a/src/OpenFOAM/db/dictionary/dictionary.C
+++ b/src/OpenFOAM/db/dictionary/dictionary.C
@@ -529,7 +529,7 @@ bool Foam::dictionary::add(entry* entryPtr, bool mergeEntry)
             }
             else
             {
-                IOWarningIn("dictionary::add(entry*)", (*this))
+                IOWarningIn("dictionary::add(entry*, bool)", (*this))
                     << "problem replacing entry "<< entryPtr->keyword()
                     << " in dictionary " << name() << endl;
 
@@ -558,7 +558,7 @@ bool Foam::dictionary::add(entry* entryPtr, bool mergeEntry)
     }
     else
     {
-        IOWarningIn("dictionary::add(entry* entryPtr)", (*this))
+        IOWarningIn("dictionary::add(entry*, bool)", (*this))
             << "attempt to add entry "<< entryPtr->keyword()
             << " which already exists in dictionary " << name()
             << endl;
@@ -574,11 +574,13 @@ void Foam::dictionary::add(const entry& e, bool mergeEntry)
     add(e.clone(*this).ptr(), mergeEntry);
 }
 
+
 void Foam::dictionary::add(const keyType& k, const word& w, bool overwrite)
 {
     add(new primitiveEntry(k, token(w)), overwrite);
 }
 
+
 void Foam::dictionary::add
 (
     const keyType& k,
@@ -589,16 +591,19 @@ void Foam::dictionary::add
     add(new primitiveEntry(k, token(s)), overwrite);
 }
 
+
 void Foam::dictionary::add(const keyType& k, const label l, bool overwrite)
 {
     add(new primitiveEntry(k, token(l)), overwrite);
 }
 
+
 void Foam::dictionary::add(const keyType& k, const scalar s, bool overwrite)
 {
     add(new primitiveEntry(k, token(s)), overwrite);
 }
 
+
 void Foam::dictionary::add
 (
     const keyType& k,
@@ -628,6 +633,7 @@ void Foam::dictionary::set(const entry& e)
     set(e.clone(*this).ptr());
 }
 
+
 void Foam::dictionary::set(const keyType& k, const dictionary& d)
 {
     set(new dictionaryEntry(k, *this, d));
@@ -770,7 +776,7 @@ bool Foam::dictionary::merge(const dictionary& dict)
 
     bool changed = false;
 
-    forAllConstIter(IDLList<entry>, *this, iter)
+    forAllConstIter(IDLList<entry>, dict, iter)
     {
         HashTable<entry*>::iterator fnd = hashedEntries_.find(iter().keyword());
 
diff --git a/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C b/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C
index aa6d0d750f..9ea1307ee1 100644
--- a/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C
+++ b/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C
@@ -67,11 +67,25 @@ Foam::interpolationTable<Type>::interpolationTable()
 
 
 template<class Type>
-Foam::interpolationTable<Type>::interpolationTable(const fileName& fn)
+Foam::interpolationTable<Type>::interpolationTable
+(
+    const List<Tuple2<scalar, Type> >& values,
+    const boundsHandling bounds,
+    const fileName& fName
+)
+:
+    List<Tuple2<scalar, Type> >(values),
+    boundsHandling_(bounds),
+    fileName_(fName)
+{}
+
+
+template<class Type>
+Foam::interpolationTable<Type>::interpolationTable(const fileName& fName)
 :
     List<Tuple2<scalar, Type> >(),
     boundsHandling_(interpolationTable::WARN),
-    fileName_(fn)
+    fileName_(fName)
 {
     readTable();
 }
diff --git a/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.H b/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.H
index 2637ccf808..884c3cabd7 100644
--- a/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.H
+++ b/src/OpenFOAM/interpolations/interpolationTable/interpolationTable.H
@@ -102,8 +102,16 @@ public:
         //- Construct null
         interpolationTable();
 
+        //- Construct from components
+        interpolationTable
+        (
+            const List<Tuple2<scalar, Type> >& values,
+            const boundsHandling bounds,
+            const fileName& fName
+        );
+
         //- Construct given the name of the file containing the table of data
-        interpolationTable(const fileName& fn);
+        interpolationTable(const fileName& fName);
 
         //- Construct by reading the fileName and boundsHandling from dictionary
         //  and read the table from that file.
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C
index 6ef8787d77..eee664b1fb 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C
@@ -236,6 +236,7 @@ void processorPointPatch::initPatchPatchPoints()
     );
 
     toNeighbProc
+        << ppmp.size()              // number of points for checking
         << patchPatchPoints
         << patchPatchPointNormals;
 
@@ -257,13 +258,33 @@ void Foam::processorPointPatch::calcPatchPatchPoints()
         neighbProcNo()
     );
 
+    label nbrNPoints(readLabel(fromNeighbProc));
     labelListList patchPatchPoints(fromNeighbProc);
     List<List<vector> > patchPatchPointNormals(fromNeighbProc);
 
     pointBoundaryMesh& pbm = const_cast<pointBoundaryMesh&>(boundaryMesh());
-
     const labelList& ppmp = meshPoints();
 
+    // Simple check for the very rare situation when not the same number
+    // of points on both sides. This can happen with decomposed cyclics.
+    // If on one side the cyclic shares a point with proc faces coming from
+    // internal faces it will have a different number of points from
+    // the situation where the cyclic and the 'normal' proc faces are fully
+    // separate.
+    if (nbrNPoints != ppmp.size())
+    {
+        WarningIn("processorPointPatch::calcPatchPatchPoints()")
+            << "Processor patch " << name()
+            << " has " << ppmp.size() << " points; coupled patch has "
+            << nbrNPoints << " points." << endl
+            << "   (usually due to decomposed cyclics)."
+            << " This might give problems" << endl
+            << "    when using point fields (interpolation, mesh motion)."
+            << endl;
+    }
+
+
+
     // Loop over the patches looking for other patches that share points
     forAll(patchPatchPoints, patchi)
     {
diff --git a/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.C b/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.C
index aff80a2ac3..d0b81bc8c2 100644
--- a/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.C
+++ b/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.C
@@ -32,6 +32,7 @@ defineTypeNameAndDebug(Foam::dynamicFvMesh, 0);
 
 defineRunTimeSelectionTable(Foam::dynamicFvMesh, IOobject);
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::dynamicFvMesh::dynamicFvMesh(const IOobject& io)
@@ -45,4 +46,5 @@ Foam::dynamicFvMesh::dynamicFvMesh(const IOobject& io)
 Foam::dynamicFvMesh::~dynamicFvMesh()
 {}
 
+
 // ************************************************************************* //
diff --git a/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.H b/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.H
index dc437ea48e..04d7a69e59 100644
--- a/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.H
+++ b/src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.H
@@ -26,7 +26,7 @@ Class
     Foam::dynamicFvMesh
 
 Description
-    The dynamicFvMesh is ..
+    Abstract base class for geometry and/or topology changing fvMesh.
 
 SourceFiles
     dynamicFvMesh.C
diff --git a/src/dynamicMesh/attachDetach/attachInterface.C b/src/dynamicMesh/attachDetach/attachInterface.C
index 21e4f882fd..4c4294f8f3 100644
--- a/src/dynamicMesh/attachDetach/attachInterface.C
+++ b/src/dynamicMesh/attachDetach/attachInterface.C
@@ -207,7 +207,20 @@ void Foam::attachDetach::attachInterface
                     mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
                 ];
         }
-            
+
+
+        label patchID = mesh.boundaryMesh().whichPatch(curFaceID);
+        label neiCell;
+        if (patchID == -1)
+        {
+            neiCell = nei[curFaceID];
+        }
+        else
+        {
+            neiCell = -1;
+        }
+
+
         // Modify the face
         ref.setAction
         (
@@ -216,9 +229,9 @@ void Foam::attachDetach::attachInterface
                 newFace,                // modified face
                 curFaceID,              // label of face being modified
                 own[curFaceID],         // owner
-                nei[curFaceID],         // neighbour
+                neiCell,                // neighbour
                 false,                  // face flip
-                mesh.boundaryMesh().whichPatch(curFaceID),// patch for face
+                patchID,                // patch for face
                 false,                  // remove from zone
                 modifiedFaceZone,       // zone for face
                 modifiedFaceZoneFlip    // face flip in zone
diff --git a/src/dynamicMesh/attachDetach/detachInterface.C b/src/dynamicMesh/attachDetach/detachInterface.C
index 7e0ac0a933..6b57ccb685 100644
--- a/src/dynamicMesh/attachDetach/detachInterface.C
+++ b/src/dynamicMesh/attachDetach/detachInterface.C
@@ -73,6 +73,37 @@ void Foam::attachDetach::detachInterface
     const polyMesh& mesh = topoChanger().mesh();
     const faceZoneMesh& zoneMesh = mesh.faceZones();
 
+    // Check that zone is in increasing order (needed since adding faces
+    // in same order - otherwise polyTopoChange face ordering will mess up
+    // correspondence)
+    if (debug)
+    {
+        const labelList& faceLabels = zoneMesh[faceZoneID_.index()];
+        if (faceLabels.size() > 0)
+        {
+            for (label i = 1; i < faceLabels.size(); i++)
+            {
+                if (faceLabels[i] <= faceLabels[i-1])
+                {
+                    FatalErrorIn
+                    (
+                        "attachDetach::detachInterface"
+                        "(polyTopoChange&) const"
+                    )   << "faceZone " << zoneMesh[faceZoneID_.index()].name()
+                        << " does not have mesh face labels in"
+                        << " increasing order." << endl
+                        << "Face label " << faceLabels[i]
+                        << " at position " << i
+                        << " is smaller than the previous value "
+                        << faceLabels[i-1]
+                        << exit(FatalError);
+                }
+            }
+        }
+    }
+
+
+
     const primitiveFacePatch& masterFaceLayer = zoneMesh[faceZoneID_.index()]();
     const pointField& points = mesh.points();
     const labelListList& meshEdgeFaces = mesh.edgeFaces();
@@ -109,14 +140,11 @@ void Foam::attachDetach::detachInterface
 
         if (edgeIsInternal)
         {
-// Pout<< "Internal edge found: (" << mp[zoneLocalEdges[curEdgeID].start()] << " " << mp[zoneLocalEdges[curEdgeID].end()] << ")" << endl;
+            const edge& e = zoneLocalEdges[curEdgeID];
 
             // Reset the point creation
-            addedPoints[zoneLocalEdges[curEdgeID].start()] =
-                mp[zoneLocalEdges[curEdgeID].start()];
-
-            addedPoints[zoneLocalEdges[curEdgeID].end()] =
-                mp[zoneLocalEdges[curEdgeID].end()];
+            addedPoints[e.start()] = mp[e.start()];
+            addedPoints[e.end()] = mp[e.end()];
         }
     }
 // Pout << "addedPoints before point creation: " << addedPoints << endl;
@@ -137,7 +165,10 @@ void Foam::attachDetach::detachInterface
                         true                       // supports a cell
                     )
                 );
-// Pout << "Adding point " << points[mp[pointI]] << " for original point " << mp[pointI] << endl;
+            //Pout<< "Adding point " << addedPoints[pointI]
+            //    << " coord1:" << points[mp[pointI]]
+            //    << " coord2:" << masterFaceLayer.localPoints()[pointI]
+            //    << " for original point " << mp[pointI] << endl;
         }
     }
 
@@ -185,6 +216,7 @@ void Foam::attachDetach::detachInterface
             );
 
             // Add renumbered face into the slave patch
+            //label addedFaceI =
             ref.setAction
             (
                 polyAddFace
@@ -201,7 +233,15 @@ void Foam::attachDetach::detachInterface
                     false                           // zone flip
                 )
             );
-// Pout << "Flip.  Modifying face: " << faces[curFaceID].reverseFace() << " next to cell: " << nei[curFaceID] << " and adding face: " << newFace << " next to cell: " << own[curFaceID] << endl;
+            //{
+            //    pointField newPts(ref.points());
+            //Pout<< "Flip.  Modifying face: " << ref.faces()[curFaceID]
+            //    << " fc:" <<  ref.faces()[curFaceID].centre(newPts)
+            //    << " next to cell: " << nei[curFaceID]
+            //    << " and adding face: " << newFace
+            //    << " fc:" << ref.faces()[addedFaceI].centre(newPts)
+            //    << " next to cell: " << own[curFaceID] << endl;
+            //}
         }
         else
         {
@@ -223,6 +263,7 @@ void Foam::attachDetach::detachInterface
             );
 
             // Add renumbered face into the slave patch
+            //label addedFaceI =
             ref.setAction
             (
                 polyAddFace
@@ -239,7 +280,15 @@ void Foam::attachDetach::detachInterface
                     false                           // face flip in zone
                 )
             );
-// Pout << "No flip.  Modifying face: " << faces[curFaceID] << " next to cell: " << own[curFaceID] << " and adding face: " << newFace << " next to cell: " << nei[curFaceID] << endl;
+            //{
+            //    pointField newPts(ref.points());
+            //Pout<< "No flip.  Modifying face: " << ref.faces()[curFaceID]
+            //    << " fc:" <<  ref.faces()[curFaceID].centre(newPts)
+            //    << " next to cell: " << own[curFaceID]
+            //    << " and adding face: " << newFace
+            //    << " fc:" << ref.faces()[addedFaceI].centre(newPts)
+            //    << " next to cell: " << nei[curFaceID] << endl;
+            //}
         }
     }
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
index f317573d8a..d2728f8b9d 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
@@ -1464,6 +1464,11 @@ void Foam::polyTopoChange::resetZones
 
             addressing[zoneI][nPoints[zoneI]++] = iter.key();
         }
+        // Sort the addressing
+        forAll(addressing, zoneI)
+        {
+            stableSort(addressing[zoneI]);
+        }
 
         // So now we both have old zones and the new addressing.
         // Invert the addressing to get pointZoneMap.
@@ -1551,6 +1556,28 @@ void Foam::polyTopoChange::resetZones
             addressing[zoneI][index] = faceI;
             flipMode[zoneI][index] = faceZoneFlip_[faceI];
         }
+        // Sort the addressing
+        forAll(addressing, zoneI)
+        {
+            labelList newToOld;
+            sortedOrder(addressing[zoneI], newToOld);
+            {
+                labelList newAddressing(addressing[zoneI].size());
+                forAll(newAddressing, i)
+                {
+                    newAddressing[i] = addressing[zoneI][newToOld[i]];
+                }
+                addressing[zoneI].transfer(newAddressing);
+            }
+            {
+                boolList newFlipMode(flipMode[zoneI].size());
+                forAll(newFlipMode, i)
+                {
+                    newFlipMode[i] = flipMode[zoneI][newToOld[i]];
+                }
+                flipMode[zoneI].transfer(newFlipMode);
+            }
+        }
 
         // So now we both have old zones and the new addressing.
         // Invert the addressing to get faceZoneFaceMap.
@@ -1644,6 +1671,11 @@ void Foam::polyTopoChange::resetZones
                 addressing[zoneI][nCells[zoneI]++] = cellI;
             }
         }
+        // Sort the addressing
+        forAll(addressing, zoneI)
+        {
+            stableSort(addressing[zoneI]);
+        }
 
         // So now we both have old zones and the new addressing.
         // Invert the addressing to get cellZoneMap.
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C
index f10681e48a..76267d4fa6 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C
@@ -229,14 +229,8 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs()
 void Foam::totalPressureFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
-    if (UName_ != "U")
-    {
-        os.writeKeyword("U") << UName_ << token::END_STATEMENT << nl;
-    }
-    if (phiName_ != "phi")
-    {
-        os.writeKeyword("phi") << phiName_ << token::END_STATEMENT << nl;
-    }
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     os.writeKeyword("rho") << rhoName_ << token::END_STATEMENT << nl;
     os.writeKeyword("psi") << psiName_ << token::END_STATEMENT << nl;
     os.writeKeyword("gamma") << gamma_ << token::END_STATEMENT << nl;
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C
index 24c3fd53d0..8eb7ca65fe 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C
@@ -259,7 +259,7 @@ void Foam::fieldAverage::readAveragingProperties()
             obr_.time().timeName(),
             "uniform",
             obr_,
-            IOobject::NO_READ,
+            IOobject::MUST_READ,
             IOobject::NO_WRITE,
             false
         );
diff --git a/src/postProcessing/functionObjects/utilities/Make/files b/src/postProcessing/functionObjects/utilities/Make/files
index 2d37c7c8f3..40a33c238d 100644
--- a/src/postProcessing/functionObjects/utilities/Make/files
+++ b/src/postProcessing/functionObjects/utilities/Make/files
@@ -1,5 +1,5 @@
-dynamicPressure/dynamicPressure.C
-dynamicPressure/dynamicPressureFunctionObject.C
+staticPressure/staticPressure.C
+staticPressure/staticPressureFunctionObject.C
 
 dsmcFields/dsmcFields.C
 dsmcFields/dsmcFieldsFunctionObject.C
diff --git a/src/postProcessing/functionObjects/utilities/dynamicPressure/IOdynamicPressure.H b/src/postProcessing/functionObjects/utilities/staticPressure/IOstaticPressure.H
similarity index 86%
rename from src/postProcessing/functionObjects/utilities/dynamicPressure/IOdynamicPressure.H
rename to src/postProcessing/functionObjects/utilities/staticPressure/IOstaticPressure.H
index 97a5c27353..8409c7f103 100644
--- a/src/postProcessing/functionObjects/utilities/dynamicPressure/IOdynamicPressure.H
+++ b/src/postProcessing/functionObjects/utilities/staticPressure/IOstaticPressure.H
@@ -23,24 +23,24 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Typedef
-    Foam::IOdynamicPressure
+    Foam::IOstaticPressure
 
 Description
-    Instance of the generic IOOutputFilter for dynamicPressure.
+    Instance of the generic IOOutputFilter for staticPressure.
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef IOdynamicPressure_H
-#define IOdynamicPressure_H
+#ifndef IOstaticPressure_H
+#define IOstaticPressure_H
 
-#include "dynamicPressure.H"
+#include "staticPressure.H"
 #include "IOOutputFilter.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    typedef IOOutputFilter<dynamicPressure> IOdynamicPressure;
+    typedef IOOutputFilter<staticPressure> IOstaticPressure;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressure.C b/src/postProcessing/functionObjects/utilities/staticPressure/staticPressure.C
similarity index 88%
rename from src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressure.C
rename to src/postProcessing/functionObjects/utilities/staticPressure/staticPressure.C
index a31c913ba7..9ca4f715df 100644
--- a/src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressure.C
+++ b/src/postProcessing/functionObjects/utilities/staticPressure/staticPressure.C
@@ -24,7 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "dynamicPressure.H"
+#include "staticPressure.H"
 #include "volFields.H"
 #include "dictionary.H"
 
@@ -32,12 +32,12 @@ License
 
 namespace Foam
 {
-    defineTypeNameAndDebug(dynamicPressure, 0);
+    defineTypeNameAndDebug(staticPressure, 0);
 }
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-bool Foam::dynamicPressure::isKinematicPressure()
+bool Foam::staticPressure::isKinematicPressure()
 {
     const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
 
@@ -47,7 +47,7 @@ bool Foam::dynamicPressure::isKinematicPressure()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::dynamicPressure::dynamicPressure
+Foam::staticPressure::staticPressure
 (
     const word& name,
     const objectRegistry& obr,
@@ -67,7 +67,7 @@ Foam::dynamicPressure::dynamicPressure
         active_ = false;
         WarningIn
         (
-            "dynamicPressure::dynamicPressure"
+            "staticPressure::staticPressure"
             "(const objectRegistry&, const dictionary&)"
         )   << "No fvMesh available, deactivating." << nl
             << endl;
@@ -80,7 +80,7 @@ Foam::dynamicPressure::dynamicPressure
             active_ = false;
             WarningIn
             (
-                "dynamicPressure::dynamicPressure"
+                "staticPressure::staticPressure"
                 "(const objectRegistry&, const dictionary&)"
             )   << "Pressure is not kinematic pressure, deactivating." << nl
                 << endl;
@@ -93,13 +93,13 @@ Foam::dynamicPressure::dynamicPressure
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::dynamicPressure::~dynamicPressure()
+Foam::staticPressure::~staticPressure()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::dynamicPressure::read(const dictionary& dict)
+void Foam::staticPressure::read(const dictionary& dict)
 {
     if (active_)
     {
@@ -109,19 +109,19 @@ void Foam::dynamicPressure::read(const dictionary& dict)
 }
 
 
-void Foam::dynamicPressure::execute()
+void Foam::staticPressure::execute()
 {
     // Do nothing - only valid on write
 }
 
 
-void Foam::dynamicPressure::end()
+void Foam::staticPressure::end()
 {
     // Do nothing - only valid on write
 }
 
 
-void Foam::dynamicPressure::write()
+void Foam::staticPressure::write()
 {
     if (active_)
     {
diff --git a/src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressure.H b/src/postProcessing/functionObjects/utilities/staticPressure/staticPressure.H
similarity index 82%
rename from src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressure.H
rename to src/postProcessing/functionObjects/utilities/staticPressure/staticPressure.H
index c3c6bdbb0f..0a7be2c620 100644
--- a/src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressure.H
+++ b/src/postProcessing/functionObjects/utilities/staticPressure/staticPressure.H
@@ -23,20 +23,22 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::dynamicPressure
+    Foam::staticPressure
 
 Description
-    Converts kinematic pressure to dynamic pressure, from the name of the
-    pressure field, and density.
+    Converts kinematic pressure to static pressure, from the name of the
+    pressure field, and density, i.e.
+
+        p_static = density*p_kinematic
 
 SourceFiles
-    dynamicPressure.C
-    IOdynamicPressure.H
+    staticPressure.C
+    IOstaticPressure.H
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef dynamicPressure_H
-#define dynamicPressure_H
+#ifndef staticPressure_H
+#define staticPressure_H
 
 #include "pointFieldFwd.H"
 
@@ -51,14 +53,14 @@ class dictionary;
 class mapPolyMesh;
 
 /*---------------------------------------------------------------------------*\
-                           Class dynamicPressure Declaration
+                           Class staticPressure Declaration
 \*---------------------------------------------------------------------------*/
 
-class dynamicPressure
+class staticPressure
 {
     // Private data
 
-        //- Name of this set of dynamicPressure objects
+        //- Name of this set of staticPressure objects
         word name_;
 
         const objectRegistry& obr_;
@@ -79,23 +81,23 @@ class dynamicPressure
         bool isKinematicPressure();
 
         //- Disallow default bitwise copy construct
-        dynamicPressure(const dynamicPressure&);
+        staticPressure(const staticPressure&);
 
         //- Disallow default bitwise assignment
-        void operator=(const dynamicPressure&);
+        void operator=(const staticPressure&);
 
 
 public:
 
     //- Runtime type information
-    TypeName("dynamicPressure");
+    TypeName("staticPressure");
 
 
     // Constructors
 
         //- Construct for given objectRegistry and dictionary.
         //  Allow the possibility to load fields from files
-        dynamicPressure
+        staticPressure
         (
             const word& name,
             const objectRegistry&,
@@ -106,18 +108,18 @@ public:
 
     // Destructor
 
-        virtual ~dynamicPressure();
+        virtual ~staticPressure();
 
 
     // Member Functions
 
-        //- Return name of the set of dynamicPressure
+        //- Return name of the set of staticPressure
         virtual const word& name() const
         {
             return name_;
         }
 
-        //- Read the dynamicPressure data
+        //- Read the staticPressure data
         virtual void read(const dictionary&);
 
         //- Execute, currently does nothing
@@ -126,7 +128,7 @@ public:
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
-        //- Calculate the dynamicPressure and write
+        //- Calculate the staticPressure and write
         virtual void write();
 
         //- Update for changes of mesh
diff --git a/src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressureFunctionObject.C b/src/postProcessing/functionObjects/utilities/staticPressure/staticPressureFunctionObject.C
similarity index 90%
rename from src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressureFunctionObject.C
rename to src/postProcessing/functionObjects/utilities/staticPressure/staticPressureFunctionObject.C
index f7c060f586..d84c802a4a 100644
--- a/src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressureFunctionObject.C
+++ b/src/postProcessing/functionObjects/utilities/staticPressure/staticPressureFunctionObject.C
@@ -24,18 +24,18 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "dynamicPressureFunctionObject.H"
+#include "staticPressureFunctionObject.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineNamedTemplateTypeNameAndDebug(dynamicPressureFunctionObject, 0);
+    defineNamedTemplateTypeNameAndDebug(staticPressureFunctionObject, 0);
 
     addToRunTimeSelectionTable
     (
         functionObject,
-        dynamicPressureFunctionObject,
+        staticPressureFunctionObject,
         dictionary
     );
 }
diff --git a/src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressureFunctionObject.H b/src/postProcessing/functionObjects/utilities/staticPressure/staticPressureFunctionObject.H
similarity index 81%
rename from src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressureFunctionObject.H
rename to src/postProcessing/functionObjects/utilities/staticPressure/staticPressureFunctionObject.H
index cca83745e8..40f425d4b1 100644
--- a/src/postProcessing/functionObjects/utilities/dynamicPressure/dynamicPressureFunctionObject.H
+++ b/src/postProcessing/functionObjects/utilities/staticPressure/staticPressureFunctionObject.H
@@ -23,29 +23,29 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Typedef
-    Foam::dynamicPressureFunctionObject
+    Foam::staticPressureFunctionObject
 
 Description
-    FunctionObject wrapper around dynamicPressure to allow it to be created via
+    FunctionObject wrapper around staticPressure to allow it to be created via
     the functions list within controlDict.
 
 SourceFiles
-    dynamicPressureFunctionObject.C
+    staticPressureFunctionObject.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef dynamicPressureFunctionObject_H
-#define dynamicPressureFunctionObject_H
+#ifndef staticPressureFunctionObject_H
+#define staticPressureFunctionObject_H
 
-#include "dynamicPressure.H"
+#include "staticPressure.H"
 #include "OutputFilterFunctionObject.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    typedef OutputFilterFunctionObject<dynamicPressure>
-        dynamicPressureFunctionObject;
+    typedef OutputFilterFunctionObject<staticPressure>
+        staticPressureFunctionObject;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/thermophysicalModels/liquidMixture/liquidMixture/liquidMixture.H b/src/thermophysicalModels/liquidMixture/liquidMixture/liquidMixture.H
index cac45b19ed..43ea408349 100644
--- a/src/thermophysicalModels/liquidMixture/liquidMixture/liquidMixture.H
+++ b/src/thermophysicalModels/liquidMixture/liquidMixture/liquidMixture.H
@@ -167,7 +167,7 @@ public:
         ) const;
 
 
-        //- Calculate the mean molecular weigth [kg/kmol]
+        //- Calculate the mean molecular weight [kg/kmol]
         // from mole fractions
         scalar W(const scalarField& x) const;
 
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
index 1f9a6ac63c..567eff9daf 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
@@ -45,9 +45,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF),
     TName_("undefinedT"),
-    emissivity_(0.0),
-    rayId_(0),
-    lambdaId_(0)
+    emissivity_(0.0)
 {
     refValue() = 0.0;
     refGrad() = 0.0;
@@ -66,9 +64,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf, p, iF, mapper),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_),
-    rayId_(ptf.rayId_),
-    lambdaId_(ptf.lambdaId_)
+    emissivity_(ptf.emissivity_)
 {}
 
 
@@ -82,9 +78,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF),
     TName_(dict.lookup("T")),
-    emissivity_(readScalar(dict.lookup("emissivity"))),
-    rayId_(-1),
-    lambdaId_(-1)
+    emissivity_(readScalar(dict.lookup("emissivity")))
 {
     const scalarField& Tp =
         patch().lookupPatchField<volScalarField, scalar>(TName_);
@@ -117,9 +111,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_),
-    rayId_(ptf.rayId_),
-    lambdaId_(ptf.lambdaId_)
+    emissivity_(ptf.emissivity_)
 {}
 
 
@@ -132,9 +124,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf, iF),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_),
-    rayId_(ptf.rayId_),
-    lambdaId_(ptf.lambdaId_)
+    emissivity_(ptf.emissivity_)
 {}
 
 
@@ -152,53 +142,31 @@ updateCoeffs()
         patch().lookupPatchField<volScalarField, scalar>(TName_);
 
     const radiationModel& radiation =
-            db().lookupObject<radiationModel>("radiationProperties");
+        db().lookupObject<radiationModel>("radiationProperties");
 
     const fvDOM& dom(refCast<const fvDOM>(radiation));
 
+    label rayId = -1;
+    label lambdaId = -1;
+    dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId);
+
     const label patchI = patch().index();
 
-    if (dom.nLambda() == 1)
-    {
-        if (rayId_ == -1)
-        {
-            for (label rayI=0; rayI < dom.nRay(); rayI++)
-            {
-                for (label lambdaI=0; lambdaI < dom.nLambda(); lambdaI++)
-                {
-                    const volScalarField& radiationField =
-                        dom.IRayLambda(rayI, lambdaI);
-                    if
-                    (
-                        &(radiationField.internalField())
-                     == &dimensionedInternalField()
-                    )
-                    {
-                        rayId_ = rayI;
-                        lambdaId_ = lambdaI;
-                        break;
-                    }
-                }
-            }
-        }
-    }
-    else
+    if (dom.nLambda() != 1)
     {
         FatalErrorIn
         (
             "Foam::radiation::"
-            "greyDiffusiveRadiationMixedFvPatchScalarField::"
-            "updateCoeffs"
-        )   << " a grey boundary condition is used with a non-grey"
-            << "absorption model"
-            << exit(FatalError);
+            "greyDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs"
+        )   << " a grey boundary condition is used with a non-grey "
+            << "absorption model" << nl << exit(FatalError);
     }
 
     scalarField& Iw = *this;
     vectorField n = patch().Sf()/patch().magSf();
 
     radiativeIntensityRay& ray =
-        const_cast<radiativeIntensityRay&>(dom.IRay(rayId_));
+        const_cast<radiativeIntensityRay&>(dom.IRay(rayId));
 
     ray.Qr().boundaryField()[patchI] += Iw*(-n & ray.dAve());
 
@@ -210,20 +178,22 @@ updateCoeffs()
         {
             const vector& d = dom.IRay(rayI).d();
 
-            const scalarField& Iface =
-                dom.IRay(rayI).ILambda(lambdaId_).boundaryField()[patchI];
+            const scalarField& IFace =
+                dom.IRay(rayI).ILambda(lambdaId).boundaryField()[patchI];
 
-            if ((-n[faceI] & d) < 0.0) // qin into the wall
+            if ((-n[faceI] & d) < 0.0)
             {
+                // q into the wall
                 const vector& dAve = dom.IRay(rayI).dAve();
-                Ir += Iface[faceI]*mag(n[faceI] & dAve);
+                Ir += IFace[faceI]*mag(n[faceI] & dAve);
             }
         }
 
-        const vector& d = dom.IRay(rayId_).d();
+        const vector& d = dom.IRay(rayId).d();
 
-        if ((-n[faceI] & d) > 0.) //direction out of the wall
+        if ((-n[faceI] & d) > 0.0)
         {
+            // direction out of the wall
             refGrad()[faceI] = 0.0;
             valueFraction()[faceI] = 1.0;
             refValue()[faceI] =
@@ -234,8 +204,9 @@ updateCoeffs()
                /mathematicalConstant::pi;
 
         }
-        else //direction into the wall
+        else
         {
+            // direction into the wall
             valueFraction()[faceI] = 0.0;
             refGrad()[faceI] = 0.0;
             refValue()[faceI] = 0.0; //not used
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H
index f236f8be15..005fe76865 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H
@@ -60,12 +60,6 @@ class greyDiffusiveRadiationMixedFvPatchScalarField
         //- Emissivity
         scalar emissivity_;
 
-        //- Ray index
-        label rayId_;
-
-        //- Wavelength index
-        label lambdaId_;
-
 
 public:
 
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
index eec4d22021..b3669bd6fa 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
@@ -45,9 +45,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF),
     TName_("undefinedT"),
-    emissivity_(0.0),
-    rayId_(0),
-    lambdaId_(0)
+    emissivity_(0.0)
 {
     refValue() = 0.0;
     refGrad() = 0.0;
@@ -66,9 +64,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf, p, iF, mapper),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_),
-    rayId_(ptf.rayId_),
-    lambdaId_(ptf.lambdaId_)
+    emissivity_(ptf.emissivity_)
 {}
 
 
@@ -82,9 +78,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF),
     TName_(dict.lookup("T")),
-    emissivity_(readScalar(dict.lookup("emissivity"))),
-    rayId_(0),
-    lambdaId_(0)
+    emissivity_(readScalar(dict.lookup("emissivity")))
 {
     const scalarField& Tp =
         patch().lookupPatchField<volScalarField, scalar>(TName_);
@@ -116,9 +110,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_),
-    rayId_(ptf.rayId_),
-    lambdaId_(ptf.lambdaId_)
+    emissivity_(ptf.emissivity_)
 {}
 
 
@@ -131,9 +123,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf, iF),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_),
-    rayId_(ptf.rayId_),
-    lambdaId_(ptf.lambdaId_)
+    emissivity_(ptf.emissivity_)
 {}
 
 
@@ -148,44 +138,23 @@ updateCoeffs()
     }
 
     const radiationModel& radiation =
-            db().lookupObject<radiationModel>("radiationProperties");
+        db().lookupObject<radiationModel>("radiationProperties");
 
     const fvDOM& dom(refCast<const fvDOM>(radiation));
 
+    label rayId = -1;
+    label lambdaId = -1;
+    dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId);
+
     const label patchI = patch().index();
 
-    if (dom.nLambda() > 1)
-    {
-        if (rayId_ == -1)
-        {
-            for (label rayI=0; rayI < dom.nRay(); rayI++)
-            {
-                for (label lambdaI=0; lambdaI < dom.nLambda(); lambdaI++)
-                {
-                    const volScalarField& radiationField =
-                        dom.IRayLambda(rayI, lambdaI);
-                    if
-                    (
-                        &(radiationField.internalField())
-                      ==&dimensionedInternalField()
-                    )
-                    {
-                        rayId_ = rayI;
-                        lambdaId_ = lambdaI;
-                        break;
-                    }
-                }
-            }
-        }
-    }
-    else
+    if (dom.nLambda() == 0)
     {
         FatalErrorIn
         (
             "Foam::radiation::"
-            "wideBandDiffusiveRadiationMixedFvPatchScalarField::"
-            "updateCoeffs"
-        )   << " a Non-grey boundary condition is used with a grey"
+            "wideBandDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs"
+        )   << " a non-grey boundary condition is used with a grey "
             << "absorption model" << nl << exit(FatalError);
     }
 
@@ -193,12 +162,12 @@ updateCoeffs()
     vectorField n = patch().Sf()/patch().magSf();
 
     radiativeIntensityRay& ray =
-        const_cast<radiativeIntensityRay&>(dom.IRay(rayId_));
+        const_cast<radiativeIntensityRay&>(dom.IRay(rayId));
 
     ray.Qr().boundaryField()[patchI] += Iw*(-n & ray.dAve());
 
     const scalarField Eb =
-        dom.blackBody().bLambda(lambdaId_).boundaryField()[patchI];
+        dom.blackBody().bLambda(lambdaId).boundaryField()[patchI];
 
     forAll(Iw, faceI)
     {
@@ -207,20 +176,21 @@ updateCoeffs()
         {
             const vector& d = dom.IRay(rayI).d();
 
-            const scalarField& Iface =
-                dom.IRay(rayI).ILambda(lambdaId_).boundaryField()[patchI];
+            const scalarField& IFace =
+                dom.IRay(rayI).ILambda(lambdaId).boundaryField()[patchI];
 
             if ((-n[faceI] & d) < 0.0) // qin into the wall
             {
                 const vector& dAve = dom.IRay(rayI).dAve();
-                Ir = Ir + Iface[faceI]*mag(n[faceI] & dAve);
+                Ir = Ir + IFace[faceI]*mag(n[faceI] & dAve);
             }
         }
 
-        const vector& d = dom.IRay(rayId_).d();
+        const vector& d = dom.IRay(rayId).d();
 
-        if ((-n[faceI] & d) > 0.0) //direction out of the wall
+        if ((-n[faceI] & d) > 0.0)
         {
+            // direction out of the wall
             refGrad()[faceI] = 0.0;
             valueFraction()[faceI] = 1.0;
             refValue()[faceI] =
@@ -230,8 +200,9 @@ updateCoeffs()
                 )
                /mathematicalConstant::pi;
         }
-        else  //direction into the wall
+        else
         {
+            // direction into the wall
             valueFraction()[faceI] = 0.0;
             refGrad()[faceI] = 0.0;
             refValue()[faceI] = 0.0; //not used
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H
index 3ff3990786..9d170131b2 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H
@@ -60,15 +60,6 @@ class wideBandDiffusiveRadiationMixedFvPatchScalarField
         //- Emissivity
         scalar emissivity_;
 
-        //- Ray index
-        label rayId_;
-
-        //- Wavelength index
-        label lambdaId_;
-
-        //- Radiative heat flux on walls.
-        scalarField qr_;
-
 
 public:
 
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C b/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C
index 39d02831bb..dee5647366 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C
@@ -27,19 +27,132 @@ License
 #include "blackBodyEmission.H"
 #include "dimensionedConstants.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::List<Foam::Tuple2<Foam::scalar, Foam::scalar> >
+Foam::radiation::blackBodyEmission::emissivePowerTable
+(
+    IStringStream
+    (
+        "("
+        "( 1000 0.00032)"
+        "( 1100 0.00091)"
+        "( 1200 0.00213)"
+        "( 1300 0.00432)"
+        "( 1400 0.00779)"
+        "( 1500 0.01280)"
+        "( 1600 0.01972)"
+        "( 1700 0.02853)"
+        "( 1800 0.03934)"
+        "( 1900 0.05210)"
+        "( 2000 0.06672)"
+        "( 2100 0.08305)"
+        "( 2200 0.10088)"
+        "( 2300 0.12002)"
+        "( 2400 0.14025)"
+        "( 2500 0.16135)"
+        "( 2600 0.18311)"
+        "( 2700 0.20535)"
+        "( 2800 0.22788)"
+        "( 2900 0.25055)"
+        "( 3000 0.27322)"
+        "( 3100 0.29576)"
+        "( 3200 0.31809)"
+        "( 3300 0.34009)"
+        "( 3400 0.36172)"
+        "( 3500 0.38290)"
+        "( 3600 0.40359)"
+        "( 3700 0.42375)"
+        "( 3800 0.44336)"
+        "( 3900 0.46240)"
+        "( 4000 0.48085)"
+        "( 4100 0.49872)"
+        "( 4200 0.51599)"
+        "( 4300 0.53267)"
+        "( 4400 0.54877)"
+        "( 4500 0.56429)"
+        "( 4600 0.57925)"
+        "( 4700 0.59366)"
+        "( 4800 0.60753)"
+        "( 4900 0.62088)"
+        "( 5000 0.63372)"
+        "( 5100 0.64606)"
+        "( 5200 0.65794)"
+        "( 5300 0.66935)"
+        "( 5400 0.68033)"
+        "( 5500 0.69087)"
+        "( 5600 0.70101)"
+        "( 5700 0.71076)"
+        "( 5800 0.72012)"
+        "( 5900 0.72913)"
+        "( 6000 0.73778)"
+        "( 6100 0.74610)"
+        "( 6200 0.75410)"
+        "( 6300 0.76180)"
+        "( 6400 0.76920)"
+        "( 6500 0.77631)"
+        "( 6600 0.78316)"
+        "( 6700 0.78975)"
+        "( 6800 0.79609)"
+        "( 6900 0.80219)"
+        "( 7000 0.80807)"
+        "( 7100 0.81373)"
+        "( 7200 0.81918)"
+        "( 7300 0.82443)"
+        "( 7400 0.82949)"
+        "( 7500 0.83436)"
+        "( 7600 0.83906)"
+        "( 7700 0.84359)"
+        "( 7800 0.84796)"
+        "( 7900 0.85218)"
+        "( 8000 0.85625)"
+        "( 8100 0.86017)"
+        "( 8200 0.86396)"
+        "( 8300 0.86762)"
+        "( 8400 0.87115)"
+        "( 8500 0.87456)"
+        "( 8600 0.87786)"
+        "( 8700 0.88105)"
+        "( 8800 0.88413)"
+        "( 8900 0.88711)"
+        "( 9000 0.88999)"
+        "( 9100 0.89277)"
+        "( 9200 0.89547)"
+        "( 9300 0.89807)"
+        "( 9400 0.90060)"
+        "( 9500 0.90304)"
+        "( 9600 0.90541)"
+        "( 9700 0.90770)"
+        "( 9800 0.90992)"
+        "( 9900 0.91207)"
+        "(10000 0.91415)"
+        "(12000 0.94505)"
+        "(15000 0.96893)"
+        "(20000 0.98555)"
+        "(30000 0.99529)"
+        "(40000 0.99792)"
+        "(50000 0.99890)"
+        ")"
+    )()
+);
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::radiation::blackBodyEmission::blackBodyEmission
 (
-    const fileName& name,
-    const word& instance,
     const label nLambda,
     const volScalarField& T
 )
 :
-    blackBodyEmissiveTable_(name, instance, T.mesh()),
-    C1_("C1",dimensionSet(1, 4, 3, 0, 0, 0, 0), 3.7419e-16),
-    C2_("C2",dimensionSet(0, 1, 0, 1, 0, 0, 0), 14.388e-6),
+    table_
+    (
+        emissivePowerTable,
+        interpolationTable<scalar>::CLAMP,
+        "blackBodyEmissivePower"
+    ),
+    C1_("C1", dimensionSet(1, 4, 3, 0, 0, 0, 0), 3.7419e-16),
+    C2_("C2", dimensionSet(0, 1, 0, 1, 0, 0, 0), 14.388e-6),
     bLambda_(nLambda),
     T_(T)
 {
@@ -79,7 +192,7 @@ Foam::scalar Foam::radiation::blackBodyEmission::fLambdaT
     const scalar lambdaT
 ) const
 {
-    return  blackBodyEmissiveTable_.lookUp(lambdaT*1.0e6)[1];
+    return  table_(lambdaT*1.0e6);
 }
 
 
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.H b/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.H
index d21582938f..01719e5076 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.H
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.H
@@ -28,6 +28,9 @@ Class
 Description
     Class black body emission
 
+    Table of black body emissive power taken from:
+        Modest, "Radiative Heat Transfer", pp.775-777, 1993
+
 SourceFiles
     blackBodyEmission.C
 
@@ -40,7 +43,7 @@ SourceFiles
 #include "dimensionedScalar.H"
 #include "mathematicalConstants.H"
 #include "radiationConstants.H"
-#include "interpolationLookUpTable.H"
+#include "interpolationTable.H"
 #include "Vector2D.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -56,9 +59,18 @@ namespace radiation
 
 class blackBodyEmission
 {
+public:
+
+    //- Static table of black body emissive power
+    static const List<Tuple2<scalar, scalar> > emissivePowerTable;
+
+
+private:
+
     // Private data
 
-        mutable interpolationLookUpTable<scalar> blackBodyEmissiveTable_;
+        //- Interpolation table of black body emissive power
+        mutable interpolationTable<scalar> table_;
 
         //- Constant C1
         const dimensionedScalar C1_;
@@ -85,8 +97,6 @@ public:
         //- Construct from components
         blackBodyEmission
         (
-            const fileName& name,
-            const word& instance,
             const label nLambda,
             const volScalarField& T
         );
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmissivePower b/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmissivePower
deleted file mode 100644
index c0151b436a..0000000000
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmissivePower
+++ /dev/null
@@ -1,174 +0,0 @@
-171
-(
-1000    0.00032
-1100    0.00091
-1200    0.00213
-1300    0.00432
-1400    0.00779
-1500    0.01285
-1600    0.01972
-1700    0.02853
-1800    0.03934
-1900    0.05210
-2000    0.06672
-2100    0.08305
-2200    0.10088
-2300    0.12002
-2400    0.14025
-2500    0.16135
-2600    0.18311
-2700    0.20535
-2800    0.22788
-2900    0.25055
-3000    0.27322
-3100    0.29576
-3200    0.31809
-3300    0.34009
-3400    0.36172
-3500    0.38290
-3600    0.40359
-3700    0.42375
-3800    0.44336
-3900    0.46240
-4000    0.48085
-4100    0.49872
-4200    0.51599
-4300    0.53267
-4400    0.54877
-4500    0.56429
-4600    0.57925
-4700    0.59366
-4800    0.60753
-4900    0.62088
-5000    0.63372
-5100    0.64606
-5200    0.65794
-5300    0.66935
-5400    0.68033
-5500    0.69087
-5600    0.70101
-5700    0.71076
-5800    0.72012
-5900    0.72913
-6000    0.73778
-6100    0.74610
-6200    0.75410
-6300    0.76180
-6400    0.76920
-6500    0.77631
-6600    0.78316
-6700    0.78975
-6800    0.79609
-6900    0.80219
-7000    0.80807
-7100    0.81373
-7200    0.81918
-7300    0.82443
-7400    0.82949
-7500    0.83436
-7600    0.83906
-7700    0.84359
-7800    0.84796
-7900    0.85218
-8000    0.85625
-8100    0.86017
-8200    0.86396
-8300    0.86762
-8400    0.87115
-8500    0.87456
-8600    0.87786
-8700    0.88105
-8800    0.88413
-8900    0.88711
-9000    0.88999
-9100   0.89277
-9200   0.89547
-9300   0.89807
-9400   0.90060
-9500   0.90304
-9600   0.90541
-9700   0.90770
-9800   0.90992
-9900   0.91207
-10000   0.91415
-10200   0.91813
-10400   0.92188
-10600   0.92540
-10800   0.92872
-11000   0.93184
-11200   0.93479
-11400   0.93758
-11600   0.94021
-11800   0.94270
-12000   0.94505
-12200   0.94728
-12400   0.94939
-12600   0.95139
-12800   0.95329
-13000   0.95509
-13200   0.95680
-13400   0.95843
-13600   0.95998
-13800   0.96145
-14000   0.96285
-14200   0.96418
-14400   0.96546
-14600   0.96667
-14800   0.96783
-15000   0.96893
-15200   0.96999
-15400   0.97100
-15600   0.97196
-15800   0.97288
-16000   0.97377
-16200   0.97461
-16400   0.97542
-16600   0.97620
-16800   0.97694
-17000   0.97765
-17200   0.97834
-17400   0.97899
-17600   0.97962
-17800   0.98023
-18000   0.98080
-18200   0.98137
-18400   0.98191
-18600   0.98243
-18900   0.98293
-19000   0.98340
-19200   0.98387
-19400   0.98431
-19600   0.98474
-19800   0.98515
-20000   0.98555
-21000   0.98735
-22000   0.98886
-23000   0.99014
-24000   0.99123
-25000   0.99217
-26000   0.99297
-27000   0.99367
-28000   0.99429
-29000   0.99482
-30000   0.99529
-31000   0.99571
-32000   0.99607
-33000   0.99640
-34000   0.99669
-35000   0.99695
-35000   0.99719
-36000   0.99740
-37000   0.99759
-38000   0.99776
-39000   0.99792
-40000   0.99806
-41000   0.99819
-42000   0.99831
-43000   0.99842
-44000   0.99851
-45000   0.99861
-46000   0.99869
-47000   0.99877
-48000   0.99884
-49000   0.99890
-)
\ No newline at end of file
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C
index af3b14ebb0..6d07cc3bea 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C
@@ -126,20 +126,15 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
     nRay_(0),
     nLambda_(absorptionEmission_->nBands()),
     aLambda_(nLambda_),
-    blackBody_
-    (
-        fileName("blackBodyEmissivePower"),
-        mesh_.time().constant(),
-        nLambda_,
-        T
-    ),
+    blackBody_(nLambda_, T),
     IRay_(0),
-    convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0))
+    convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
+    maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50))
 {
     if (mesh_.nSolutionD() == 3)    //3D
     {
-        IRay_.setSize(4*nPhi_*nTheta_);
         nRay_ = 4*nPhi_*nTheta_;
+        IRay_.setSize(nRay_);
         scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
         scalar deltaTheta = mathematicalConstant::pi/nTheta_;
         label i = 0;
@@ -175,8 +170,8 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         {
             scalar thetai = mathematicalConstant::pi/2.0;
             scalar deltaTheta = mathematicalConstant::pi;
-            IRay_.setSize(4*nPhi_);
             nRay_ = 4*nPhi_;
+            IRay_.setSize(nRay_);
             scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
             label i = 0;
             for (label m = 1; m <= 4*nPhi_; m++)
@@ -205,8 +200,8 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         {
             scalar thetai = mathematicalConstant::pi/2.0;
             scalar deltaTheta = mathematicalConstant::pi;
-            IRay_.setSize(2);
             nRay_ = 2;
+            IRay_.setSize(nRay_);
             scalar deltaPhi = mathematicalConstant::pi;
             label i = 0;
             for (label m = 1; m <= 2; m++)
@@ -302,9 +297,9 @@ void Foam::radiation::fvDOM::calculate()
             maxResidual = max(maxBandResidual, maxResidual);
         }
 
-        Info << "Radiation solver iter: " <<  radIter << endl;
+        Info << "Radiation solver iter: " << radIter << endl;
 
-    } while(maxResidual > convergence_);
+    } while(maxResidual > convergence_ && radIter < maxIter_);
 
     updateG();
 }
@@ -369,4 +364,20 @@ void Foam::radiation::fvDOM::updateG()
 }
 
 
+void Foam::radiation::fvDOM::setRayIdLambdaId
+(
+    const word& name,
+    label& rayId,
+    label& lambdaId
+) const
+{
+    // assuming name is in the form: CHARS_rayId_lambdaId
+    size_type i1 = name.find_first_of("_");
+    size_type i2 = name.find_last_of("_");
+
+    rayId = readLabel(IStringStream(name.substr(i1+1, i2-1))());
+    lambdaId = readLabel(IStringStream(name.substr(i2+1, name.size()-1))());
+}
+
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.H b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.H
index 3cbd921406..f91163805a 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.H
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.H
@@ -37,19 +37,18 @@ Description
     i.e. dictionary
     fvDOMCoeffs
     {
-        Nphi    1;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
-        Ntheta  2;          // polar angles in P1 (from Z to X-Y plane)
+        nPhi    1;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
+        nTheta  2;          // polar angles in PI (from Z to X-Y plane)
         convergence 1e-4;   // convergence criteria for radiation iteration
     }
 
-    nFlowIterPerRadIter   1; // Number of flow iterations per radiation
-                                iteration
+    solverFreq   1; // Number of flow iterations per radiation iteration
 
-    The total number of solid angles is  4*Nphi*Ntheta.
+    The total number of solid angles is  4*nPhi*nTheta.
 
-    In 1D the direction of the rays is X (Nphi and Ntheta are ignored)
-    In 2D the direction of the rays is on X-Y plane (only Nphi is considered)
-    In 3D (Nphi and Ntheta are considered)
+    In 1D the direction of the rays is X (nPhi and nTheta are ignored)
+    In 2D the direction of the rays is on X-Y plane (only nPhi is considered)
+    In 3D (nPhi and nTheta are considered)
 
 SourceFiles
     fvDOM.C
@@ -118,6 +117,9 @@ class fvDOM
         //- Convergence criterion
         scalar convergence_;
 
+        //- Maximum number of iterations
+        scalar maxIter_;
+
 
     // Private member functions
 
@@ -127,9 +129,6 @@ class fvDOM
         //- Disallow default bitwise assignment
         void operator=(const fvDOM&);
 
-        //- Update Absorption Coefficients
-//        void updateAbsorptionCoeffs(void);
-
         //- Update nlack body emission
         void updateBlackBodyEmission();
 
@@ -163,6 +162,15 @@ public:
             //- Update G and calculate total heat flux on boundary
             void updateG();
 
+            //- Set the rayId and lambdaId from by decomposing an intensity
+            //  field name
+            void setRayIdLambdaId
+            (
+                const word& name,
+                label& rayId,
+                label& lambdaId
+            ) const;
+
             //- Source term component (for power of T^4)
             virtual tmp<volScalarField> Rp() const;
 
@@ -194,19 +202,19 @@ public:
             //- Number of wavelengths
             inline label nLambda() const;
 
-            // Const access to total absorption coefficient
+            //- Const access to total absorption coefficient
             inline const volScalarField& a() const;
 
-            // Const access to wavelength total absorption coefficient
+            //- Const access to wavelength total absorption coefficient
             inline const volScalarField& aLambda(const label lambdaI) const;
 
-            // Const access to incident radiation field
+            //- Const access to incident radiation field
             inline const volScalarField& G() const;
 
-            //  Const access to total radiative heat flux field
+            //- Const access to total radiative heat flux field
             inline const volScalarField& Qr() const;
 
-            //  Const access to black body
+            //- Const access to black body
             inline const blackBodyEmission& blackBody() const;
 };
 
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.C b/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.C
index e9e28658c4..f000e43fa3 100755
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.C
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.C
@@ -29,23 +29,23 @@ License
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 template <class Type>
-Foam::label Foam::interpolationLookUpTable <Type>::index
+Foam::label Foam::interpolationLookUpTable<Type>::index
 (
     const List<scalar>& indices,
     const bool lastDim
 ) const
 {
-    label totalindex = 0;
+    label totalIndex = 0;
 
-    for (int i = 0; i < dim_.size() - 1; i++)
+    forAll(dim_, i)
     {
         label dim = 1;
         for (int j = i + 1; j < dim_.size(); j++)
         {
-            dim *=(dim_[j]+1);
+            dim *= dim_[j] + 1;
         }
 
-        totalindex +=
+        totalIndex +=
             dim
            *min
             (
@@ -57,7 +57,7 @@ Foam::label Foam::interpolationLookUpTable <Type>::index
     if (lastDim)
     {
         label iLastdim = dim_.size() - 1;
-        totalindex += Foam::min
+        totalIndex += Foam::min
         (
             max
             (
@@ -68,12 +68,12 @@ Foam::label Foam::interpolationLookUpTable <Type>::index
         );
     }
 
-    return totalindex;
+    return totalIndex;
 }
 
 
 template <class Type>
-Foam::label Foam::interpolationLookUpTable <Type>::index
+Foam::label Foam::interpolationLookUpTable<Type>::index
 (
     const scalar indice
 ) const
@@ -171,15 +171,15 @@ void Foam::interpolationLookUpTable<Type>::dimensionTable()
         max_[i] = readScalar(entries_[i].lookup("max"));
         min_[i] = readScalar(entries_[i].lookup("min"));
         delta_[i] = (max_[i] - min_[i])/dim_[i];
-        tableDim *= (dim_[i] + 1);
-        fieldIndices_.insert(entries_[i].lookup("name"),index);
+        tableDim *= dim_[i] + 1;
+        fieldIndices_.insert(entries_[i].lookup("name"), index);
         entryIndices_[i] = index;
         index++;
     }
 
     forAll(output_,i)
     {
-        fieldIndices_.insert(output_[i].lookup("name"),index);
+        fieldIndices_.insert(output_[i].lookup("name"), index);
         outputIndices_[i] = index;
         index++;
     }
@@ -229,8 +229,7 @@ void Foam::interpolationLookUpTable<Type>::readTable
         FatalErrorIn
         (
             "Foam::interpolationLookUpTable<Type>::readTable()"
-        )   << "table is empty" << nl
-            << exit(FatalError);
+        )   << "table is empty" << nl << exit(FatalError);
     }
 }
 
@@ -240,7 +239,7 @@ void Foam::interpolationLookUpTable<Type>::readTable
 template<class Type>
 Foam::interpolationLookUpTable<Type>::interpolationLookUpTable()
 :
-    List<scalarField >(),
+    List<scalarField>(),
     fileName_("fileNameIsUndefined")
 {}
 
@@ -251,7 +250,7 @@ Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
     const fileName& fn, const word& instance, const fvMesh& mesh
 )
 :
-    List<scalarField >(),
+    List<scalarField>(),
     fileName_(fn),
     dim_(0),
     min_(0),
@@ -273,7 +272,7 @@ Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
      const interpolationLookUpTable& interpTable
 )
 :
-    List<scalarField >(interpTable),
+    List<scalarField>(interpTable),
     fileName_(interpTable.fileName_),
     entryIndices_(interpTable.entryIndices_),
     outputIndices_(interpTable.outputIndices_),
@@ -293,7 +292,7 @@ Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
     const dictionary& dict
 )
 :
-    List<scalarField >(),
+    List<scalarField>(),
     fileName_(fileName(dict.lookup("fileName")).expand()),
     dim_(0),
     min_(0.0),
@@ -315,20 +314,19 @@ Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
 template<class Type>
 void Foam::interpolationLookUpTable<Type>::check() const
 {
-//  check order in the first dimension.
-
-    scalar prevValue = List<scalarField >::operator[](0).operator[](0);
-    label dim = 1 ;
+    // check order in the first dimension.
+    scalar prevValue = List<scalarField>::operator[](0).operator[](0);
+    label dim = 1;
     for (int j = 1; j < dim_.size(); j++)
     {
-        dim *=(dim_[j]+1);
+        dim *= dim_[j] + 1;
     }
 
     for (label i = 1; i < dim_[0]; i++)
     {
         label index = i*dim;
         const scalar currValue =
-            List<scalarField >::operator[](0).operator[](index);
+            List<scalarField>::operator[](0).operator[](index);
 
         // avoid duplicate values (divide-by-zero error)
         if (currValue <= prevValue)
@@ -336,9 +334,8 @@ void Foam::interpolationLookUpTable<Type>::check() const
             FatalErrorIn
             (
                 "Foam::interpolationLookUpTable<Type>::checkOrder() const"
-            )   << "out-of-order value: "
-                << currValue << " at index " << index << nl
-                << exit(FatalError);
+            )   << "out-of-order value: " << currValue
+                << " at index " << index << nl << exit(FatalError);
         }
         prevValue = currValue;
     }
@@ -379,8 +376,7 @@ void Foam::interpolationLookUpTable<Type>::write
         FatalErrorIn
         (
             "Foam::interpolationTable<Type>::write()"
-        )   << "table is empty" << nl
-            << exit(FatalError);
+        )   << "table is empty" << nl << exit(FatalError);
     }
     os.writeKeyword("values");
     os << *this << token::END_STATEMENT << nl;
@@ -400,31 +396,25 @@ Foam::interpolationLookUpTable<Type>::operator[](const label i)
     {
         FatalErrorIn
         (
-            "Foam::interpolationLookUpTable<Type>::operator[]"
-            "(const label) const"
-        )   << "table has (" << n << ") columns" << nl
-            << exit(FatalError);
+            "Foam::interpolationLookUpTable<Type>::operator[](const label)"
+        )   << "table has (" << n << ") columns" << nl << exit(FatalError);
     }
     else if (ii < 0)
     {
         FatalErrorIn
         (
-            "Foam::interpolationLookUpTable<Type>::operator[]"
-            "(const label) const"
-        )   << "index (" << ii << ") underflow" << nl
-            << exit(FatalError);
+            "Foam::interpolationLookUpTable<Type>::operator[](const label)"
+        )   << "index (" << ii << ") underflow" << nl << exit(FatalError);
     }
     else if (ii > n)
     {
         FatalErrorIn
         (
-            "Foam::interpolationLookUpTable<Type>::operator[]"
-            "(const label) const"
-        )   << "index (" << ii << ") overflow" << nl
-            << exit(FatalError);
+            "Foam::interpolationLookUpTable<Type>::operator[](const label)"
+        )   << "index (" << ii << ") overflow" << nl << exit(FatalError);
     }
 
-    return List<scalarField >::operator[](ii);
+    return List<scalarField>::operator[](ii);
 }
 
 
@@ -441,8 +431,7 @@ Foam::interpolationLookUpTable<Type>::operator[](const label i) const
         (
             "Foam::interpolationLookUpTable<Type>::operator[]"
             "(const label) const"
-        )   << "table has (" << n << ") columns" << nl
-            << exit(FatalError);
+        )   << "table has (" << n << ") columns" << nl << exit(FatalError);
     }
     else if (ii < 0)
     {
@@ -450,8 +439,7 @@ Foam::interpolationLookUpTable<Type>::operator[](const label i) const
         (
             "Foam::interpolationLookUpTable<Type>::operator[]"
             "(const label) const"
-        )   << "index (" << ii << ") underflow" << nl
-            << exit(FatalError);
+        )   << "index (" << ii << ") underflow" << nl << exit(FatalError);
     }
 
     else if (ii > n)
@@ -464,7 +452,7 @@ Foam::interpolationLookUpTable<Type>::operator[](const label i) const
             << exit(FatalError);
     }
 
-    return List<scalarField >::operator[](ii);
+    return List<scalarField>::operator[](ii);
 }
 
 
@@ -500,7 +488,7 @@ void Foam::interpolationLookUpTable<Type>::findHi
 
         forAll(entryIndices_,i)
         {
-            if (checkRange(retvals,entryIndices_[i]))
+            if (checkRange(retvals, entryIndices_[i]))
             {
                 label dim = 1;
 
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.H b/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.H
index d2518d4e02..6140ae4bc3 100755
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.H
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.H
@@ -54,7 +54,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                    Class interpolationLookUpTable Declaration
+                  Class interpolationLookUpTable Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type>
@@ -64,22 +64,22 @@ class interpolationLookUpTable
 {
 private:
 
-    // Private data
+    // Privsate data
 
         //- File name
         fileName fileName_;
 
         //- Table dimensions
-        List<label>  dim_;
+        List<label> dim_;
 
         //- Min on each dimension
-        List<scalar>  min_;
+        List<scalar> min_;
 
         //- Deltas on each dimension
-        List<scalar>  delta_;
+        List<scalar> delta_;
 
         //- Maximum on each dimension
-        List<scalar>  max_;
+        List<scalar> max_;
 
         //- Dictionary entries
         List<dictionary> entries_;
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C b/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
index 74ec95827c..6af99c86a4 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
@@ -37,7 +37,10 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-Foam::label Foam::radiation::radiativeIntensityRay::rayId = 0;
+Foam::label Foam::radiation::radiativeIntensityRay::rayId(0);
+
+const Foam::word
+Foam::radiation::radiativeIntensityRay::intensityPrefix("ILambda");
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -113,11 +116,14 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
         0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
     );
 
-    forAll(ILambda_, i)
+
+    autoPtr<volScalarField> IDefaultPtr;
+
+    forAll(ILambda_, lambdaI)
     {
         IOobject IHeader
         (
-            "ILambda_" + name(rayId) + "_" + name(i),
+            intensityPrefix + "_" + name(rayId) + "_" + name(lambdaI),
             mesh_.time().timeName(),
             mesh_,
             IOobject::MUST_READ,
@@ -129,29 +135,40 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
         {
             ILambda_.set
             (
-                i,
+                lambdaI,
                 new volScalarField(IHeader, mesh_)
             );
         }
         else
         {
-            volScalarField IDefault
-            (
-                IOobject
+            // Demand driven load the IDefault field
+            if (!IDefaultPtr.valid())
+            {
+                IDefaultPtr.reset
                 (
-                    "IDefault",
-                    mesh_.time().timeName(),
-                    mesh_,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh_
-            );
+                    new volScalarField
+                    (
+                        IOobject
+                        (
+                            "IDefault",
+                            mesh_.time().timeName(),
+                            mesh_,
+                            IOobject::MUST_READ,
+                            IOobject::NO_WRITE
+                        ),
+                        mesh_
+                    )
+                );
+            }
+
+            // Reset the MUST_READ flag
+            IOobject noReadHeader(IHeader);
+            noReadHeader.readOpt() = IOobject::NO_READ;
 
             ILambda_.set
             (
-                i,
-                new volScalarField(IHeader, IDefault)
+                lambdaI,
+                new volScalarField(noReadHeader, IDefaultPtr())
             );
         }
     }
@@ -201,7 +218,6 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
         ).initialResidual();
 
         maxResidual = max(eqnResidual, maxResidual);
-
     }
 
     return maxResidual;
diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H b/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
index 27fbf1bc71..1e558c8f3c 100644
--- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
+++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
@@ -55,6 +55,13 @@ class fvDOM;
 
 class radiativeIntensityRay
 {
+public:
+
+    static const word intensityPrefix;
+
+
+private:
+
     // Private data
 
         //- Refence to the owner fvDOM object
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
index 96eb031fdd..dccbbbad25 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
@@ -203,7 +203,7 @@ Foam::radiation::greyMeanAbsorptionEmission::aCont(const label bandI) const
             {
                 Ti = 1./T[i];
             }
-            a[i]+=
+            a[i] +=
                 Yipi
                *(
                     ((((b[5]*Ti + b[4])*Ti + b[3])*Ti + b[2])*Ti + b[1])*Ti
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.H b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.H
index 7eaf936ec5..2c24af42c4 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.H
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.H
@@ -100,7 +100,7 @@ namespace radiation
 {
 
 /*---------------------------------------------------------------------------*\
-                  Class constantAbsorptionEmission Declaration
+                 Class greyMeanAbsorptionEmission Declaration
 \*---------------------------------------------------------------------------*/
 
 class greyMeanAbsorptionEmission
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
index 6a3f290706..a96deeb0b4 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
@@ -80,7 +80,7 @@ Foam::radiation::wideBandAbsorptionEmission::wideBandAbsorptionEmission
         const dictionary& dict = iter().dict();
         dict.lookup("bandLimits") >> iBands_[nBand];
         dict.lookup("EhrrCoeff") >> iEhrrCoeffs_[nBand];
-        totalWaveLength_ += (iBands_[nBand][1] - iBands_[nBand][0]);
+        totalWaveLength_ += iBands_[nBand][1] - iBands_[nBand][0];
 
         label nSpec = 0;
         const dictionary& specDicts = dict.subDict("species");
@@ -111,7 +111,7 @@ Foam::radiation::wideBandAbsorptionEmission::wideBandAbsorptionEmission
     nBands_ = nBand;
 
     // Check that all the species on the dictionary are present in the
-    // look-up table  and save the corresponding indexes of the look-up table
+    // look-up table  and save the corresponding indices of the look-up table
 
     label j = 0;
     forAllConstIter(HashTable<label>, speciesNames_, iter)
@@ -180,7 +180,7 @@ Foam::radiation::wideBandAbsorptionEmission::aCont(const label bandI) const
                 IOobject::NO_WRITE
             ),
             mesh(),
-            dimensionedScalar("a",dimless/dimLength, 0.0)
+            dimensionedScalar("a", dimless/dimLength, 0.0)
         )
     );
 
@@ -193,7 +193,7 @@ Foam::radiation::wideBandAbsorptionEmission::aCont(const label bandI) const
         for (label n=0; n<nSpecies; n++)
         {
             label l = 0;
-            scalar Yipi = 0;
+            scalar Yipi = 0.0;
             if (specieIndex_[n] != 0)
             {
                 // moles x pressure [atm]
@@ -213,7 +213,7 @@ Foam::radiation::wideBandAbsorptionEmission::aCont(const label bandI) const
 
             if (coeffs_[n][bandI].invTemp())
             {
-                Ti = 1./T[i];
+                Ti = 1.0/T[i];
             }
 
             a[i]+=
@@ -245,7 +245,7 @@ Foam::radiation::wideBandAbsorptionEmission::eCont(const label bandI) const
                 IOobject::NO_WRITE
             ),
             mesh(),
-            dimensionedScalar("e",dimless/dimLength, 0.0)
+            dimensionedScalar("e", dimless/dimLength, 0.0)
         )
     );
 
@@ -275,7 +275,7 @@ Foam::radiation::wideBandAbsorptionEmission::ECont(const label bandI) const
 
     if (mesh().foundObject<volScalarField>("hrr"))
     {
-        const volScalarField& hrr =  mesh().lookupObject<volScalarField>("hrr");
+        const volScalarField& hrr = mesh().lookupObject<volScalarField>("hrr");
         E().internalField() =
             iEhrrCoeffs_[bandI]
            *hrr.internalField()
@@ -305,7 +305,7 @@ void Foam::radiation::wideBandAbsorptionEmission::correct
 {
     a = dimensionedScalar("zero", dimless/dimLength, 0.0);
 
-    for (label j = 0; j < nBands_; j++)
+    for (label j=0; j<nBands_; j++)
     {
         Info<< "Calculating absorption in band: " << j << endl;
         aLambda[j].internalField() = this->a(j);
@@ -317,4 +317,6 @@ void Foam::radiation::wideBandAbsorptionEmission::correct
     }
 
 }
+
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.H b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.H
index d0763595ea..21e10c400b 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.H
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.H
@@ -36,9 +36,6 @@ Description
     There is no check of continuity of the bands. They should not ovelap or
     have gaps.
 
-    The black body emission power table(constant/blackBodyEmissivePower) range
-    of lambda * T = [1000; 10000] x 10E-6 (90% of the total emission).
-
     The emission constant proportionality is specified per band (EhrrCoeff).
 
     The coefficients for the species in the LookUpTable have to be specified
@@ -153,13 +150,13 @@ private:
         HashTable<label> speciesNames_;
 
         //- Indices of species in the look-up table
-        FixedList<label,nSpecies_> specieIndex_;
+        FixedList<label, nSpecies_> specieIndex_;
 
         //- Bands
-        FixedList<Vector2D<scalar>, maxBands_ > iBands_;
+        FixedList<Vector2D<scalar>, maxBands_> iBands_;
 
         //- Proportion of the heat released rate emitted
-        FixedList<scalar, maxBands_ > iEhrrCoeffs_;
+        FixedList<scalar, maxBands_> iEhrrCoeffs_;
 
         //- Look-up table of species related to ft
         mutable interpolationLookUpTable<scalar> lookUpTable_;
@@ -168,7 +165,7 @@ private:
         const basicThermo& thermo_;
 
         //- Bands
-        label nBands_ ;
+        label nBands_;
 
         //- Pointer list of species being solved involved in the absorption
         UPtrList<volScalarField> Yj_;
diff --git a/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C b/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C
index 2f3804afde..dedf725e24 100644
--- a/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C
+++ b/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C
@@ -343,7 +343,7 @@ bool Foam::movingConeTopoFvMesh::update()
             //    << endl;
 
             {
-                OFstream str(thisDb().path()/"meshPoints.obj");
+                OFstream str(time().timePath()/"meshPoints.obj");
                 Pout<< "Writing mesh with meshPoints to " << str.name()
                     << endl;
 
@@ -361,7 +361,7 @@ bool Foam::movingConeTopoFvMesh::update()
                 }
             }
             {
-                OFstream str(thisDb().path()/"preMotionPoints.obj");
+                OFstream str(time().timePath()/"preMotionPoints.obj");
                 Pout<< "Writing mesh with preMotionPoints to " << str.name()
                     << endl;
 
diff --git a/tutorials/molecularDynamics/gnemdFoam/Allclean b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/Allclean
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/Allclean
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/Allclean
diff --git a/tutorials/molecularDynamics/gnemdFoam/Allrun b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/Allrun
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/Allrun
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/Allrun
diff --git a/tutorials/molecularDynamics/gnemdFoam/constrictedChannel/constant/polyMesh/blockMeshDict b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/constant/polyMesh/blockMeshDict
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/constrictedChannel/constant/polyMesh/blockMeshDict
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/constant/polyMesh/blockMeshDict
diff --git a/tutorials/molecularDynamics/gnemdFoam/constrictedChannel/constant/polyMesh/boundary b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/constant/polyMesh/boundary
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/constrictedChannel/constant/polyMesh/boundary
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/constant/polyMesh/boundary
diff --git a/tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/controlDict b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/controlDict
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/controlDict
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/controlDict
diff --git a/tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/fvSchemes b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/fvSchemes
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/fvSchemes
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/fvSchemes
diff --git a/tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/fvSolution b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/fvSolution
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/fvSolution
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/fvSolution
diff --git a/tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/mdSolution b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/mdSolution
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/mdSolution
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/mdSolution
diff --git a/tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/molConfigDict b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/molConfigDict
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/molConfigDict
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/molConfigDict
diff --git a/tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/potentialDict b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/potentialDict
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/constrictedChannel/system/potentialDict
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/constrictedChannel/system/potentialDict
diff --git a/tutorials/molecularDynamics/gnemdFoam/nanoNozzle/constant/polyMesh/blockMeshDict b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/constant/polyMesh/blockMeshDict
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/nanoNozzle/constant/polyMesh/blockMeshDict
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/constant/polyMesh/blockMeshDict
diff --git a/tutorials/molecularDynamics/gnemdFoam/nanoNozzle/constant/polyMesh/boundary b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/constant/polyMesh/boundary
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/nanoNozzle/constant/polyMesh/boundary
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/constant/polyMesh/boundary
diff --git a/tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/controlDict b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/controlDict
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/controlDict
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/controlDict
diff --git a/tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/decomposeParDict b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/decomposeParDict
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/decomposeParDict
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/decomposeParDict
diff --git a/tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/fvSchemes b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/fvSchemes
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/fvSchemes
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/fvSchemes
diff --git a/tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/fvSolution b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/fvSolution
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/fvSolution
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/fvSolution
diff --git a/tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/mdSolution b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/mdSolution
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/mdSolution
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/mdSolution
diff --git a/tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/molConfigDict b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/molConfigDict
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/molConfigDict
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/molConfigDict
diff --git a/tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/potentialDict b/tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/potentialDict
similarity index 100%
rename from tutorials/molecularDynamics/gnemdFoam/nanoNozzle/system/potentialDict
rename to tutorials/discreteMethods/molecularDynamics/gnemdFoam/nanoNozzle/system/potentialDict
diff --git a/tutorials/molecularDynamics/mdEquilibrationFoam/Allclean b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/Allclean
similarity index 100%
rename from tutorials/molecularDynamics/mdEquilibrationFoam/Allclean
rename to tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/Allclean
diff --git a/tutorials/molecularDynamics/mdEquilibrationFoam/Allrun b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/Allrun
similarity index 100%
rename from tutorials/molecularDynamics/mdEquilibrationFoam/Allrun
rename to tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/Allrun
diff --git a/tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/constant/polyMesh/blockMeshDict b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/constant/polyMesh/blockMeshDict
similarity index 100%
rename from tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/constant/polyMesh/blockMeshDict
rename to tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/constant/polyMesh/blockMeshDict
diff --git a/tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/constant/polyMesh/boundary b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/constant/polyMesh/boundary
similarity index 100%
rename from tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/constant/polyMesh/boundary
rename to tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/constant/polyMesh/boundary
diff --git a/tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/controlDict b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/controlDict
similarity index 100%
rename from tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/controlDict
rename to tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/controlDict
diff --git a/tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/fvSchemes b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/fvSchemes
similarity index 100%
rename from tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/fvSchemes
rename to tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/fvSchemes
diff --git a/tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/fvSolution b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/fvSolution
similarity index 100%
rename from tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/fvSolution
rename to tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/fvSolution
diff --git a/tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/mdEquilibrationDict b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/mdEquilibrationDict
similarity index 100%
rename from tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/mdEquilibrationDict
rename to tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/mdEquilibrationDict
diff --git a/tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/mdSolution b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/mdSolution
similarity index 100%
rename from tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/mdSolution
rename to tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/mdSolution
diff --git a/tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/molConfigDict b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/molConfigDict
similarity index 100%
rename from tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/molConfigDict
rename to tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/molConfigDict
diff --git a/tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/potentialDict b/tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/potentialDict
similarity index 100%
rename from tutorials/molecularDynamics/mdEquilibrationFoam/periodicCube/system/potentialDict
rename to tutorials/discreteMethods/molecularDynamics/mdEquilibrationFoam/periodicCube/system/potentialDict
diff --git a/tutorials/incompressible/pimpleFoam/t-junction/0/epsilon b/tutorials/incompressible/pimpleFoam/t-junction/0/epsilon
index 2c4cfab1c0..5b446a2d8a 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction/0/epsilon
+++ b/tutorials/incompressible/pimpleFoam/t-junction/0/epsilon
@@ -24,7 +24,7 @@ boundaryField
     inlet
     {
         type            turbulentMixingLengthDissipationRateInlet;
-        mixingLength    0.01;
+        mixingLength    0.01;       // 1cm - half channel height
         value           uniform 1;
     }
 
diff --git a/tutorials/incompressible/pimpleFoam/t-junction/0/k b/tutorials/incompressible/pimpleFoam/t-junction/0/k
index 51c0ab7cd3..ad217bc9e4 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction/0/k
+++ b/tutorials/incompressible/pimpleFoam/t-junction/0/k
@@ -24,7 +24,7 @@ boundaryField
     inlet
     {
         type            turbulentIntensityKineticEnergyInlet;
-        intensity       0.05;
+        intensity       0.05;       // 5% turbulent intensity
         value           uniform 1;
     }
 
diff --git a/tutorials/incompressible/pimpleFoam/t-junction/0/p b/tutorials/incompressible/pimpleFoam/t-junction/0/p
index b3d4646058..2446eb2f7c 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction/0/p
+++ b/tutorials/incompressible/pimpleFoam/t-junction/0/p
@@ -22,8 +22,14 @@ boundaryField
 {
     inlet      
     {
-        type            totalPressure;
-        p0              uniform 100040;
+        //type            totalPressure;
+        //p0              uniform 100040;
+
+        type            timeVaryingTotalPressure;
+        p0              100040; // only used for restarts
+        outOfBounds     clamp;
+        fileName        "$FOAM_CASE/constant/p0vsTime";
+
         U               U;
         phi             phi;
         rho             none;
diff --git a/tutorials/incompressible/pimpleFoam/t-junction/constant/p0vsTime b/tutorials/incompressible/pimpleFoam/t-junction/constant/p0vsTime
new file mode 100644
index 0000000000..a9ee6bb3e7
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/t-junction/constant/p0vsTime
@@ -0,0 +1,4 @@
+(
+    (0 100010)
+    (1 100040)
+)
diff --git a/tutorials/incompressible/pimpleFoam/t-junction/constant/turbulenceProperties b/tutorials/incompressible/pimpleFoam/t-junction/constant/turbulenceProperties
new file mode 100644
index 0000000000..a18cc53b5b
--- /dev/null
+++ b/tutorials/incompressible/pimpleFoam/t-junction/constant/turbulenceProperties
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType RASModel;
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/t-junction/system/controlDict b/tutorials/incompressible/pimpleFoam/t-junction/system/controlDict
index a5ba827eec..2d7aaa81da 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction/system/controlDict
+++ b/tutorials/incompressible/pimpleFoam/t-junction/system/controlDict
@@ -52,16 +52,24 @@ functions
     probes
     {
         type        probes;
-        name        probes;
+        // Where to load it from
         functionObjectLibs ( "libsampling.so" );
+        // Name of the directory for probe data
+        name        probes;
         probeLocations
         (
-            ( 1e-06 0 0.01 )
-            ( 0.21 -0.20999 0.01 )
-            ( 0.21 0.20999 0.01 )
-            ( 0.21 0 0.01 )
+            ( 1e-06 0 0.01 )        // at inlet
+            ( 0.21 -0.20999 0.01 )  // at outlet1
+            ( 0.21 0.20999 0.01 )   // at outlet2
+            ( 0.21 0 0.01 )         // at central block
         );
+
+        // Fields to be probed
         fields ( p U );
+
+        // Write at same frequency as fields
+        outputControl   outputTime;
+        outputInterval  1;
     }
 }
 
diff --git a/tutorials/multiphase/interFoam/laminar/damBreak/system/setFieldsDict b/tutorials/multiphase/interFoam/laminar/damBreak/system/setFieldsDict
index 21e4416d79..d719028fc5 100644
--- a/tutorials/multiphase/interFoam/laminar/damBreak/system/setFieldsDict
+++ b/tutorials/multiphase/interFoam/laminar/damBreak/system/setFieldsDict
@@ -15,9 +15,23 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-defaultFieldValues ( volScalarFieldValue alpha1 0 volVectorFieldValue U ( 0 0 0 ) );
+defaultFieldValues
+(
+    volScalarFieldValue alpha1 0
+    volVectorFieldValue U (0 0 0)
+);
 
-regions         ( boxToCell { box ( 0 0 -1 ) ( 0.1461 0.292 1 ) ; fieldValues ( volScalarFieldValue alpha1 1 ) ; } );
+regions
+(
+    boxToCell
+    {
+        box (0 0 -1) (0.1461 0.292 1);
+        fieldValues
+        (
+            volScalarFieldValue alpha1 1
+        );
+    }
+);
 
 
 // ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/ras/damBreak/system/setFieldsDict b/tutorials/multiphase/interFoam/ras/damBreak/system/setFieldsDict
index bafb3a32c0..d719028fc5 100644
--- a/tutorials/multiphase/interFoam/ras/damBreak/system/setFieldsDict
+++ b/tutorials/multiphase/interFoam/ras/damBreak/system/setFieldsDict
@@ -15,9 +15,23 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-defaultFieldValues ( volScalarFieldValue gamma 0 volVectorFieldValue U ( 0 0 0 ) );
+defaultFieldValues
+(
+    volScalarFieldValue alpha1 0
+    volVectorFieldValue U (0 0 0)
+);
 
-regions         ( boxToCell { box ( 0 0 -1 ) ( 0.1461 0.292 1 ) ; fieldValues ( volScalarFieldValue gamma 1 ) ; } );
+regions
+(
+    boxToCell
+    {
+        box (0 0 -1) (0.1461 0.292 1);
+        fieldValues
+        (
+            volScalarFieldValue alpha1 1
+        );
+    }
+);
 
 
 // ************************************************************************* //