diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index 5270e3d8..cb732902 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -40,8 +40,6 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        fvOptions.constrain(pEqn);
-
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -74,8 +72,6 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        fvOptions.constrain(pEqn);
-
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
index e633ebee..fcfc2aa7 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,9 @@ License
 Application
     rhoPimpleFoam
 
+Group
+    grpCompressibleSolvers grpMovingMeshSolvers
+
 Description
     Transient solver for laminar or turbulent flow of compressible fluids
     for HVAC and similar applications.
@@ -71,17 +74,21 @@ int main(int argc, char *argv[])
 
         #include "setDeltaT.H"
 
-        runTime++;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
         {
+            // Store divrhoU from the previous time-step/mesh for the correctPhi
+            volScalarField divrhoU
+            (
+                "divrhoU",
+                fvc::div(fvc::absolute(phi, rho, U))
+            );
+
+            runTime++;
+
+            Info<< "Time = " << runTime.timeName() << nl << endl;
+
             // Store momentum to set rhoUf for introduced faces.
             volVectorField rhoU("rhoU", rho*U);
 
-            // Store divrhoU from the previous time-step/mesh for the correctPhi
-            volScalarField divrhoU(fvc::div(fvc::absolute(phi, rho, U)));
-
             // Do any mesh changes
             mesh.update();
 
@@ -102,12 +109,9 @@ int main(int argc, char *argv[])
             #include "meshCourantNo.H"
         }
 
-        if (pimple.nCorrPIMPLE() <= 1)
-        {
-            #include "rhoEqn.H"
-            Info<< "rhoEqn max/min : " << max(rho).value()
-                << " " << min(rho).value() << endl;
-        }
+        #include "rhoEqn.H"
+        Info<< "rhoEqn max/min : " << max(rho).value()
+            << " " << min(rho).value() << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
index f2c091e4..e1387379 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,6 +47,7 @@ Description
 #include "fvMesh.H"
 #include "MeshedSurfaces.H"
 #include "globalIndex.H"
+#include "cellSet.H"
 
 #include "extrudedMesh.H"
 #include "extrudeModel.H"
@@ -192,6 +193,23 @@ void updateFaceLabels(const mapPolyMesh& map, labelList& faceLabels)
 }
 
 
+void updateCellSet(const mapPolyMesh& map, labelHashSet& cellLabels)
+{
+    const labelList& reverseMap = map.reverseCellMap();
+
+    labelHashSet newCellLabels(2*cellLabels.size());
+
+    forAll(cellLabels, i)
+    {
+        label oldCellI = cellLabels[i];
+
+        if (reverseMap[oldCellI] >= 0)
+        {
+            newCellLabels.insert(reverseMap[oldCellI]);
+        }
+    }
+    cellLabels.transfer(newCellLabels);
+}
 
 
 int main(int argc, char *argv[])
@@ -270,6 +288,9 @@ int main(int argc, char *argv[])
     word backPatchName;
     labelList backPatchFaces;
 
+    // Optional added cells (get written to cellSet)
+    labelHashSet addedCellsSet;
+
     if (mode == PATCH || mode == MESH)
     {
         if (flipNormals && mode == MESH)
@@ -670,11 +691,33 @@ int main(int argc, char *argv[])
             map().reverseFaceMap(),
             backPatchFaces
         );
