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/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/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/direct/directMethod.C b/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/direct/directMethod.C
index f8d5598f..a67cbe4a 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/direct/directMethod.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/calcMethod/direct/directMethod.C
@@ -154,16 +154,14 @@ void Foam::directMethod::calculateAddressing
     // transfer addressing into persistent storage
     forAll(srcToTgtCellAddr, i)
     {
-        scalar v = srcVc[i];
+        srcToTgtCellWght[i] = scalarList(srcToTgt[i].size(), srcVc[i]);
         srcToTgtCellAddr[i].transfer(srcToTgt[i]);
-        srcToTgtCellWght[i] = scalarList(1, v);
     }
 
     forAll(tgtToSrcCellAddr, i)
     {
-        scalar v = tgtVc[i];
+        tgtToSrcCellWght[i] = scalarList(tgtToSrc[i].size(), tgtVc[i]);
         tgtToSrcCellAddr[i].transfer(tgtToSrc[i]);
-        tgtToSrcCellWght[i] = scalarList(1, v);
     }
 }
 
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_
     );