+
+        // Store added cells
+        if (mode == MESH)
+        {
+            const labelListList addedCells
+            (
+                layerExtrude.addedCells
+                (
+                    meshFromMesh,
+                    layerExtrude.layerFaces()
+                )
+            );
+            forAll(addedCells, faceI)
+            {
+                const labelList& aCells = addedCells[faceI];
+                forAll(aCells, i)
+                {
+                    addedCellsSet.insert(aCells[i]);
+                }
+            }
+        }
     }
     else
     {
         // Read from surface
         fileName surfName(dict.lookup("surface"));
+        surfName.expand();
 
         Info<< "Extruding surfaceMesh read from file " << surfName << nl
             << endl;
@@ -810,6 +853,7 @@ int main(int argc, char *argv[])
             // Update stored data
             updateFaceLabels(map(), frontPatchFaces);
             updateFaceLabels(map(), backPatchFaces);
+            updateCellSet(map(), addedCellsSet);
 
             // Move mesh (if inflation used)
             if (map().hasMotionPoints())
@@ -913,6 +957,9 @@ int main(int argc, char *argv[])
         // Update fields
         mesh.updateMesh(map);
 
+        // Update local data
+        updateCellSet(map(), addedCellsSet);
+
         // Move mesh (if inflation used)
         if (map().hasMotionPoints())
         {
@@ -929,6 +976,21 @@ int main(int argc, char *argv[])
             << exit(FatalError);
     }
 
+    // Need writing cellSet
+    label nAdded = returnReduce(addedCellsSet.size(), sumOp<label>());
+    if (nAdded > 0)
+    {
+        cellSet addedCells(mesh, "addedCells", addedCellsSet);
+        Info<< "Writing added cells to cellSet " << addedCells.name()
+            << nl << endl;
+        if (!addedCells.write())
+        {
+            FatalErrorIn(args.executable()) << "Failed writing cellSet"
+                << addedCells.name()
+                << exit(FatalError);
+        }
+    }
+
     Info<< "End\n" << endl;
 
     return 0;
diff --git a/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectList.H b/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectList.H
index 56121a22..4409dac8 100644
--- a/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectList.H
+++ b/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectList.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -128,6 +128,21 @@ public:
         //- Type of values the UList contains.
         typedef T value_type;
 
+        //- Type that can be used for storing into
+        //  UList::value_type objects.
+        typedef T& reference;
+
+        //- Type that can be used for storing into
+        //  constant UList::value_type objects
+        typedef const T& const_reference;
+
+        //- The type that can represent the difference between any two
+        //  UList iterator objects.
+        typedef label difference_type;
+
+        //- The type that can represent the size of a UList.
+        typedef label size_type;
+
 
     // Ostream operator
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
index 643e06a4..d7508c3d 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,8 +61,8 @@ void Foam::polyPatch::movePoints(PstreamBuffers&, const pointField& p)
 
 void Foam::polyPatch::updateMesh(PstreamBuffers&)
 {
+    primitivePatch::clearGeom();
     clearAddressing();
-    primitivePatch::clearOut();
 }
 
 
@@ -348,6 +348,8 @@ const Foam::labelList& Foam::polyPatch::meshEdges() const
 
 void Foam::polyPatch::clearAddressing()
 {
+    primitivePatch::clearTopology();
+    primitivePatch::clearPatchMeshAddr();
     deleteDemandDrivenData(faceCellsPtr_);
     deleteDemandDrivenData(mePtr_);
 }
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchClear.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchClear.C
index 664e5a19..d4aa873d 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchClear.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchClear.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,6 +92,7 @@ clearTopology()
     deleteDemandDrivenData(pointEdgesPtr_);
     deleteDemandDrivenData(pointFacesPtr_);
     deleteDemandDrivenData(edgeLoopsPtr_);
+    deleteDemandDrivenData(localPointOrderPtr_);
 }
 
 
diff --git a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
index d0d3b30b..7e68450a 100644
--- a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
+++ b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
@@ -150,7 +150,7 @@ Foam::vector Foam::eigenValues(const tensor& t)
             if (R2 < Q3)
             {
                 scalar sqrtQ = sqrt(Q);
-                scalar theta = acos(R/(Q*sqrtQ));
+                scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ))));
 
                 scalar m2SqrtQ = -2*sqrtQ;
                 scalar aBy3 = a/3;
@@ -344,7 +344,7 @@ Foam::vector Foam::eigenValues(const symmTensor& t)
             if (R2 < Q3)
             {
                 scalar sqrtQ = sqrt(Q);
-                scalar theta = acos(R/(Q*sqrtQ));
+                scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ))));
 
                 scalar m2SqrtQ = -2*sqrtQ;
                 scalar aBy3 = a/3;
diff --git a/src/Pstream/mpi/UPstream.C b/src/Pstream/mpi/UPstream.C
index 72f1e4b0..610c6fa9 100644
--- a/src/Pstream/mpi/UPstream.C
+++ b/src/Pstream/mpi/UPstream.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -286,7 +286,7 @@ void Foam::reduce
         Pout<< "UPstream::allocateRequest for non-blocking reduce"
             << " : request:" << requestID
             << endl;
-
+    }
 #else
     // Non-blocking not yet implemented in mpi
     reduce(Value, bop, tag, communicator);
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
index 8607d49f..2ae2b910 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
@@ -843,7 +843,9 @@ tmp<surfaceScalarField> CoEulerDdtScheme<Type>::meshPhi
             (
                 "meshPhi",
                 mesh().time().timeName(),
-                mesh()
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
             ),
             mesh(),
             dimensionedScalar("0", dimVolume/dimTime, 0.0)
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
index c9b095fa..d866f8c5 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
@@ -978,7 +978,7 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
 
         fvm.source() =
         (
-            rDtCoef*rho.internalField()*vf.oldTime().internalField()
+            rDtCoef*rho.oldTime().internalField()*vf.oldTime().internalField()
           + offCentre_(ddt0.internalField())
         )*mesh().V0();
     }
@@ -1079,8 +1079,8 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
         fvm.source() =
         (
             rDtCoef
-           *alpha.internalField()
-           *rho.internalField()
+           *alpha.oldTime().internalField()
+           *rho.oldTime().internalField()
            *vf.oldTime().internalField()
           + offCentre_(ddt0.internalField())
         )*mesh().V0();
@@ -1426,7 +1426,22 @@ tmp<surfaceScalarField> CrankNicolsonDdtScheme<Type>::meshPhi
             coef0_(meshPhi0)*mesh().phi().oldTime() - offCentre_(meshPhi0());
     }
 
-    return coef_(meshPhi0)*mesh().phi() - offCentre_(meshPhi0());
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                mesh().phi().name(),
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            coef_(meshPhi0)*mesh().phi() - offCentre_(meshPhi0())
+        )
+    );
 }
 
 
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
index 955bba08..a53301f5 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
@@ -845,7 +845,10 @@ tmp<surfaceScalarField> SLTSDdtScheme<Type>::meshPhi
             (
                 "meshPhi",
                 mesh().time().timeName(),
-                mesh()
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
             ),
             mesh(),
             dimensionedScalar("0", dimVolume/dimTime, 0.0)
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
index 1474b8c8..3490fd38 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
@@ -976,7 +976,22 @@ tmp<surfaceScalarField> backwardDdtScheme<Type>::meshPhi
     // Coefficient for t-1/2 (between times n and 0)
     scalar coefftn_0 = 1 + coefft0_00;
 
-    return coefftn_0*mesh().phi() - coefft0_00*mesh().phi().oldTime();
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                mesh().phi().name(),
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            coefftn_0*mesh().phi() - coefft0_00*mesh().phi().oldTime()
+        )
+    );
 }
 
 
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
index e1c0ea5e..f308714e 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C
@@ -740,7 +740,10 @@ tmp<surfaceScalarField> localEulerDdtScheme<Type>::meshPhi
             (
                 "meshPhi",
                 mesh().time().timeName(),
-                mesh()
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
             ),
             mesh(),
             dimensionedScalar("0", dimVolume/dimTime, 0.0)
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C
index fb695f1b..97f4b40a 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C
@@ -412,7 +412,10 @@ tmp<surfaceScalarField> steadyStateDdtScheme<Type>::meshPhi
             (
                 "meshPhi",
                 mesh().time().timeName(),
-                mesh()
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
             ),
             mesh(),
             dimensionedScalar("0", dimVolume/dimTime, 0.0)
diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.C
index ddc60991..f5ce1243 100644
--- a/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.C
+++ b/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.C
@@ -46,15 +46,6 @@ wedgeFvPatch::wedgeFvPatch(const polyPatch& patch, const fvBoundaryMesh& bm)
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::vectorField> Foam::wedgeFvPatch::delta() const
-{
-    const vectorField nHat(nf());
-    return nHat*(nHat & (Cf() - Cn()));
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.H
index 16fab3e4..25be556b 100644
--- a/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.H
+++ b/src/finiteVolume/fvMesh/fvPatches/constraint/wedge/wedgeFvPatch.H
@@ -83,9 +83,6 @@ public:
             {
                 return wedgePolyPatch_.cellT();
             }
-
-            //- Return cell-centre to face normal vector
-            virtual tmp<vectorField> delta() const;
 };
 
 
diff --git a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C
index c55ad272..e7dd4120 100644
--- a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C
+++ b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.C
@@ -140,7 +140,9 @@ const Foam::scalarField& Foam::fvPatch::magSf() const
 
 Foam::tmp<Foam::vectorField> Foam::fvPatch::delta() const
 {
-    return Cf() - Cn();
+    // Use patch-normal delta for all non-coupled BCs
+    const vectorField nHat(nf());
+    return nHat*(nHat & (Cf() - Cn()));
 }
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C b/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C
index 61424852..a883c39d 100644
--- a/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C
+++ b/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C
@@ -103,8 +103,12 @@ bool Foam::trackedParticle::move
     scalar tEnd = (1.0 - stepFraction())*trackTime;
     scalar dtMax = tEnd;
 
-    if (tEnd <= SMALL)
+    if (tEnd <= SMALL && onBoundary())
     {
+        // This is a hack to handle particles reaching their endpoint
+        // on a processor boundary. If the endpoint is on a processor face
+        // it currently gets transferred backwards and forwards infinitely.
+
         // Remove the particle
         td.keepParticle = false;
     }
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
index 70673d31..b63e0f7a 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
@@ -30,7 +30,7 @@ License
 template<class SourcePatch, class TargetPatch>
 void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
 (
-    boolList& mapFlag,
+    labelList& mapFlag,
     labelList& srcTgtSeed,
     DynamicList<label>& srcSeeds,
     DynamicList<label>& nonOverlapFaces,
@@ -50,19 +50,32 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
     {
         label srcI = srcNbr[i];
 
-        if (mapFlag[srcI] && srcTgtSeed[srcI] == -1)
+        if ((mapFlag[srcI] == 0) && (srcTgtSeed[srcI] == -1))
         {
+            // first attempt: match by comparing face centres
             const face& srcF = this->srcPatch_[srcI];
-            const vector srcN = srcF.normal(srcPoints);
+            const point& srcC = srcCf[srcI];
+
+            scalar tol = GREAT;
+            forAll(srcF, fpI)
+            {
+                const point& p = srcPoints[srcF[fpI]];
+                scalar d2 = magSqr(p - srcC);
+                if (d2 < tol)
+                {
+                    tol = d2;
+                }
+            }
+            tol = max(SMALL, 0.0001*sqrt(tol));
 
             bool found = false;
             forAll(tgtNbr, j)
             {
                 label tgtI = tgtNbr[j];
                 const face& tgtF = this->tgtPatch_[tgtI];
-                pointHit ray = tgtF.ray(srcCf[srcI], srcN, tgtPoints);
+                const point tgtC = tgtF.centre(tgtPoints);
 
-                if (ray.hit())
+                if (mag(srcC - tgtC) < tol)
                 {
                     // new match - append to lists
                     found = true;
@@ -74,11 +87,57 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
                 }
             }
 
+            // second attempt: match by shooting a ray into the tgt face
             if (!found)
             {
-                // no match available for source face srcI
-                mapFlag[srcI] = false;
+                const vector srcN = srcF.normal(srcPoints);
+
+                forAll(tgtNbr, j)
+                {
+                    label tgtI = tgtNbr[j];
+                    const face& tgtF = this->tgtPatch_[tgtI];
+                    pointHit ray = tgtF.ray(srcCf[srcI], srcN, tgtPoints);
+
+                    if (ray.hit())
+                    {
+                        // new match - append to lists
+                        found = true;
+
+                        srcTgtSeed[srcI] = tgtI;
+                        srcSeeds.append(srcI);
+
+                        break;
+                    }
+                }
+            }
+
+            // no match available for source face srcI
+            if (!found)
+            {
+                mapFlag[srcI] = -1;
                 nonOverlapFaces.append(srcI);
+
+                if (debug)
+                {
+                    Pout<< "source face not found: id=" << srcI
+                        << " centre=" << srcCf[srcI]
+                        << " normal=" << srcF.normal(srcPoints)
+                        << " points=" << srcF.points(srcPoints)
+                        << endl;
+
+                    Pout<< "target neighbours:" << nl;
+                    forAll(tgtNbr, j)
+                    {
+                        label tgtI = tgtNbr[j];
+                        const face& tgtF = this->tgtPatch_[tgtI];
+
+                        Pout<< "face id: " << tgtI
+                            << " centre=" << tgtF.centre(tgtPoints)
+                            << " normal=" << tgtF.normal(tgtPoints)
+                            << " points=" << tgtF.points(tgtPoints)
+                            << endl;
+                    }
+                }
             }
         }
     }
@@ -96,6 +155,36 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
 }
 
 
+template<class SourcePatch, class TargetPatch>
+void Foam::directAMI<SourcePatch, TargetPatch>::restartAdvancingFront
+(
+    labelList& mapFlag,
+    DynamicList<label>& nonOverlapFaces,
+    label& srcFaceI,
+    label& tgtFaceI
+) const
+{
+    forAll(mapFlag, faceI)
+    {
+        if (mapFlag[faceI] == 0)
+        {
+            tgtFaceI = this->findTargetFace(faceI);
+
+            if (tgtFaceI < 0)
+            {
+                mapFlag[faceI] = -1;
+                nonOverlapFaces.append(faceI);
+            }
+            else
+            {
+                srcFaceI = faceI;
+                break;
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class SourcePatch, class TargetPatch>
@@ -176,15 +265,19 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
     srcTgtSeed[srcFaceI] = tgtFaceI;
 
     // list to keep track of whether src face can be mapped
-    boolList mapFlag(srcAddr.size(), true);
+    // 1 = mapped, 0 = untested, -1 = cannot map
+    labelList mapFlag(srcAddr.size(), 0);
 
+    label nTested = 0;
     DynamicList<label> nonOverlapFaces;
     do
     {
         srcAddr[srcFaceI].append(tgtFaceI);
         tgtAddr[tgtFaceI].append(srcFaceI);
 
-        mapFlag[srcFaceI] = false;
+        mapFlag[srcFaceI] = 1;
+
+        nTested++;
 
         // Do advancing front starting from srcFaceI, tgtFaceI
         appendToDirectSeeds
@@ -196,6 +289,12 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
             srcFaceI,
             tgtFaceI
         );
+
+        if (srcFaceI < 0 && nTested < this->srcPatch_.size())
+        {
+            restartAdvancingFront(mapFlag, nonOverlapFaces, srcFaceI, tgtFaceI);
+        }
+
     } while (srcFaceI >= 0);
 
     if (nonOverlapFaces.size() != 0)
@@ -211,14 +310,16 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
     forAll(srcAddr, i)
     {
         scalar magSf = this->srcMagSf_[i];
-        srcAddress[i].transfer(srcAddr[i]);
+//        srcWeights[i] = scalarList(srcAddr[i].size(), magSf);
         srcWeights[i] = scalarList(1, magSf);
+        srcAddress[i].transfer(srcAddr[i]);
     }
     forAll(tgtAddr, i)
     {
         scalar magSf = this->tgtMagSf_[i];
-        tgtAddress[i].transfer(tgtAddr[i]);
+//        tgtWeights[i] = scalarList(tgtAddr[i].size(), magSf);
         tgtWeights[i] = scalarList(1, magSf);
+        tgtAddress[i].transfer(tgtAddr[i]);
     }
 }
 
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H
index dc273c79..d74094d0 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H
@@ -67,7 +67,7 @@ private:
             //- Append to list of src face seed indices
             void appendToDirectSeeds
             (
-                boolList& mapFlag,
+                labelList& mapFlag,
                 labelList& srcTgtSeed,
                 DynamicList<label>& srcSeeds,
                 DynamicList<label>& nonOverlapFaces,
@@ -75,6 +75,16 @@ private:
                 label& tgtFaceI
             ) const;
 
+            //- Restart the advancing front - typically happens for
+            //  disconnected regions
+            void restartAdvancingFront
+            (
+                labelList& mapFlag,
+                DynamicList<label>& nonOverlapFaces,
+                label& srcFaceI,
+                label& tgtFaceI
+            ) const;
+
 
         // Evaluation
 
diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
index ba1e110c..5a862378 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
@@ -99,11 +99,10 @@ void Foam::fieldValues::faceSource::setFaceZoneFaces()
 
     const faceZone& fZone = mesh().faceZones()[zoneId];
 
-    faceId_.setSize(fZone.size());
-    facePatchId_.setSize(fZone.size());
-    faceSign_.setSize(fZone.size());
+    DynamicList<label> faceIds(fZone.size());
+    DynamicList<label> facePatchIds(fZone.size());
+    DynamicList<label> faceSigns(fZone.size());
 
-    label count = 0;
     forAll(fZone, i)
     {
         label faceI = fZone[i];
@@ -145,27 +144,26 @@ void Foam::fieldValues::faceSource::setFaceZoneFaces()
         {
             if (fZone.flipMap()[i])
             {
-                faceSign_[count] = -1;
+                faceSigns.append(-1);
             }
             else
             {
-                faceSign_[count] = 1;
+                faceSigns.append(1);
             }
-            faceId_[count] = faceId;
-            facePatchId_[count] = facePatchId;
-            count++;
+            faceIds.append(faceId);
+            facePatchIds.append(facePatchId);
         }
     }
 
-    faceId_.setSize(count);
-    facePatchId_.setSize(count);
-    faceSign_.setSize(count);
+    faceId_.transfer(faceIds);
+    facePatchId_.transfer(facePatchIds);
+    faceSign_.transfer(faceSigns);
     nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
 
     if (debug)
     {
         Pout<< "Original face zone size = " << fZone.size()
-            << ", new size = " << count << endl;
+            << ", new size = " << faceId_.size() << endl;
     }
 }
 
@@ -445,12 +443,43 @@ void Foam::fieldValues::faceSource::initialise(const dictionary& dict)
 
     if (dict.readIfPresent("weightField", weightFieldName_))
     {
-        Info<< "    weight field = " << weightFieldName_;
+        Info<< "    weight field = " << weightFieldName_ << nl;
+    }
+
+    if (dict.found("orientedWeightField"))
+    {
+        if (weightFieldName_ == "none")
+        {
+            dict.lookup("orientedWeightField") >>  weightFieldName_;
+            Info<< "    weight field = " << weightFieldName_ << nl;
+            orientWeightField_ = true;
+        }
+        else
+        {
+            FatalIOErrorIn
+            (
+                "void Foam::fieldValues::faceSource::initialise"
+                "("
+                    "const dictionary&"
+                ")",
+                dict
+            )
+                << "Either weightField or orientedWeightField can be supplied, "
+                << "but not both"
+                << exit(FatalIOError);
+        }
+    }
+
+    List<word> orientedFields;
+    if (dict.readIfPresent("orientedFields", orientedFields))
+    {
+        orientedFieldsStart_ = fields_.size();
+        fields_.append(orientedFields);
     }
 
     if (dict.readIfPresent("scaleFactor", scaleFactor_))
     {
-        Info<< "    scale factor = " << scaleFactor_;
+        Info<< "    scale factor = " << scaleFactor_ << nl;
     }
 
     Info<< nl << endl;
@@ -581,6 +610,8 @@ Foam::fieldValues::faceSource::faceSource
     source_(sourceTypeNames_.read(dict.lookup("source"))),
     operation_(operationTypeNames_.read(dict.lookup("operation"))),
     weightFieldName_("none"),
+    orientWeightField_(false),
+    orientedFieldsStart_(labelMax),
     scaleFactor_(1.0),
     nFaces_(0),
     faceId_(),
@@ -628,24 +659,42 @@ void Foam::fieldValues::faceSource::write()
             totalArea = gSum(filterField(mesh().magSf(), false));
         }
 
-
         if (Pstream::master())
         {
             file() << obr_.time().value() << tab << totalArea;
         }
 
+        // construct weight field
+        scalarField weightField(faceId_.size(), 1.0);
+        if (weightFieldName_ != "none")
+        {
+            weightField =
+                getFieldValues<scalar>
+                (
+                    weightFieldName_,
+                    true,
+                    orientWeightField_
+                );
+        }
+
+        // Combine onto master
+        combineFields(weightField);
+
+        // process the fields
         forAll(fields_, i)
         {
             const word& fieldName = fields_[i];
-            bool processed = false;
+            bool ok = false;
 
-            processed = processed || writeValues<scalar>(fieldName);
-            processed = processed || writeValues<vector>(fieldName);
-            processed = processed || writeValues<sphericalTensor>(fieldName);
-            processed = processed || writeValues<symmTensor>(fieldName);
-            processed = processed || writeValues<tensor>(fieldName);
+            bool orient = i >= orientedFieldsStart_;
+            ok = ok || writeValues<scalar>(fieldName, weightField, orient);
+            ok = ok || writeValues<vector>(fieldName, weightField, orient);
+            ok = ok
+              || writeValues<sphericalTensor>(fieldName, weightField, orient);
+            ok = ok || writeValues<symmTensor>(fieldName, weightField, orient);
+            ok = ok || writeValues<tensor>(fieldName, weightField, orient);
 
-            if (!processed)
+            if (!ok)
             {
                 WarningIn("void Foam::fieldValues::faceSource::write()")
                     << "Requested field " << fieldName
diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
index d58f1d9a..9589304c 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,8 +71,10 @@ Description
         sourceName   | name of face source if required  | no |
         operation    | operation to perform    | yes         |
         weightField  | name of field to apply weighting | no |
+        orientedWeightField  | name of oriented field to apply weighting | no |
         scaleFactor  | scale factor            | no          | 1
         fields       | list of fields to operate on | yes    |
+        orientedFields | list of oriented fields to operate on | no |
     \endtable
 
     \linebreak
@@ -109,8 +111,8 @@ Note
     - faces on empty patches get ignored
     - if the field is a volField the \c faceZone can only consist of boundary
       faces
-    - all fields are oriented according to the \c faceZone (so you might
-      e.g. see negative pressure)
+    - the `oriented' entries relate to mesh-oriented fields, such as the
+      flux, phi.  These fields will be oriented according to the face normals.
     - using \c sampledSurfaces:
         - not available for surface fields
         - if interpolate=true they use \c interpolationCellPoint
@@ -128,6 +130,7 @@ SeeAlso
 
 SourceFiles
     faceSource.C
+    faceSourceTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -242,6 +245,12 @@ protected:
         //- Weight field name - optional
         word weightFieldName_;
 
+        //- Flag to indicate if flipMap should be applied to the weight field
+        bool orientWeightField_;
+
+        //- Start index of fields that require application of flipMap
+        label orientedFieldsStart_;
+
         //- Scale factor - optional
         scalar scaleFactor_;
 
@@ -282,7 +291,8 @@ protected:
         tmp<Field<Type> > getFieldValues
         (
             const word& fieldName,
-            const bool mustGet = false
+            const bool mustGet = false,
+            const bool applyOrientation = false
         ) const;
 
         //- Apply the 'operation' to the values. Operation has to
@@ -356,7 +366,12 @@ public:
 
             //- Templated helper function to output field values
             template<class Type>
-            bool writeValues(const word& fieldName);
+            bool writeValues
+            (
+                const word& fieldName,
+                const scalarField& weightField,
+                const bool orient
+            );
 
             //- Filter a surface field according to faceIds
             template<class Type>
diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
index 6def30c5..c9e69771 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
@@ -54,7 +54,8 @@ template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::getFieldValues
 (
     const word& fieldName,
-    const bool mustGet
+    const bool mustGet,
+    const bool applyOrientation
 ) const
 {
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> sf;
@@ -62,7 +63,7 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::getFieldValues
 
     if (source_ != stSampledSurface && obr_.foundObject<sf>(fieldName))
     {
-        return filterField(obr_.lookupObject<sf>(fieldName), true);
+        return filterField(obr_.lookupObject<sf>(fieldName), applyOrientation);
     }
     else if (obr_.foundObject<vf>(fieldName))
     {
@@ -103,7 +104,7 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::getFieldValues
         }
         else
         {
-            return filterField(fld, true);
+            return filterField(fld, applyOrientation);
         }
     }
 
@@ -115,6 +116,7 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::getFieldValues
             "Foam::fieldValues::faceSource::getFieldValues"
             "("
                 "const word&, "
+                "const bool, "
                 "const bool"
             ") const"
         )   << "Field " << fieldName << " not found in database"
@@ -267,22 +269,20 @@ Type Foam::fieldValues::faceSource::processValues
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-bool Foam::fieldValues::faceSource::writeValues(const word& fieldName)
+bool Foam::fieldValues::faceSource::writeValues
+(
+    const word& fieldName,
+    const scalarField& weightField,
+    const bool orient
+)
 {
     const bool ok = validField<Type>(fieldName);
 
     if (ok)
     {
-        Field<Type> values(getFieldValues<Type>(fieldName));
-        scalarField weightField(values.size(), 1.0);
-
-        if (weightFieldName_ != "none")
-        {
-            weightField = getFieldValues<scalar>(weightFieldName_, true);
-        }
+        Field<Type> values(getFieldValues<Type>(fieldName, true, orient));
 
         vectorField Sf;
-
         if (surfacePtr_.valid())
         {
             // Get oriented Sf
@@ -291,13 +291,12 @@ bool Foam::fieldValues::faceSource::writeValues(const word& fieldName)
         else
         {
             // Get oriented Sf
-            Sf = filterField(mesh().Sf(), false);
+            Sf = filterField(mesh().Sf(), true);
         }
 
         // Combine onto master
         combineFields(values);
         combineFields(Sf);
-        combineFields(weightField);
 
         // Write raw values on surface if specified
         if (surfaceWriterPtr_.valid())
@@ -378,7 +377,8 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::filterField
             (
                 "fieldValues::faceSource::filterField"
                 "("
-                    "const GeometricField<Type, fvPatchField, volMesh>&"
+                    "const GeometricField<Type, fvPatchField, volMesh>&, "
+                    "const bool"
                 ") const"
             )   << type() << " " << name_ << ": "
                 << sourceTypeNames_[source_] << "(" << sourceName_ << "):"
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
index 4da01ef2..b7266df9 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
@@ -52,7 +52,7 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
 
     os.writeKeyword("centreOfMass")
         << initialCentreOfMass_ << token::END_STATEMENT << nl;
-    os.writeKeyword("orientation")
+    os.writeKeyword("initialOrientation")
         << initialQ_ << token::END_STATEMENT << nl;
     os.writeKeyword("mass")
         << mass_ << token::END_STATEMENT << nl;
diff --git a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.C b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.C
index e2b99d08..c3fefab7 100644
--- a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.C
+++ b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.C
@@ -33,6 +33,8 @@ Foam::PengRobinsonGas<Specie>::PengRobinsonGas(Istream& is)
 :
     Specie(is),
     Tc_(readScalar(is)),
+    Vc_(readScalar(is)),
+    Zc_(readScalar(is)),
     Pc_(readScalar(is)),
     omega_(readScalar(is))
 {
@@ -48,9 +50,13 @@ Foam::PengRobinsonGas<Specie>::PengRobinsonGas
 :
     Specie(dict),
     Tc_(readScalar(dict.subDict("equationOfState").lookup("Tc"))),
+    Vc_(readScalar(dict.subDict("equationOfState").lookup("Vc"))),
+    Zc_(1.0),
     Pc_(readScalar(dict.subDict("equationOfState").lookup("Pc"))),
     omega_(readScalar(dict.subDict("equationOfState").lookup("omega")))
-{}
+{
+    Zc_ = Pc_*Vc_/(specie::RR*Tc_);
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -74,6 +80,8 @@ Foam::Ostream& Foam::operator<<
 {
     os  << static_cast<const Specie&>(pg)
         << token::SPACE << pg.Tc_
+        << token::SPACE << pg.Vc_
+        << token::SPACE << pg.Zc_
         << token::SPACE << pg.Pc_
         << token::SPACE << pg.omega_;
 
diff --git a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H
index 7d1df112..568b35a3 100644
--- a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H
+++ b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H
@@ -98,10 +98,16 @@ class PengRobinsonGas
         //- Critical Temperature [K]
         scalar Tc_;
 
+        //- Critical volume [m^3/kmol]
+        scalar Vc_;
+
+        //- Critical compression factor [-]
+        scalar Zc_;
+
         //- Critical Pressure [Pa]
         scalar Pc_;
 
-        //- Accentric factor [-]
+        //- Acentric factor [-]
         scalar omega_;
 
 
@@ -114,6 +120,8 @@ public:
         (
             const Specie& sp,
             const scalar& Tc,
+            const scalar& Vc,
+            const scalar& Zc,
             const scalar& Pc,
             const scalar& omega
         );
@@ -163,7 +171,7 @@ public:
             //- Return compressibility rho/p [s^2/m^2]
             inline scalar psi(scalar p, scalar T) const;
 
-            //- Return compression factor []
+            //- Return compression factor [-]
             inline scalar Z(scalar p, scalar T) const;
 
             //- Return (cp - cv) [J/(kmol K]
diff --git a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H
index aae0b92f..d5b0c475 100644
--- a/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H
+++ b/src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H
@@ -21,7 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-
 \*---------------------------------------------------------------------------*/
 
 #include "PengRobinsonGas.H"
@@ -34,12 +33,16 @@ inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
 (
     const Specie& sp,
     const scalar& Tc,
+    const scalar& Vc,
+    const scalar& Zc,
     const scalar& Pc,
     const scalar& omega
 )
 :
     Specie(sp),
     Tc_(Tc),
+    Vc_(Vc),
+    Zc_(Zc),
     Pc_(Pc),
     omega_(omega)
 {}
@@ -55,9 +58,11 @@ inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
 )
 :
     Specie(name, pg),
-    Tc_(pg.Tc),
-    Pc_(pg.Pc),
-    omega_(pg.omega)
+    Tc_(pg.Tc_),
+    Pc_(pg.Pc_),
+    Vc_(pg.Vc_),
+    Zc_(pg.Zc_),
+    omega_(pg.omega_)
 {}
 
 
@@ -96,6 +101,7 @@ Foam::PengRobinsonGas<Specie>::New
     );
 }
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Specie>
@@ -214,7 +220,9 @@ inline void Foam::PengRobinsonGas<Specie>::operator+=
     scalar molr2 = pg.nMoles()/this->nMoles();
 
     Tc_ = molr1*Tc_ + molr2*pg.Tc_;
-    Pc_ = molr1*Pc_ + molr2*pg.Pc_;
+    Vc_ = molr1*Vc_ + molr2*pg.Vc_;
+    Zc_ = molr1*Zc_ + molr2*pg.Zc_;
+    Pc_ = specie::RR*Zc_*Tc_/Vc_;
     omega_ = molr1*omega_ + molr2*pg.omega_;
 }
 
@@ -233,7 +241,9 @@ inline void Foam::PengRobinsonGas<Specie>::operator-=
     scalar molr2 = pg.nMoles()/this->nMoles();
 
     Tc_ = molr1*Tc_ - molr2*pg.Tc_;
-    Pc_ = molr1*Pc_ - molr2*pg.Pc_;
+    Vc_ = molr1*Vc_ - molr2*pg.Vc_;
+    Zc_ = molr1*Zc_ - molr2*pg.Zc_;
+    Pc_ = specie::RR*Zc_*Tc_/Vc_;
     omega_ = molr1*omega_ - molr2*pg.omega_;
 }
 
@@ -259,12 +269,18 @@ Foam::PengRobinsonGas<Specie> Foam::operator+
     scalar molr1 = pg1.nMoles()/nMoles;
     scalar molr2 = pg2.nMoles()/nMoles;
 
+    const scalar Tc = molr1*pg1.Tc_ + molr2*pg2.Tc_;
+    const scalar Vc = molr1*pg1.Vc_ + molr2*pg2.Vc_;
+    const scalar Zc = molr1*pg1.Zc_ + molr2*pg2.Zc_;
+
     return PengRobinsonGas<Specie>
     (
         static_cast<const Specie&>(pg1)
       + static_cast<const Specie&>(pg2),
-        molr1*pg1.Tc_ + molr2*pg2.Tc_,
-        molr1*pg1.Pc_ + molr2*pg2.Pc_,
+        Tc,
+        Vc,
+        Zc,
+        specie::RR*Zc*Tc/Vc,
         molr1*pg1.omega_ + molr2*pg2.omega_
     );
 }
@@ -281,12 +297,18 @@ Foam::PengRobinsonGas<Specie> Foam::operator-
     scalar molr1 = pg1.nMoles()/nMoles;
     scalar molr2 = pg2.nMoles()/nMoles;
 
+    const scalar Tc = molr1*pg1.Tc_ + molr2*pg2.Tc_;
+    const scalar Vc = molr1*pg1.Vc_ + molr2*pg2.Vc_;
+    const scalar Zc = molr1*pg1.Zc_ + molr2*pg2.Zc_;
+
     return PengRobinsonGas<Specie>
     (
         static_cast<const Specie&>(pg1)
       - static_cast<const Specie&>(pg2),
-        molr1*pg1.Tc_ - molr2*pg2.Tc_,
-        molr1*pg1.Pc_ - molr2*pg2.Pc_,
+        Tc,
+        Vc,
+        Zc,
+        specie::RR*Zc*Tc/Vc,
         molr1*pg1.omega_ - molr2*pg2.omega_
     );
 }
@@ -303,6 +325,8 @@ Foam::PengRobinsonGas<Specie> Foam::operator*
     (
         s*static_cast<const Specie&>(pg),
         pg.Tc_,
+        pg.Vc_,
+        pg.Zc_,
         pg.Pc_,
         pg.omega_
     );