diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.C b/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.C
index fe96d0bfe7b8db5ad8bd93301f4fa891f6c69988..45050a171b67130246a28f074ba591103392087e 100644
--- a/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.C
@@ -70,7 +70,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiGModels::basicSubGrid::G() const
     const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
 
     tmp<volScalarField> tGtot = XiGModel_->G();
-    volScalarField& Gtot = tGtot();
+    volScalarField& Gtot = tGtot.ref();
 
     const scalarField Cw = pow(Su_.mesh().V(), 2.0/3.0);
     scalarField N(Nv.internalField()*Cw);
diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.C b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.C
index 7454b2f7eecb7e340ade262d2849737143db20d4..0de82a82e4c2498b75c26e215116e7d165731570 100644
--- a/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.C
@@ -113,7 +113,7 @@ Foam::tmp<Foam::volSymmTensorField> Foam::PDRDragModels::basic::Dcu() const
         )
     );
 
-    volSymmTensorField& DragDcu = tDragDcu();
+    volSymmTensorField& DragDcu = tDragDcu.ref();
 
     if (on_)
     {
@@ -147,7 +147,7 @@ Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basic::Gk() const
         )
     );
 
-    volScalarField& Gk = tGk();
+    volScalarField& Gk = tGk.ref();
 
     if (on_)
     {
diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
index ebc92684a440794e4ab6b90d0ea0365d9d2a558a..afe97f7f7f38c406a67fd0d66c9c42e38735e783 100644
--- a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
@@ -161,9 +161,9 @@ void PDRkEpsilon::correct()
       - fvm::Sp(C2_*betav*rho_*epsilon_/k_, epsilon_)
     );
 
-    epsEqn().relax();
+    epsEqn.ref().relax();
 
-    epsEqn().boundaryManipulate(epsilon_.boundaryField());
+    epsEqn.ref().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
     bound(epsilon_, epsilonMin_);
@@ -182,7 +182,7 @@ void PDRkEpsilon::correct()
       - fvm::Sp(betav*rho_*epsilon_/k_, k_)
     );
 
-    kEqn().relax();
+    kEqn.ref().relax();
     solve(kEqn);
     bound(k_, kMin_);
 
diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
index ad8452f5a43795be9bb00501fa71de0165207b62..0208818ee0034cb202e0a4e1a5d213c8c2a17972 100644
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -108,7 +108,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
             dimensionedScalar("XiEq", dimless, 0.0)
         )
     );
-    volScalarField& xieq = tXiEq();
+    volScalarField& xieq = tXiEq.ref();
 
     forAll(xieq, celli)
     {
diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C
index 59f5c81bd340c4cedd9684a1e5adcb9d9fbc39a3..a84aca32108953ab41f906b96e4cd4926117018e 100644
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C
@@ -113,7 +113,7 @@ Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
             dimensionedScalar("zero", Nv.dimensions(), 0.0)
         )
     );
-    volScalarField& N = tN();
+    volScalarField& N = tN.ref();
     N.internalField() = Nv.internalField()*pow(mesh.V(), 2.0/3.0);
 
     volSymmTensorField ns
diff --git a/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C b/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
index d4d77822878c37ba04401c25b552cb9da77af6eb..fd07b0b3cbfcada4c87eda806251c708cdaf73b8 100644
--- a/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
+++ b/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -259,7 +259,7 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
         )
     );
 
-    volScalarField& Su0 = tSu0();
+    volScalarField& Su0 = tSu0.ref();
 
     forAll(Su0, celli)
     {
@@ -306,7 +306,7 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
         )
     );
 
-    volScalarField& Su0 = tSu0();
+    volScalarField& Su0 = tSu0.ref();
 
     forAll(Su0, celli)
     {
@@ -358,7 +358,7 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
         )
     );
 
-    volScalarField& ma = tMa();
+    volScalarField& ma = tMa.ref();
 
     forAll(ma, celli)
     {
diff --git a/applications/solvers/combustion/reactingFoam/UEqn.H b/applications/solvers/combustion/reactingFoam/UEqn.H
index 907c9934e2b27ff86d3b862fbaa1d3fc804ff5c4..e4caa1cea4acd908e2a31004bfe1cf67278bbff7 100644
--- a/applications/solvers/combustion/reactingFoam/UEqn.H
+++ b/applications/solvers/combustion/reactingFoam/UEqn.H
@@ -2,7 +2,7 @@
 
 MRF.correctBoundaryVelocity(U);
 
-tmp<fvVectorMatrix> UEqn
+tmp<fvVectorMatrix> tUEqn
 (
     fvm::ddt(rho, U) + fvm::div(phi, U)
   + MRF.DDt(rho, U)
@@ -10,14 +10,15 @@ tmp<fvVectorMatrix> UEqn
  ==
     fvOptions(rho, U)
 );
+fvVectorMatrix& UEqn = tUEqn.ref();
 
-UEqn().relax();
+UEqn.relax();
 
-fvOptions.constrain(UEqn());
+fvOptions.constrain(UEqn);
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index ae5478c5b6ea0e202fd7d9e1afc64f5a17ba9dfd..f5b71412247bef3a384e7e7f7bc0713268305c0a 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -3,13 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
diff --git a/applications/solvers/combustion/reactingFoam/pcEqn.H b/applications/solvers/combustion/reactingFoam/pcEqn.H
index 8bef07a5bbbb1c72e37c5e312784464748216611..bf8eb899aff9b436e5a72bb976d0868d1815f531 100644
--- a/applications/solvers/combustion/reactingFoam/pcEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pcEqn.H
@@ -3,13 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
-volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
-volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
+volScalarField rAU(1.0/UEqn.A());
+volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
index cd191534625765c8be13a231e046d79a07b3d999..3d02c63381e0d423fb4a9c036f4f4885a3b3f344 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
@@ -5,9 +5,9 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p;
 
-    volScalarField rAU(1.0/UEqn().A());
+    volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
     if (pimple.transonic())
     {
diff --git a/applications/solvers/combustion/reactingFoam/setRDeltaT.H b/applications/solvers/combustion/reactingFoam/setRDeltaT.H
index b84be5a66d9356296b446d99dd902e1b73b04cce..12bdceb7e521e0cd809f708697cfb3d2a86e8091 100644
--- a/applications/solvers/combustion/reactingFoam/setRDeltaT.H
+++ b/applications/solvers/combustion/reactingFoam/setRDeltaT.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 {
-    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rDeltaT = trDeltaT.ref();
 
     const dictionary& pimpleDict = pimple.dict();
 
diff --git a/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H b/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H
index 7d361ac9f84d8305524c8de91d2c23d13f765953..bc94840d1f357782b70046ac3cec0bf74f5b7324 100644
--- a/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H
+++ b/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H
@@ -22,7 +22,7 @@ tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
         )
     );
 
-    GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf();
+    GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf.ref();
 
     sf.rename(vf.name() + '_' + dir.name());
 
diff --git a/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H b/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H
index e469bafa346c78f4a94aadcba842626841204650..806ddd3832433c1a0eaca7995bee31f630a6b2d2 100644
--- a/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H
+++ b/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H
@@ -1,5 +1,5 @@
 {
-    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rDeltaT = trDeltaT.ref();
 
     scalar rDeltaTSmoothingCoeff
     (
diff --git a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
index 907c9934e2b27ff86d3b862fbaa1d3fc804ff5c4..e4caa1cea4acd908e2a31004bfe1cf67278bbff7 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/UEqn.H
@@ -2,7 +2,7 @@
 
 MRF.correctBoundaryVelocity(U);
 
-tmp<fvVectorMatrix> UEqn
+tmp<fvVectorMatrix> tUEqn
 (
     fvm::ddt(rho, U) + fvm::div(phi, U)
   + MRF.DDt(rho, U)
@@ -10,14 +10,15 @@ tmp<fvVectorMatrix> UEqn
  ==
     fvOptions(rho, U)
 );
+fvVectorMatrix& UEqn = tUEqn.ref();
 
-UEqn().relax();
+UEqn.relax();
 
-fvOptions.constrain(UEqn());
+fvOptions.constrain(UEqn);
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index ae5478c5b6ea0e202fd7d9e1afc64f5a17ba9dfd..f5b71412247bef3a384e7e7f7bc0713268305c0a 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -3,13 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index 8bef07a5bbbb1c72e37c5e312784464748216611..bf8eb899aff9b436e5a72bb976d0868d1815f531 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -3,13 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
-volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
-volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
+volScalarField rAU(1.0/UEqn.A());
+volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index 615682db12d608ce8a393c7bfd2ff86c73816535..4abc5cd4aadc7639da4f6ec347091db20109c2c2 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -3,13 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H b/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
index 30a3908b05c3c77b809ddc8341393b0d5de901c8..ba2db0d04a42fd01a1275f3b03135a2f7274c179 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/setRDeltaT.H
@@ -1,5 +1,5 @@
 {
-    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rDeltaT = trDeltaT.ref();
 
     const dictionary& pimpleDict = pimple.dict();
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
index 396b923494ceec1c3d112ca925ad55343ee9e503..a6f410250bd53abff192ee5aef7a0db525486d92 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(rho, U)
@@ -10,11 +10,12 @@
      ==
         fvOptions(rho, U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 4c6276711c195fb4e864a58fccc45999acc7b77c..b085a2b9d1c3077e214600758ca9fbbb8911d0ec 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -1,8 +1,8 @@
 {
-    volScalarField rAU(1.0/UEqn().A());
+    volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
-    UEqn.clear();
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+    tUEqn.clear();
 
     bool closedVolume = false;
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
index 7b77048a999b0bdfd322144f5450f8796b120b8a..51b3865165718a89eba453ddb781ee2bf5991bb1 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
@@ -1,7 +1,7 @@
-volScalarField rAU(1.0/UEqn().A());
-volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
-volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
-UEqn.clear();
+volScalarField rAU(1.0/UEqn.A());
+volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+tUEqn.clear();
 
 bool closedVolume = false;
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/UEqn.H
index fcb8f713d6aade67647af6085455e4b30e4cc702..3ec153c5679eb92661ea038284a0231edba525c8 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/UEqn.H
@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(rho, U)
@@ -10,8 +10,9 @@
      ==
         fvOptions(rho, U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
     // Include the porous media resistance and solve the momentum equation
     // either implicit in the tensorial resistance or transport using by
@@ -22,18 +23,18 @@
 
     if (pressureImplicitPorosity)
     {
-        tmp<volTensorField> tTU = tensor(I)*UEqn().A();
-        pZones.addResistance(UEqn(), tTU());
+        tmp<volTensorField> tTU = tensor(I)*UEqn.A();
+        pZones.addResistance(UEqn, tTU.ref());
         trTU = inv(tTU());
-        trTU().rename("rAU");
+        trTU.ref().rename("rAU");
 
-        fvOptions.constrain(UEqn());
+        fvOptions.constrain(UEqn);
 
         volVectorField gradp(fvc::grad(p));
 
         for (int UCorr=0; UCorr<nUCorr; UCorr++)
         {
-            U = trTU() & (UEqn().H() - gradp);
+            U = trTU() & (UEqn.H() - gradp);
         }
         U.correctBoundaryConditions();
 
@@ -41,14 +42,14 @@
     }
     else
     {
-        pZones.addResistance(UEqn());
+        pZones.addResistance(UEqn);
 
-        fvOptions.constrain(UEqn());
+        fvOptions.constrain(UEqn);
 
-        solve(UEqn() == -fvc::grad(p));
+        solve(UEqn == -fvc::grad(p));
 
         fvOptions.correct(U);
 
-        trAU = 1.0/UEqn().A();
-        trAU().rename("rAU");
+        trAU = 1.0/UEqn.A();
+        trAU.ref().rename("rAU");
     }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
index d54f7eec248587ece2b20ae169935196bd9689cc..b586cd15c0d0715449b2effff47c58107483bc04 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
@@ -4,15 +4,15 @@
     tmp<volVectorField> tHbyA;
     if (pressureImplicitPorosity)
     {
-        tHbyA = constrainHbyA(trTU()&UEqn().H(), U, p);
+        tHbyA = constrainHbyA(trTU()&UEqn.H(), U, p);
     }
     else
     {
-        tHbyA = constrainHbyA(trAU()*UEqn().H(), U, p);
+        tHbyA = constrainHbyA(trAU()*UEqn.H(), U, p);
     }
-    volVectorField& HbyA = tHbyA();
+    volVectorField& HbyA = tHbyA.ref();
 
-    UEqn.clear();
+    tUEqn.clear();
 
     bool closedVolume = false;
 
@@ -51,13 +51,15 @@
             );
         }
 
-        tpEqn().setReference(pRefCell, pRefValue);
+        fvScalarMatrix& pEqn = tpEqn.ref();
 
-        tpEqn().solve();
+        pEqn.setReference(pRefCell, pRefValue);
+
+        pEqn.solve();
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi = phiHbyA - tpEqn().flux();
+            phi = phiHbyA - pEqn.flux();
         }
     }
 
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
index 72bf78c7fb113ed244743a5faf7d01813a21dd66..d8a347f87d73f3524052bc597bdb416fb0b50cf9 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(U)
@@ -10,16 +10,17 @@
      ==
         fvOptions(U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
     if (simple.momentumPredictor())
     {
         solve
         (
-            UEqn()
+            UEqn
           ==
             fvc::reconstruct
             (
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
index 13200035d85a48dccd9882d832dada6a1d898c05..7df0ac46d290683e2a53a8d085bc1dad26fad5c0 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
@@ -1,9 +1,9 @@
 {
-    volScalarField rAU("rAU", 1.0/UEqn().A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p_rgh));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
 
-    UEqn.clear();
+    tUEqn.clear();
 
     surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
index 04796b475430a3ed53bbb72a47de0a01e66af4d6..1a47d06a8e8255150a6fc9d6887a1a7fbc833900 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(rho, U)
@@ -10,16 +10,17 @@
      ==
         fvOptions(rho, U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
     if (simple.momentumPredictor())
     {
         solve
         (
-            UEqn()
+            UEqn
          ==
             fvc::reconstruct
             (
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index 132e440bcd25146947fa38fb0cab2480457f479b..a744b9d19fbc52fc502dba3ac275c6d27a8ab7e4 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -2,10 +2,10 @@
     rho = thermo.rho();
     rho.relax();
 
-    volScalarField rAU("rAU", 1.0/UEqn().A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p_rgh));
-    UEqn.clear();
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+    tUEqn.clear();
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H
index b65ef6bfdb2831ec71b687021dc754a20b77107b..ed0402c2d0287d24fa0395ac801038899b43cc3c 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H
@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(rho, U)
@@ -10,14 +10,15 @@
      ==
         fvOptions(rho, U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
     solve
     (
-        UEqn()
+        UEqn
      ==
         fvc::reconstruct
         (
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index 341de6b13e9bc9bc3e0e37273007838b2a304d38..3fdadc042a2684935f04b5563e8a87a54181d352 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -4,10 +4,10 @@
     rho = min(rho, rhoMax[i]);
     rho.relax();
 
-    volScalarField rAU("rAU", 1.0/UEqn().A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p_rgh));
-    UEqn.clear();
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+    tUEqn.clear();
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
index 563ad65eba2e0d9c0083de8d4cf79831d3804c5e..a12d0e8146703889de6ac56ff3963b44b32d3d21 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::ddt(rho, U) + fvm::div(phi, U)
       + MRF.DDt(rho, U)
@@ -10,16 +10,17 @@
      ==
         fvOptions(rho, U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
     if (momentumPredictor)
     {
         solve
         (
-            UEqn()
+            UEqn
           ==
             fvc::reconstruct
             (
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 0b82c4406860b4355d971445a7c6df4835027301..358535efe97f9ed21a6dab5ba5cf6675f80c63f6 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -5,9 +5,9 @@
 
     rho = thermo.rho();
 
-    volScalarField rAU("rAU", 1.0/UEqn().A());
+    volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p_rgh));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
 
     surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
index 49fd39e45b78e2a7ad389897076677352afdc1ee..30a670768fff800ad521fc9e03c4d990584691d3 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
@@ -6,7 +6,7 @@ if (finalIter)
 {
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        tmp<fvScalarMatrix> hEqn
+        fvScalarMatrix hEqn
         (
             fvm::ddt(betav*rho, h)
           - (
@@ -18,11 +18,11 @@ if (finalIter)
             fvOptions(rho, h)
         );
 
-        hEqn().relax();
+        hEqn.relax();
 
-        fvOptions.constrain(hEqn());
+        fvOptions.constrain(hEqn);
 
-        hEqn().solve(mesh.solver(h.select(finalIter)));
+        hEqn.solve(mesh.solver(h.select(finalIter)));
 
         fvOptions.correct(h);
     }
diff --git a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
index fbd8e5f15077722b6d2be224a37e3d47b26af0a8..cdc1616e2b4ebbb4afdb581aefd4646891e64adf 100644
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
@@ -107,7 +107,7 @@ int main(int argc, char *argv[])
         {
             // Momentum predictor
 
-            tmp<fvVectorMatrix> UEqn
+            tmp<fvVectorMatrix> tUEqn
             (
                 fvm::div(phi, U)
               + turbulence->divDevReff(U)
@@ -115,18 +115,19 @@ int main(int argc, char *argv[])
              ==
                 fvOptions(U)
             );
+            fvVectorMatrix& UEqn = tUEqn.ref();
 
-            UEqn().relax();
+            UEqn.relax();
 
-            fvOptions.constrain(UEqn());
+            fvOptions.constrain(UEqn);
 
-            solve(UEqn() == -fvc::grad(p));
+            solve(UEqn == -fvc::grad(p));
 
             fvOptions.correct(U);
 
-            volScalarField rAU(1.0/UEqn().A());
-            volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
-            UEqn.clear();
+            volScalarField rAU(1.0/UEqn.A());
+            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+            tUEqn.clear();
             surfaceScalarField phiHbyA
             (
                 "phiHbyA",
@@ -180,7 +181,7 @@ int main(int argc, char *argv[])
 
             zeroCells(adjointTransposeConvection, inletCells);
 
-            tmp<fvVectorMatrix> UaEqn
+            tmp<fvVectorMatrix> tUaEqn
             (
                 fvm::div(-phi, Ua)
               - adjointTransposeConvection
@@ -189,19 +190,20 @@ int main(int argc, char *argv[])
              ==
                 fvOptions(Ua)
             );
+            fvVectorMatrix& UaEqn = tUaEqn.ref();
 
-            UaEqn().relax();
+            UaEqn.relax();
 
-            fvOptions.constrain(UaEqn());
+            fvOptions.constrain(UaEqn);
 
-            solve(UaEqn() == -fvc::grad(pa));
+            solve(UaEqn == -fvc::grad(pa));
 
             fvOptions.correct(Ua);
 
-            volScalarField rAUa(1.0/UaEqn().A());
+            volScalarField rAUa(1.0/UaEqn.A());
             volVectorField HbyAa("HbyAa", Ua);
-            HbyAa = rAUa*UaEqn().H();
-            UaEqn.clear();
+            HbyAa = rAUa*UaEqn.H();
+            tUaEqn.clear();
             surfaceScalarField phiHbyAa
             (
                 "phiHbyAa",
diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/UrelEqn.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/UrelEqn.H
index 5aa6cef0b147a6126a32a94c420f4bceb88d9ccd..0d55edd5a0892235ac482e9cfc360968d92d6873 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/UrelEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/UrelEqn.H
@@ -1,5 +1,5 @@
     // Relative momentum predictor
-    tmp<fvVectorMatrix> UrelEqn
+    tmp<fvVectorMatrix> tUrelEqn
     (
         fvm::ddt(Urel)
       + fvm::div(phi, Urel)
@@ -8,11 +8,12 @@
      ==
         fvOptions(Urel)
     );
+    fvVectorMatrix& UrelEqn = tUrelEqn.ref();
 
-    UrelEqn().relax();
+    UrelEqn.relax();
 
-    fvOptions.constrain(UrelEqn());
+    fvOptions.constrain(UrelEqn);
 
-    solve(UrelEqn() == -fvc::grad(p));
+    solve(UrelEqn == -fvc::grad(p));
 
     fvOptions.correct(Urel);
diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
index 0e3936bc5d77ddb4c7a0ce55c063f3054bdf5017..aaeb12ae0671f60aa9f28c1b7347b666252ce873 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
@@ -1,6 +1,6 @@
-volScalarField rAUrel(1.0/UrelEqn().A());
+volScalarField rAUrel(1.0/UrelEqn.A());
 volVectorField HbyA("HbyA", Urel);
-HbyA = rAUrel*UrelEqn().H();
+HbyA = rAUrel*UrelEqn.H();
 
 surfaceScalarField phiHbyA
 (
@@ -15,7 +15,7 @@ tmp<volScalarField> rAtUrel(rAUrel);
 
 if (pimple.consistent())
 {
-    rAtUrel = 1.0/max(1.0/rAUrel - UrelEqn().H1(), 0.1/rAUrel);
+    rAtUrel = 1.0/max(1.0/rAUrel - UrelEqn.H1(), 0.1/rAUrel);
     phiHbyA +=
         fvc::interpolate(rAtUrel() - rAUrel)*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAUrel - rAtUrel())*fvc::grad(p);
@@ -23,7 +23,7 @@ if (pimple.consistent())
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UrelEqn.clear();
+    tUrelEqn.clear();
 }
 
 // Update the pressure BCs to ensure flux consistency
diff --git a/applications/solvers/incompressible/pimpleFoam/UEqn.H b/applications/solvers/incompressible/pimpleFoam/UEqn.H
index 880c0066e363ad1dee9c52b4b86cf1d844fad65a..1a6b404885a6f8648ae02cae50e3f57f7b4d6c99 100644
--- a/applications/solvers/incompressible/pimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/UEqn.H
@@ -2,7 +2,7 @@
 
 MRF.correctBoundaryVelocity(U);
 
-tmp<fvVectorMatrix> UEqn
+tmp<fvVectorMatrix> tUEqn
 (
     fvm::ddt(U) + fvm::div(phi, U)
   + MRF.DDt(U)
@@ -10,14 +10,15 @@ tmp<fvVectorMatrix> UEqn
  ==
     fvOptions(U)
 );
+fvVectorMatrix& UEqn = tUEqn.ref();
 
-UEqn().relax();
+UEqn.relax();
 
-fvOptions.constrain(UEqn());
+fvOptions.constrain(UEqn);
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
 }
diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H
index 6d4eb9033471c91752c86defe72170512d8dc854..ffce1373263fc703f13b32a16e742c5595a22830 100644
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@@ -1,5 +1,5 @@
-volScalarField rAU(1.0/UEqn().A());
-volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
+volScalarField rAU(1.0/UEqn.A());
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
@@ -15,7 +15,7 @@ tmp<volScalarField> rAtU(rAU);
 
 if (pimple.consistent())
 {
-    rAtU = 1.0/max(1.0/rAU - UEqn().H1(), 0.1/rAU);
+    rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
     phiHbyA +=
         fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU())*fvc::grad(p);
@@ -23,7 +23,7 @@ if (pimple.consistent())
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 // Update the pressure BCs to ensure flux consistency
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index f4106997c607297612909e73090a6e24cd57ca29..452f13fb1f126e7af325f9bad7f17a7d40485a4e 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
@@ -1,5 +1,5 @@
-volScalarField rAU(1.0/UEqn().A());
-volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
+volScalarField rAU(1.0/UEqn.A());
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
@@ -20,7 +20,7 @@ tmp<volScalarField> rAtU(rAU);
 
 if (pimple.consistent())
 {
-    rAtU = 1.0/max(1.0/rAU - UEqn().H1(), 0.1/rAU);
+    rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
     phiHbyA +=
         fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU())*fvc::grad(p);
@@ -28,7 +28,7 @@ if (pimple.consistent())
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 // Update the pressure BCs to ensure flux consistency
diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H
index 75df41d00b62a2a023c2972e014dfba5b78e355f..6a2070eb4e4102197c044d34f67af7e620f77834 100644
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/UrelEqn.H
@@ -1,6 +1,6 @@
     // Relative momentum predictor
 
-    tmp<fvVectorMatrix> UrelEqn
+    tmp<fvVectorMatrix> tUrelEqn
     (
         fvm::div(phi, Urel)
       + turbulence->divDevReff(Urel)
@@ -8,11 +8,12 @@
      ==
         fvOptions(Urel)
     );
+    fvVectorMatrix& UrelEqn = tUrelEqn.ref();
 
-    UrelEqn().relax();
+    UrelEqn.relax();
 
-    fvOptions.constrain(UrelEqn());
+    fvOptions.constrain(UrelEqn);
 
-    solve(UrelEqn() == -fvc::grad(p));
+    solve(UrelEqn == -fvc::grad(p));
 
     fvOptions.correct(Urel);
diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
index 636a4c571de381fc5c8cff3d69ffe929ea54d764..9982ddf546798fbb979a59af5cc7efd99a6477cf 100644
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
@@ -1,7 +1,7 @@
 {
-    volScalarField rAUrel(1.0/UrelEqn().A());
+    volScalarField rAUrel(1.0/UrelEqn.A());
     volVectorField HbyA("HbyA", Urel);
-    HbyA = rAUrel*UrelEqn().H();
+    HbyA = rAUrel*UrelEqn.H();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
     adjustPhi(phiHbyA, Urel, p);
@@ -10,13 +10,13 @@
 
     if (simple.consistent())
     {
-        rAtUrel = 1.0/(1.0/rAUrel - UrelEqn().H1());
+        rAtUrel = 1.0/(1.0/rAUrel - UrelEqn.H1());
         phiHbyA +=
             fvc::interpolate(rAtUrel() - rAUrel)*fvc::snGrad(p)*mesh.magSf();
         HbyA -= (rAUrel - rAtUrel())*fvc::grad(p);
     }
 
-    UrelEqn.clear();
+    tUrelEqn.clear();
 
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, Urel, phiHbyA, rAtUrel());
diff --git a/applications/solvers/incompressible/simpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/UEqn.H
index 12b5621636669765be3e4919fdaab45771e94b11..31c081792c1dd27fa0b0854fbbcc5f661ee25602 100644
--- a/applications/solvers/incompressible/simpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/UEqn.H
@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(U)
@@ -10,11 +10,12 @@
      ==
         fvOptions(U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H
index b506bf43cfeb94fc53e5722d610b28c9f382a82c..5d5dde848c7d83355d24659f3e21ef9c5e55eba3 100644
--- a/applications/solvers/incompressible/simpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
-    volScalarField rAU(1.0/UEqn().A());
-    volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
+    volScalarField rAU(1.0/UEqn.A());
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
     MRF.makeRelative(phiHbyA);
     adjustPhi(phiHbyA, U, p);
@@ -9,13 +9,13 @@
 
     if (simple.consistent())
     {
-        rAtU = 1.0/(1.0/rAU - UEqn().H1());
+        rAtU = 1.0/(1.0/rAU - UEqn.H1());
         phiHbyA +=
             fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
         HbyA -= (rAU - rAtU())*fvc::grad(p);
     }
 
-    UEqn.clear();
+    tUEqn.clear();
 
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, U, phiHbyA, rAtU(), MRF);
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H
index 035b83ebd1bb2d416c82c1287db90b1a669f0f5b..e2c9917d5ae19cbd85ed692c3620b65b643e72fe 100644
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/UEqn.H
@@ -2,7 +2,7 @@
 
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(U)
@@ -10,8 +10,9 @@
       ==
         fvOptions(U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
     // Include the porous media resistance and solve the momentum equation
     // either implicit in the tensorial resistance or transport using by
@@ -22,18 +23,18 @@
 
     if (pressureImplicitPorosity)
     {
-        tmp<volTensorField> tTU = tensor(I)*UEqn().A();
-        pZones.addResistance(UEqn(), tTU());
+        tmp<volTensorField> tTU = tensor(I)*UEqn.A();
+        pZones.addResistance(UEqn, tTU.ref());
         trTU = inv(tTU());
-        trTU().rename("rAU");
+        trTU.ref().rename("rAU");
 
-        fvOptions.constrain(UEqn());
+        fvOptions.constrain(UEqn);
 
         volVectorField gradp(fvc::grad(p));
 
         for (int UCorr=0; UCorr<nUCorr; UCorr++)
         {
-            U = trTU() & (UEqn().H() - gradp);
+            U = trTU() & (UEqn.H() - gradp);
         }
         U.correctBoundaryConditions();
 
@@ -41,14 +42,14 @@
     }
     else
     {
-        pZones.addResistance(UEqn());
+        pZones.addResistance(UEqn);
 
-        fvOptions.constrain(UEqn());
+        fvOptions.constrain(UEqn);
 
-        solve(UEqn() == -fvc::grad(p));
+        solve(UEqn == -fvc::grad(p));
 
         fvOptions.correct(U);
 
-        trAU = 1.0/UEqn().A();
-        trAU().rename("rAU");
+        trAU = 1.0/UEqn.A();
+        trAU.ref().rename("rAU");
     }
diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
index e4ad23aabc2c7ce80b5117025519e3c0bba979a3..68bdcc77e0ee43efa5b914150f70b0fd47e6cbfd 100644
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
@@ -1,15 +1,15 @@
 tmp<volVectorField> tHbyA;
 if (pressureImplicitPorosity)
 {
-    tHbyA = constrainHbyA(trTU()&UEqn().H(), U, p);
+    tHbyA = constrainHbyA(trTU()&UEqn.H(), U, p);
 }
 else
 {
-    tHbyA = constrainHbyA(trAU()*UEqn().H(), U, p);
+    tHbyA = constrainHbyA(trAU()*UEqn.H(), U, p);
 }
-volVectorField& HbyA = tHbyA();
+volVectorField& HbyA = tHbyA.ref();
 
-UEqn.clear();
+tUEqn.clear();
 surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
 
 MRF.makeRelative(phiHbyA);
@@ -29,13 +29,15 @@ while (simple.correctNonOrthogonal())
         tpEqn = (fvm::laplacian(trAU(), p) == fvc::div(phiHbyA));
     }
 
-    tpEqn().setReference(pRefCell, pRefValue);
+    fvScalarMatrix pEqn = tpEqn.ref();
 
-    tpEqn().solve();
+    pEqn.setReference(pRefCell, pRefValue);
+
+    pEqn.solve();
 
     if (simple.finalNonOrthogonalIter())
     {
-        phi = phiHbyA - tpEqn().flux();
+        phi = phiHbyA - pEqn.flux();
     }
 }
 
diff --git a/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options
index bd9cbcb6133ea6dd7c9b41509c231bdd930b8b8b..902afe59e4f2e3d14c15180e41b3820174b4a584 100644
--- a/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels \
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H b/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H
index 1c7d7cb5e293faa2a4171a694a2d6f27eeda9289..65b0388e12dcba4447824a85f2dbf740fd1086f8 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 {
-    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rDeltaT = trDeltaT.ref();
 
     const dictionary& pimpleDict = pimple.dict();
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H b/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H
index d599d491f83d784219429597018e59a96ffef74d..f51eeb49bf3ba7f3fa78028e895226faec408b8c 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 {
-    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rDeltaT = trDeltaT.ref();
 
     const dictionary& pimpleDict = pimple.dict();
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H
index 2f5a74d45b017c572361d9e6cfac8801802e4934..86cf5362f81ac1ee4c31eec659ff2a7481affe43 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/UEqn.H
@@ -1,6 +1,6 @@
     MRF.correctBoundaryVelocity(U);
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(rho, U)
@@ -10,11 +10,12 @@
       + parcels.SU(U)
       + fvOptions(rho, U)
     );
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
-    UEqn().relax();
+    UEqn.relax();
 
-    fvOptions.constrain(UEqn());
+    fvOptions.constrain(UEqn);
 
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
index aa775bc0d0ff0b2b73cbef8e9fffee1b84c3dc05..9b48bb939df2b32b8cfc2ab770cade94b13dadd0 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
@@ -3,10 +3,10 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p;
 
-    volScalarField rAU(1.0/UEqn().A());
+    volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
-    UEqn.clear();
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+    tUEqn.clear();
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
diff --git a/applications/solvers/lagrangian/sprayFoam/UEqn.H b/applications/solvers/lagrangian/sprayFoam/UEqn.H
index c26bf1cf65e887ed40b326784b9c29da05a45893..97a83f7d5b6fc5309e8e0bb9e71a72fdf31622ae 100644
--- a/applications/solvers/lagrangian/sprayFoam/UEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/UEqn.H
@@ -2,7 +2,7 @@
 
 MRF.correctBoundaryVelocity(U);
 
-tmp<fvVectorMatrix> UEqn
+tmp<fvVectorMatrix> tUEqn
 (
     fvm::ddt(rho, U) + fvm::div(phi, U)
   + MRF.DDt(rho, U)
@@ -12,14 +12,15 @@ tmp<fvVectorMatrix> UEqn
   + parcels.SU(U)
   + fvOptions(rho, U)
 );
+fvVectorMatrix& UEqn = tUEqn.ref();
 
-UEqn().relax();
+UEqn.relax();
 
-fvOptions.constrain(UEqn());
+fvOptions.constrain(UEqn);
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p));
+    solve(UEqn == -fvc::grad(p));
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);
diff --git a/applications/solvers/lagrangian/sprayFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/pEqn.H
index 5d28412683f5f9b9ec797aa0a3021c4a4fd021b0..1447f5b6be5b71db0523889e4349295c8aac40a7 100644
--- a/applications/solvers/lagrangian/sprayFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/pEqn.H
@@ -3,13 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
index 406128feed655457c3cdbfc9c44a9ad1bb8dd160..8522e609185a8673b0c96c54e5de1ce103fae54b 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
@@ -3,13 +3,13 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 if (pimple.transonic())
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index e1ac1288101f2aa02c764290790a2387a23b6ad1..c326215512ee705e7c82819d9879d6da9da3814e 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -40,8 +40,8 @@
                 psi1*fvm::ddt(p_rgh)
               + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
             );
-        deleteDemandDrivenData(p_rghEqnComp1().faceFluxCorrectionPtr());
-        p_rghEqnComp1().relax();
+        deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr());
+        p_rghEqnComp1.ref().relax();
 
         p_rghEqnComp2 =
             fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
@@ -50,8 +50,8 @@
                 psi2*fvm::ddt(p_rgh)
               + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
             );
-        deleteDemandDrivenData(p_rghEqnComp2().faceFluxCorrectionPtr());
-        p_rghEqnComp2().relax();
+        deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr());
+        p_rghEqnComp2.ref().relax();
     }
     else
     {
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index fb23042a87c50d5f249d81e9e25cc052ffb05677..a0a91c3255451b09dcfc9b0b8408a144760759cd 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -37,8 +37,8 @@
                 psi1*fvm::ddt(p_rgh)
               + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
             );
-        deleteDemandDrivenData(p_rghEqnComp1().faceFluxCorrectionPtr());
-        p_rghEqnComp1().relax();
+        deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr());
+        p_rghEqnComp1.ref().relax();
 
         p_rghEqnComp2 =
             fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
@@ -47,8 +47,8 @@
                 psi2*fvm::ddt(p_rgh)
               + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
             );
-        deleteDemandDrivenData(p_rghEqnComp2().faceFluxCorrectionPtr());
-        p_rghEqnComp2().relax();
+        deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr());
+        p_rghEqnComp2.ref().relax();
     }
     else
     {
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
index 5730273383fb92b2efabb3b0a3a40f4e7f114c56..14bbbe8708650e3aaab85a355e759832387ce386 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
@@ -190,7 +190,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::he
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        the() += phasei()*phasei().thermo().he(p, T);
+        the.ref() += phasei()*phasei().thermo().he(p, T);
     }
 
     return the;
@@ -213,7 +213,8 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::he
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        the() += scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells);
+        the.ref() +=
+            scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells);
     }
 
     return the;
@@ -236,7 +237,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::he
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        the() +=
+        the.ref() +=
             phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi);
     }
 
@@ -252,7 +253,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::hc() const
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        thc() += phasei()*phasei().thermo().hc();
+        thc.ref() += phasei()*phasei().thermo().hc();
     }
 
     return thc;
@@ -293,7 +294,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::rho() const
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        trho() += phasei()*phasei().thermo().rho();
+        trho.ref() += phasei()*phasei().thermo().rho();
     }
 
     return trho;
@@ -314,7 +315,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::rho
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        trho() +=
+        trho.ref() +=
             phasei().boundaryField()[patchi]*phasei().thermo().rho(patchi);
     }
 
@@ -330,7 +331,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cp() const
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tCp() += phasei()*phasei().thermo().Cp();
+        tCp.ref() += phasei()*phasei().thermo().Cp();
     }
 
     return tCp;
@@ -353,7 +354,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cp
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tCp() +=
+        tCp.ref() +=
             phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi);
     }
 
@@ -369,7 +370,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cv() const
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tCv() += phasei()*phasei().thermo().Cv();
+        tCv.ref() += phasei()*phasei().thermo().Cv();
     }
 
     return tCv;
@@ -392,7 +393,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cv
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tCv() +=
+        tCv.ref() +=
             phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi);
     }
 
@@ -408,7 +409,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::gamma() const
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tgamma() += phasei()*phasei().thermo().gamma();
+        tgamma.ref() += phasei()*phasei().thermo().gamma();
     }
 
     return tgamma;
@@ -431,7 +432,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::gamma
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tgamma() +=
+        tgamma.ref() +=
             phasei().boundaryField()[patchi]
            *phasei().thermo().gamma(p, T, patchi);
     }
@@ -448,7 +449,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cpv() const
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tCpv() += phasei()*phasei().thermo().Cpv();
+        tCpv.ref() += phasei()*phasei().thermo().Cpv();
     }
 
     return tCpv;
@@ -471,7 +472,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cpv
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tCpv() +=
+        tCpv.ref() +=
             phasei().boundaryField()[patchi]
            *phasei().thermo().Cpv(p, T, patchi);
     }
@@ -488,7 +489,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::CpByCpv() const
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tCpByCpv() += phasei()*phasei().thermo().CpByCpv();
+        tCpByCpv.ref() += phasei()*phasei().thermo().CpByCpv();
     }
 
     return tCpByCpv;
@@ -511,7 +512,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::CpByCpv
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tCpByCpv() +=
+        tCpByCpv.ref() +=
             phasei().boundaryField()[patchi]
            *phasei().thermo().CpByCpv(p, T, patchi);
     }
@@ -543,7 +544,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::kappa() const
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tkappa() += phasei()*phasei().thermo().kappa();
+        tkappa.ref() += phasei()*phasei().thermo().kappa();
     }
 
     return tkappa;
@@ -564,7 +565,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::kappa
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tkappa() +=
+        tkappa.ref() +=
             phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi);
     }
 
@@ -583,7 +584,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::kappaEff
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tkappaEff() += phasei()*phasei().thermo().kappaEff(alphat);
+        tkappaEff.ref() += phasei()*phasei().thermo().kappaEff(alphat);
     }
 
     return tkappaEff;
@@ -606,7 +607,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::kappaEff
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        tkappaEff() +=
+        tkappaEff.ref() +=
             phasei().boundaryField()[patchi]
            *phasei().thermo().kappaEff(alphat, patchi);
     }
@@ -626,7 +627,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::alphaEff
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        talphaEff() += phasei()*phasei().thermo().alphaEff(alphat);
+        talphaEff.ref() += phasei()*phasei().thermo().alphaEff(alphat);
     }
 
     return talphaEff;
@@ -649,7 +650,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::alphaEff
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        talphaEff() +=
+        talphaEff.ref() +=
             phasei().boundaryField()[patchi]
            *phasei().thermo().alphaEff(alphat, patchi);
     }
@@ -666,7 +667,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::rCv() const
 
     for (++phasei; phasei != phases_.end(); ++phasei)
     {
-        trCv() += phasei()/phasei().thermo().Cv();
+        trCv.ref() += phasei()/phasei().thermo().Cv();
     }
 
     return trCv;
@@ -696,7 +697,7 @@ Foam::multiphaseMixtureThermo::surfaceTensionForce() const
         )
     );
 
-    surfaceScalarField& stf = tstf();
+    surfaceScalarField& stf = tstf.ref();
 
     forAllConstIter(PtrDictionary<phaseModel>, phases_, phase1)
     {
@@ -922,7 +923,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::K
 {
     tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
 
-    correctContactAngle(alpha1, alpha2, tnHatfv().boundaryField());
+    correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryField());
 
     // Simple expression for curvature
     return -fvc::div(tnHatfv & mesh_.Sf());
@@ -949,7 +950,8 @@ Foam::multiphaseMixtureThermo::nearInterface() const
 
     forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
     {
-        tnearInt() = max(tnearInt(), pos(phase() - 0.01)*pos(0.99 - phase()));
+        tnearInt.ref() =
+            max(tnearInt(), pos(phase() - 0.01)*pos(0.99 - phase()));
     }
 
     return tnearInt;
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
index ccf4bc9794e598ddd30a1726b44290895fd1390c..26ffef825a5927cec9358ff1299df282e85e64d9 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
@@ -80,7 +80,7 @@
             }
             else
             {
-                p_rghEqnComp() += hmm;
+                p_rghEqnComp.ref() += hmm;
             }
 
             phasei++;
diff --git a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
index 68f027cd2463e750c63e09b01966799816634b27..3938fde30f0e5121081ddad3154f7414686461c2 100644
--- a/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
@@ -34,7 +34,7 @@
             (
                 alpha1,
                 alphaPhi,
-                talphaPhiCorr0(),
+                talphaPhiCorr0.ref(),
                 mixture.alphaMax(),
                 0
             );
@@ -73,7 +73,7 @@
             (
                 alpha1,
                 talphaPhiUn(),
-                talphaPhiCorr(),
+                talphaPhiCorr.ref(),
                 mixture.alphaMax(),
                 0
             );
diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H
index 055aac39669ca8993aa7683d5d12cf5c10a09639..9d6dd2dc18b68f8e2adbd1b422ad2a54c997d915 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqn.H
@@ -32,7 +32,8 @@
         }
 
         ocCoeff =
-            refCast<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha()).ocCoeff();
+            refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha())
+           .ocCoeff();
     }
     else
     {
@@ -104,7 +105,7 @@
         if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
         {
             Info<< "Applying the previous iteration compression flux" << endl;
-            MULES::correct(alpha1, alphaPhi, talphaPhiCorr0(), 1, 0);
+            MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 1, 0);
 
             alphaPhi += talphaPhiCorr0();
         }
@@ -150,7 +151,7 @@
             tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
             volScalarField alpha10("alpha10", alpha1);
 
-            MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr(), 1, 0);
+            MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr.ref(), 1, 0);
 
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
diff --git a/applications/solvers/multiphase/interFoam/setRDeltaT.H b/applications/solvers/multiphase/interFoam/setRDeltaT.H
index 01bb8aa411ce23424216641173a206e932ade3e1..dd0e6a4da45a946c9fbd7512b7c972da58d7043b 100644
--- a/applications/solvers/multiphase/interFoam/setRDeltaT.H
+++ b/applications/solvers/multiphase/interFoam/setRDeltaT.H
@@ -1,5 +1,5 @@
 {
-    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rDeltaT = trDeltaT.ref();
 
     const dictionary& pimpleDict = pimple.dict();
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
index 9cc90a9ebe4edd47d1c122e05def876089346c2e..150ccd1f787e98885e24a3c0a9fbec536737b776 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
@@ -62,7 +62,7 @@
 
         if (MULESCorr)
         {
-            talphaPhiCorr() -= talphaPhi();
+            talphaPhiCorr.ref() -= talphaPhi();
 
             volScalarField alpha100("alpha100", alpha10);
             alpha10 = alpha1;
@@ -72,7 +72,7 @@
                 geometricOneField(),
                 alpha1,
                 talphaPhi(),
-                talphaPhiCorr(),
+                talphaPhiCorr.ref(),
                 vDotvmcAlphal,
                 (
                     divU*(alpha10 - alpha100)
@@ -85,12 +85,12 @@
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
             {
-                talphaPhi() += talphaPhiCorr();
+                talphaPhi.ref() += talphaPhiCorr();
             }
             else
             {
                 alpha1 = 0.5*alpha1 + 0.5*alpha10;
-                talphaPhi() += 0.5*talphaPhiCorr();
+                talphaPhi.ref() += 0.5*talphaPhiCorr();
             }
         }
         else
@@ -100,7 +100,7 @@
                 geometricOneField(),
                 alpha1,
                 phi,
-                talphaPhiCorr(),
+                talphaPhiCorr.ref(),
                 vDotvmcAlphal,
                 (divU*alpha1 + vDotcAlphal)(),
                 1,
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 6755e96ecb31740748a439155ec4dd4dd2491115..27a5067edd74b322a2210aeb1aba996cfcdf9088 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -372,7 +372,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::K
 {
     tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
 
-    correctContactAngle(phase1, phase2, tnHatfv().boundaryField());
+    correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryField());
 
     // Simple expression for curvature
     return -fvc::div(tnHatfv & mesh_.Sf());
@@ -493,10 +493,11 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::rho() const
     PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
 
     tmp<volScalarField> trho = iter()*iter().rho();
+    volScalarField& rho = trho.ref();
 
     for (++iter; iter != phases_.end(); ++iter)
     {
-        trho() += iter()*iter().rho();
+        rho += iter()*iter().rho();
     }
 
     return trho;
@@ -509,10 +510,11 @@ Foam::multiphaseSystem::rho(const label patchi) const
     PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
 
     tmp<scalarField> trho = iter().boundaryField()[patchi]*iter().rho().value();
+    scalarField& rho = trho.ref();
 
     for (++iter; iter != phases_.end(); ++iter)
     {
-        trho() += iter().boundaryField()[patchi]*iter().rho().value();
+        rho += iter().boundaryField()[patchi]*iter().rho().value();
     }
 
     return trho;
@@ -524,10 +526,11 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::nu() const
     PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
 
     tmp<volScalarField> tmu = iter()*(iter().rho()*iter().nu());
+    volScalarField& mu = tmu.ref();
 
     for (++iter; iter != phases_.end(); ++iter)
     {
-        tmu() += iter()*(iter().rho()*iter().nu());
+        mu += iter()*(iter().rho()*iter().nu());
     }
 
     return tmu/rho();
@@ -542,10 +545,11 @@ Foam::multiphaseSystem::nu(const label patchi) const
     tmp<scalarField> tmu =
         iter().boundaryField()[patchi]
        *(iter().rho().value()*iter().nu().value());
+    scalarField& mu = tmu.ref();
 
     for (++iter; iter != phases_.end(); ++iter)
     {
-        tmu() +=
+        mu +=
             iter().boundaryField()[patchi]
            *(iter().rho().value()*iter().nu().value());
     }
@@ -592,7 +596,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::Cvm
 
             if (Cvm != Cvms_.end())
             {
-                tCvm() += Cvm()*phase2.rho()*phase2;
+                tCvm.ref() += Cvm()*phase2.rho()*phase2;
             }
             else
             {
@@ -600,7 +604,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::Cvm
 
                 if (Cvm != Cvms_.end())
                 {
-                    tCvm() += Cvm()*phase.rho()*phase2;
+                    tCvm.ref() += Cvm()*phase.rho()*phase2;
                 }
             }
         }
@@ -648,7 +652,7 @@ Foam::tmp<Foam::volVectorField> Foam::multiphaseSystem::Svm
 
             if (Cvm != Cvms_.end())
             {
-                tSvm() += Cvm()*phase2.rho()*phase2*phase2.DDtU();
+                tSvm.ref() += Cvm()*phase2.rho()*phase2*phase2.DDtU();
             }
             else
             {
@@ -656,7 +660,7 @@ Foam::tmp<Foam::volVectorField> Foam::multiphaseSystem::Svm
 
                 if (Cvm != Cvms_.end())
                 {
-                    tSvm() += Cvm()*phase.rho()*phase2*phase2.DDtU();
+                    tSvm.ref() += Cvm()*phase.rho()*phase2*phase2.DDtU();
                 }
             }
         }
@@ -673,7 +677,7 @@ Foam::tmp<Foam::volVectorField> Foam::multiphaseSystem::Svm
             )
         )
         {
-            tSvm().boundaryField()[patchi] = vector::zero;
+            tSvm.ref().boundaryField()[patchi] = vector::zero;
         }
     }
 
@@ -772,7 +776,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::dragCoeff
          || &phase == &dmIter()->phase2()
         )
         {
-            tdragCoeff() += *dcIter();
+            tdragCoeff.ref() += *dcIter();
         }
     }
 
@@ -818,7 +822,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
 
             if (sigma != sigmas_.end())
             {
-                tSurfaceTension() +=
+                tSurfaceTension.ref() +=
                     dimensionedScalar("sigma", dimSigma_, sigma())
                    *fvc::interpolate(K(phase1, phase2))*
                     (
@@ -853,7 +857,7 @@ Foam::multiphaseSystem::nearInterface() const
 
     forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
     {
-        tnearInt() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
+        tnearInt.ref() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
     }
 
     return tnearInt;
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index 27f0f51c35119d0c64846255e3fcd29ea965a76b..ea08ba29865801976b32ffb25801466fad2b786c 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -141,10 +141,11 @@ Foam::multiphaseMixture::rho() const
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
     tmp<volScalarField> trho = iter()*iter().rho();
+    volScalarField& rho = trho.ref();
 
     for (++iter; iter != phases_.end(); ++iter)
     {
-        trho() += iter()*iter().rho();
+        rho += iter()*iter().rho();
     }
 
     return trho;
@@ -157,10 +158,11 @@ Foam::multiphaseMixture::rho(const label patchi) const
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
     tmp<scalarField> trho = iter().boundaryField()[patchi]*iter().rho().value();
+    scalarField& rho = trho.ref();
 
     for (++iter; iter != phases_.end(); ++iter)
     {
-        trho() += iter().boundaryField()[patchi]*iter().rho().value();
+        rho += iter().boundaryField()[patchi]*iter().rho().value();
     }
 
     return trho;
@@ -173,10 +175,11 @@ Foam::multiphaseMixture::mu() const
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
     tmp<volScalarField> tmu = iter()*iter().rho()*iter().nu();
+    volScalarField& mu = tmu.ref();
 
     for (++iter; iter != phases_.end(); ++iter)
     {
-        tmu() += iter()*iter().rho()*iter().nu();
+        mu += iter()*iter().rho()*iter().nu();
     }
 
     return tmu;
@@ -192,10 +195,11 @@ Foam::multiphaseMixture::mu(const label patchi) const
         iter().boundaryField()[patchi]
        *iter().rho().value()
        *iter().nu(patchi);
+    scalarField& mu = tmu.ref();
 
     for (++iter; iter != phases_.end(); ++iter)
     {
-        tmu() +=
+        mu +=
             iter().boundaryField()[patchi]
            *iter().rho().value()
            *iter().nu(patchi);
@@ -212,10 +216,11 @@ Foam::multiphaseMixture::muf() const
 
     tmp<surfaceScalarField> tmuf =
         fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
+    surfaceScalarField& muf = tmuf.ref();
 
     for (++iter; iter != phases_.end(); ++iter)
     {
-        tmuf() +=
+        muf +=
             fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
     }
 
@@ -267,7 +272,7 @@ Foam::multiphaseMixture::surfaceTensionForce() const
         )
     );
 
-    surfaceScalarField& stf = tstf();
+    surfaceScalarField& stf = tstf.ref();
 
     forAllConstIter(PtrDictionary<phase>, phases_, iter1)
     {
@@ -518,7 +523,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
 {
     tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
 
-    correctContactAngle(alpha1, alpha2, tnHatfv().boundaryField());
+    correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryField());
 
     // Simple expression for curvature
     return -fvc::div(tnHatfv & mesh_.Sf());
@@ -545,7 +550,7 @@ Foam::multiphaseMixture::nearInterface() const
 
     forAllConstIter(PtrDictionary<phase>, phases_, iter)
     {
-        tnearInt() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
+        tnearInt.ref() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
     }
 
     return tnearInt;
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/UEqn.H b/applications/solvers/multiphase/potentialFreeSurfaceFoam/UEqn.H
index 794c0aca520f1df2e0c7a871125be70687bbbdae..421e45ae7da2d560c7e16341a78ed2b922fb21b5 100644
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/UEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/UEqn.H
@@ -1,6 +1,6 @@
 MRF.correctBoundaryVelocity(U);
 
-tmp<fvVectorMatrix> UEqn
+tmp<fvVectorMatrix> tUEqn
 (
     fvm::ddt(U) + fvm::div(phi, U)
   + MRF.DDt(U)
@@ -8,14 +8,15 @@ tmp<fvVectorMatrix> UEqn
  ==
     fvOptions(U)
 );
+fvVectorMatrix& UEqn = tUEqn.ref();
 
-UEqn().relax();
+UEqn.relax();
 
-fvOptions.constrain(UEqn());
+fvOptions.constrain(UEqn);
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p_gh));
+    solve(UEqn == -fvc::grad(p_gh));
 
     fvOptions.correct(U);
 }
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
index e182350114a40ef36d5208b0a7d79271f953e1ae..7d4c5f9039d24e6053b0a98181d94975f78e7458 100644
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
@@ -1,10 +1,10 @@
-volScalarField rAU(1.0/UEqn().A());
+volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p_gh));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_gh));
 
 if (pimple.nCorrPISO() <= 1)
 {
-    UEqn.clear();
+    tUEqn.clear();
 }
 
 surfaceScalarField phiHbyA
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H
index 04b187cf8509593f382f6d581bba3dc9cc87ddeb..a74f6c2664b39a9b971f7499dfb621773a1dca2e 100644
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/pEqn.H
@@ -1,11 +1,11 @@
 {
-    rAU = 1.0/UEqn().A();
+    rAU = 1.0/UEqn.A();
     surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn().H(), U, p_gh));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_gh));
 
     if (pimple.nCorrPISO() <= 1)
     {
-        UEqn.clear();
+        tUEqn.clear();
     }
 
     surfaceScalarField phiHbyA
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.C
index 5fd646c5ceeda5899f4107a19ae9d98538dc1329..ff5d3e5a6ed8e9b9fa762ed84dec9f7b7ab717e8 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -153,7 +153,7 @@ Foam::InterfaceCompositionModel<Thermo, OtherThermo>::D
         )
     );
 
-    volScalarField& D(tmpD());
+    volScalarField& D(tmpD.ref());
 
     forAll(p, cellI)
     {
@@ -207,7 +207,7 @@ Foam::InterfaceCompositionModel<Thermo, OtherThermo>::L
         )
     );
 
-    volScalarField& L(tmpL());
+    volScalarField& L(tmpL.ref());
 
     forAll(p, cellI)
     {
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C
index 886fefe8e89925d539a5030846f12f8f9b89fb0a..a2fbcd0bbcc61d9068ccfbd42e33fbd5c59e82af 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -111,7 +111,7 @@ Foam::saturationModels::polynomial::Tsat
         )
     );
 
-    volScalarField& Tsat = tTsat();
+    volScalarField& Tsat = tTsat.ref();
 
     forAll(Tsat,celli)
     {
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
index fb755c63471c55dbdf54197ccac15df79eb7c39d..c17b649212cb423493c743dfc893c9940b8103f9 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
@@ -46,7 +46,7 @@ Foam::tmp<Foam::volVectorField> Foam::wallLubricationModel::zeroGradWalls
     tmp<volVectorField> tFi
 ) const
 {
-    volVectorField& Fi = tFi();
+    volVectorField& Fi = tFi.ref();
     const fvPatchList& patches =  Fi.mesh().boundary();
 
     forAll(patches, patchi)
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
index 3cb567273f3cf5e20ce2335c41281819296b425b..78eb01cf1c3d1785a938a2c681ce2e395b5b5c80 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
@@ -199,15 +199,15 @@ Foam::BlendedInterfacialModel<ModelType>::K() const
 
     if (model_.valid())
     {
-        x() += model_->K()*(scalar(1) - f1() - f2());
+        x.ref() += model_->K()*(scalar(1) - f1() - f2());
     }
     if (model1In2_.valid())
     {
-        x() += model1In2_->K()*f1;
+        x.ref() += model1In2_->K()*f1;
     }
     if (model2In1_.valid())
     {
-        x() += model2In1_->K()*f2;
+        x.ref() += model2In1_->K()*f2;
     }
 
     if
@@ -216,7 +216,7 @@ Foam::BlendedInterfacialModel<ModelType>::K() const
      && (model_.valid() || model1In2_.valid() || model2In1_.valid())
     )
     {
-        correctFixedFluxBCs(x());
+        correctFixedFluxBCs(x.ref());
     }
 
     return x;
@@ -259,15 +259,15 @@ Foam::BlendedInterfacialModel<ModelType>::K(const scalar residualAlpha) const
 
     if (model_.valid())
     {
-        x() += model_->K(residualAlpha)*(scalar(1) - f1() - f2());
+        x.ref() += model_->K(residualAlpha)*(scalar(1) - f1() - f2());
     }
     if (model1In2_.valid())
     {
-        x() += model1In2_->K(residualAlpha)*f1;
+        x.ref() += model1In2_->K(residualAlpha)*f1;
     }
     if (model2In1_.valid())
     {
-        x() += model2In1_->K(residualAlpha)*f2;
+        x.ref() += model2In1_->K(residualAlpha)*f2;
     }
 
     if
@@ -276,7 +276,7 @@ Foam::BlendedInterfacialModel<ModelType>::K(const scalar residualAlpha) const
      && (model_.valid() || model1In2_.valid() || model2In1_.valid())
     )
     {
-        correctFixedFluxBCs(x());
+        correctFixedFluxBCs(x.ref());
     }
 
     return x;
@@ -325,17 +325,17 @@ Foam::BlendedInterfacialModel<ModelType>::Kf() const
 
     if (model_.valid())
     {
-        x() += model_->Kf()*(scalar(1) - f1() - f2());
+        x.ref() += model_->Kf()*(scalar(1) - f1() - f2());
     }
 
     if (model1In2_.valid())
     {
-        x() += model1In2_->Kf()*f1;
+        x.ref() += model1In2_->Kf()*f1;
     }
 
     if (model2In1_.valid())
     {
-        x() += model2In1_->Kf()*f2;
+        x.ref() += model2In1_->Kf()*f2;
     }
 
     if
@@ -344,7 +344,7 @@ Foam::BlendedInterfacialModel<ModelType>::Kf() const
      && (model_.valid() || model1In2_.valid() || model2In1_.valid())
     )
     {
-        correctFixedFluxBCs(x());
+        correctFixedFluxBCs(x.ref());
     }
 
     return x;
@@ -388,17 +388,17 @@ Foam::BlendedInterfacialModel<ModelType>::F() const
 
     if (model_.valid())
     {
-        x() += model_->F()*(scalar(1) - f1() - f2());
+        x.ref() += model_->F()*(scalar(1) - f1() - f2());
     }
 
     if (model1In2_.valid())
     {
-        x() += model1In2_->F()*f1;
+        x.ref() += model1In2_->F()*f1;
     }
 
     if (model2In1_.valid())
     {
-        x() -= model2In1_->F()*f2; // note : subtraction
+        x.ref() -= model2In1_->F()*f2; // note : subtraction
     }
 
     if
@@ -407,7 +407,7 @@ Foam::BlendedInterfacialModel<ModelType>::F() const
      && (model_.valid() || model1In2_.valid() || model2In1_.valid())
     )
     {
-        correctFixedFluxBCs(x());
+        correctFixedFluxBCs(x.ref());
     }
 
     return x;
@@ -456,17 +456,17 @@ Foam::BlendedInterfacialModel<ModelType>::Ff() const
 
     if (model_.valid())
     {
-        x() += model_->Ff()*(scalar(1) - f1() - f2());
+        x.ref() += model_->Ff()*(scalar(1) - f1() - f2());
     }
 
     if (model1In2_.valid())
     {
-        x() += model1In2_->Ff()*f1;
+        x.ref() += model1In2_->Ff()*f1;
     }
 
     if (model2In1_.valid())
     {
-        x() -= model2In1_->Ff()*f2; // note : subtraction
+        x.ref() -= model2In1_->Ff()*f2; // note : subtraction
     }
 
     if
@@ -475,7 +475,7 @@ Foam::BlendedInterfacialModel<ModelType>::Ff() const
      && (model_.valid() || model1In2_.valid() || model2In1_.valid())
     )
     {
-        correctFixedFluxBCs(x());
+        correctFixedFluxBCs(x.ref());
     }
 
     return x;
@@ -518,15 +518,15 @@ Foam::BlendedInterfacialModel<ModelType>::D() const
 
     if (model_.valid())
     {
-        x() += model_->D()*(scalar(1) - f1() - f2());
+        x.ref() += model_->D()*(scalar(1) - f1() - f2());
     }
     if (model1In2_.valid())
     {
-        x() += model1In2_->D()*f1;
+        x.ref() += model1In2_->D()*f1;
     }
     if (model2In1_.valid())
     {
-        x() += model2In1_->D()*f2;
+        x.ref() += model2In1_->D()*f2;
     }
 
     if
@@ -535,7 +535,7 @@ Foam::BlendedInterfacialModel<ModelType>::D() const
      && (model_.valid() || model1In2_.valid() || model2In1_.valid())
     )
     {
-        correctFixedFluxBCs(x());
+        correctFixedFluxBCs(x.ref());
     }
 
     return x;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
index 03c26c89dcbb3d6972d965b3efe295f798fd50d3..c63f6bfd447ba2b4d91e4d33cfb019ff9dc5438a 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
@@ -218,7 +218,7 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt
         {
             if (phase1 == &phase)
             {
-                tdmdt() += this->dmdt(pair);
+                tdmdt.ref() += this->dmdt(pair);
             }
 
             Swap(phase1, phase2);
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
index b737cd0322276dfe26acc9e9f026085cd38c0272..6ca75c6ed49838560e37e30cdb4daa8ff72e3688 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -112,13 +112,13 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
     // Add the appropriate pressure-work term
     if (he.name() == this->thermo_->phasePropertyName("e"))
     {
-        tEEqn() +=
+        tEEqn.ref() +=
             fvc::ddt(alpha)*this->thermo().p()
           + fvc::div(alphaPhi, this->thermo().p());
     }
     else if (this->thermo_->dpdt())
     {
-        tEEqn() -= alpha*this->fluid().dpdt();
+        tEEqn.ref() -= alpha*this->fluid().dpdt();
     }
 
     return tEEqn;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
index 4202907a2e1db20a63a74b02a8ca75621f647ccc..b9abe549a0bf1e7068d462f77fc8dce79df18f00 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,7 +57,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::phaseSystem::calcPhi
 
     for (label phasei=1; phasei<phaseModels.size(); phasei++)
     {
-        tmpPhi() +=
+        tmpPhi.ref() +=
             fvc::interpolate(phaseModels[phasei])*phaseModels[phasei].phi();
     }
 
@@ -197,7 +197,7 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rho() const
 
     for (label phasei=1; phasei<phaseModels_.size(); phasei++)
     {
-        tmpRho() += phaseModels_[phasei]*phaseModels_[phasei].rho();
+        tmpRho.ref() += phaseModels_[phasei]*phaseModels_[phasei].rho();
     }
 
     return tmpRho;
@@ -213,7 +213,7 @@ Foam::tmp<Foam::volVectorField> Foam::phaseSystem::U() const
 
     for (label phasei=1; phasei<phaseModels_.size(); phasei++)
     {
-        tmpU() += phaseModels_[phasei]*phaseModels_[phasei].U();
+        tmpU.ref() += phaseModels_[phasei]*phaseModels_[phasei].U();
     }
 
     return tmpU;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H
index 422f6324b1db3a15caeca4b1d7e8af961310751e..c1d98bcb4397ae50b6587e88f9bce133e31e80f4 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H
@@ -29,7 +29,7 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
             );
 
             EEqn->relax();
-            fvOptions.constrain(EEqn());
+            fvOptions.constrain(EEqn.ref());
             EEqn->solve();
         }
     }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseCompressibleTurbulenceModels/Make/options b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseCompressibleTurbulenceModels/Make/options
index ec5ca1596d4fbdb85cbc6c97e91694472e5b3299..5bdb4cea14078d765f89ceaefb7db8452f217dfd 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseCompressibleTurbulenceModels/Make/options
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseCompressibleTurbulenceModels/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I../multiphaseSystem/lnInclude \
     -I../../phaseSystems/lnInclude \
     -I../../interfacialModels/lnInclude\
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index a9de2e6d558f4128c73e24ede2c52399eac8f515..b81e174db878bfe6fbd061764179cc07e0c7ed19 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -477,7 +477,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::K
 {
     tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
 
-    correctContactAngle(phase1, phase2, tnHatfv().boundaryField());
+    correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryField());
 
     // Simple expression for curvature
     return -fvc::div(tnHatfv & mesh_.Sf());
@@ -568,7 +568,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
 
             if (cAlpha != cAlphas_.end())
             {
-                tSurfaceTension() +=
+                tSurfaceTension.ref() +=
                     fvc::interpolate(sigma(key12)*K(phase1, phase2))
                    *(
                         fvc::interpolate(phase2)*fvc::snGrad(phase1)
@@ -602,7 +602,7 @@ Foam::multiphaseSystem::nearInterface() const
 
     forAll(phases(), phasei)
     {
-        tnearInt() = max
+        tnearInt.ref() = max
         (
             tnearInt(),
             pos(phases()[phasei] - 0.01)*pos(0.99 - phases()[phasei])
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index 74e29d3df241e675165adc9230715bb593d44d07..42e0fe70e49adad0b8fab7715f7bce9a1977b1c7 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -457,7 +457,7 @@ while (pimple.correct())
     forAll(phases, phasei)
     {
         phaseModel& phase = phases[phasei];
-        phase.rho()() += phase.thermo().psi()*(p_rgh - p_rgh_0);
+        phase.thermo().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
     }
 
     // Correct p_rgh for consistency with p and the updated densities
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/setRDeltaT.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/setRDeltaT.H
index c0ca68835cee98c1ad49aaef2c6e24e9e8080418..fec017ea30462cb759075c6691600626fea1575d 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/setRDeltaT.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/setRDeltaT.H
@@ -1,5 +1,5 @@
 {
-    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rDeltaT = trDeltaT.ref();
 
     const dictionary& pimpleDict = pimple.dict();
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H
index 3060c444de311c1ca2c9ee9023e436ef897bd26f..0ac8e841297a47eec2873f24010e2e6641c20092 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H
@@ -23,7 +23,7 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
             );
 
             E1Eqn->relax();
-            fvOptions.constrain(E1Eqn());
+            fvOptions.constrain(E1Eqn.ref());
             E1Eqn->solve();
         }
     }
@@ -43,7 +43,7 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
             );
 
             E2eqn->relax();
-            fvOptions.constrain(E2eqn());
+            fvOptions.constrain(E2eqn.ref());
             E2eqn->solve();
         }
     }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
index d6ba189b129f00ded23082d0fb02739047cf433b..104bc0047ff7094bf8d5c83bbd36ca299ff76ea6 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -261,8 +261,8 @@ while (pimple.correct())
                     )
                 );
 
-            deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
-            pEqnComp1().relax();
+            deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
+            pEqnComp1.ref().relax();
         }
 
         if (phase2.compressible())
@@ -286,8 +286,8 @@ while (pimple.correct())
                       + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
                     )
                 );
-            deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
-            pEqnComp2().relax();
+            deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
+            pEqnComp2.ref().relax();
         }
     }
     else
@@ -317,7 +317,7 @@ while (pimple.correct())
     {
         if (pEqnComp1.valid())
         {
-            pEqnComp1() -= fluid.dmdt()/rho1;
+            pEqnComp1.ref() -= fluid.dmdt()/rho1;
         }
         else
         {
@@ -326,7 +326,7 @@ while (pimple.correct())
 
         if (pEqnComp2.valid())
         {
-            pEqnComp2() += fluid.dmdt()/rho2;
+            pEqnComp2.ref() += fluid.dmdt()/rho2;
         }
         else
         {
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index 57f491480e39845e31c03b008488814bf7fa3b7d..d618226411a105084d34428bb5d6628373c6249c 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -242,8 +242,8 @@ while (pimple.correct())
                   + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
                 )
             );
-            deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
-            pEqnComp1().relax();
+            deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
+            pEqnComp1.ref().relax();
         }
 
         if (phase2.compressible())
@@ -261,8 +261,8 @@ while (pimple.correct())
                   + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
                 )
             );
-            deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
-            pEqnComp2().relax();
+            deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
+            pEqnComp2.ref().relax();
         }
     }
     else
@@ -292,7 +292,7 @@ while (pimple.correct())
     {
         if (pEqnComp1.valid())
         {
-            pEqnComp1() -= fluid.dmdt()/rho1;
+            pEqnComp1.ref() -= fluid.dmdt()/rho1;
         }
         else
         {
@@ -301,7 +301,7 @@ while (pimple.correct())
 
         if (pEqnComp2.valid())
         {
-            pEqnComp2() += fluid.dmdt()/rho2;
+            pEqnComp2.ref() += fluid.dmdt()/rho2;
         }
         else
         {
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaT.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaT.H
index ff07dc4a32220dad7329838e4b87f17d5c94e762..4bef17a6bb4a851edbea09e23e8794ad3ff6cbd1 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaT.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaT.H
@@ -1,5 +1,5 @@
 {
-    volScalarField& rDeltaT = trDeltaT();
+    volScalarField& rDeltaT = trDeltaT.ref();
 
     const dictionary& pimpleDict = pimple.dict();
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/options b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/options
index 46e46e4599516204a0e140aa63b4ec4c0fd0e904..3c3db948ddad8755f7a8faceec780fbe92b2b843 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/options
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I../twoPhaseSystem/lnInclude \
     -I../../phaseSystems/lnInclude \
     -I../../interfacialModels/lnInclude\
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
index 924f8209a4c5c1affcc47ef04d81fb445c407068..4aac384117822f4242c75cb44306d621dcad60e7 100755
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -299,7 +299,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
     const scalarField liquidw(liquid.boundaryField()[patchi]);
 
     // Damp boiling at high void fractions.
-    const scalarField W(min(liquidw/0.2, scalar(0.1)));
+    const scalarField W(min(liquidw/0.2, scalar(1)));
 
     const scalarField A2(W*min(M_PI*sqr(Ddep)*N*Al/4, scalar(1)));
     const scalarField A1(max(1 - A2, scalar(1e-4)));
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 36e73662cdc05a3043055e09f18c3b96d213a6bd..597d15e7ff04f702098af467fc8830cfde9eff43 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -284,7 +284,7 @@ void Foam::twoPhaseSystem::solve()
 
         if (tdgdt.valid())
         {
-            scalarField& dgdt = tdgdt();
+            scalarField& dgdt = tdgdt.ref();
 
             forAll(dgdt, celli)
             {
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
index 2d4baf024ac9e9f95164eedfa1773551f32cdc53..c2621ae66af286cf3bff0ba0bcae87a72d650160 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
@@ -260,8 +260,8 @@ while (pimple.correct())
                   + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
-        pEqnComp1().relax();
+        deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
+        pEqnComp1.ref().relax();
 
         pEqnComp2 =
             (
@@ -276,8 +276,8 @@ while (pimple.correct())
                   + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
-        pEqnComp2().relax();
+        deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
+        pEqnComp2.ref().relax();
     }
     else
     {
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
index 9fd243f2758ac37d3acedf66c32352ff9af85a64..63321bab7d8c8835d3ea958f27986e52e3731489 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
@@ -239,8 +239,8 @@ while (pimple.correct())
                   + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
-        pEqnComp1().relax();
+        deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
+        pEqnComp1.ref().relax();
 
         pEqnComp2 =
             (
@@ -255,8 +255,8 @@ while (pimple.correct())
                   + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
                 )
             );
-        deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
-        pEqnComp2().relax();
+        deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
+        pEqnComp2.ref().relax();
     }
     else
     {
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options
index 58f3051b4d7e61bbedb05494fd571178510984b1..84786bdc5dbe4ee4ceccd265fd150dfbb5ced98b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I../twoPhaseSystem/lnInclude \
     -I../interfacialModels/lnInclude\
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
index 1cac6d7176b3229556a6345f811b8351c99d814f..1b2a492782545d0bbe21ee721507e98739d49997 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
@@ -155,17 +155,17 @@ Foam::BlendedInterfacialModel<modelType>::K() const
 
     if (model_.valid())
     {
-        x() += model_->K()*(f1() - f2());
+        x.ref() += model_->K()*(f1() - f2());
     }
 
     if (model1In2_.valid())
     {
-        x() += model1In2_->K()*(1 - f1);
+        x.ref() += model1In2_->K()*(1 - f1);
     }
 
     if (model2In1_.valid())
     {
-        x() += model2In1_->K()*f2;
+        x.ref() += model2In1_->K()*f2;
     }
 
     if
@@ -174,7 +174,7 @@ Foam::BlendedInterfacialModel<modelType>::K() const
      && (model_.valid() || model1In2_.valid() || model2In1_.valid())
     )
     {
-        correctFixedFluxBCs(x());
+        correctFixedFluxBCs(x.ref());
     }
 
     return x;
@@ -223,17 +223,17 @@ Foam::BlendedInterfacialModel<modelType>::Kf() const
 
     if (model_.valid())
     {
-        x() += model_->Kf()*(f1() - f2());
+        x.ref() += model_->Kf()*(f1() - f2());
     }
 
     if (model1In2_.valid())
     {
-        x() += model1In2_->Kf()*(1 - f1);
+        x.ref() += model1In2_->Kf()*(1 - f1);
     }
 
     if (model2In1_.valid())
     {
-        x() += model2In1_->Kf()*f2;
+        x.ref() += model2In1_->Kf()*f2;
     }
 
     if
@@ -242,7 +242,7 @@ Foam::BlendedInterfacialModel<modelType>::Kf() const
      && (model_.valid() || model1In2_.valid() || model2In1_.valid())
     )
     {
-        correctFixedFluxBCs(x());
+        correctFixedFluxBCs(x.ref());
     }
 
     return x;
@@ -286,17 +286,17 @@ Foam::BlendedInterfacialModel<modelType>::F() const
 
     if (model_.valid())
     {
-        x() += model_->F()*(f1() - f2());
+        x.ref() += model_->F()*(f1() - f2());
     }
 
     if (model1In2_.valid())
     {
-        x() += model1In2_->F()*(1 - f1);
+        x.ref() += model1In2_->F()*(1 - f1);
     }
 
     if (model2In1_.valid())
     {
-        x() -= model2In1_->F()*f2; // note : subtraction
+        x.ref() -= model2In1_->F()*f2; // note : subtraction
     }
 
     if
@@ -305,7 +305,7 @@ Foam::BlendedInterfacialModel<modelType>::F() const
      && (model_.valid() || model1In2_.valid() || model2In1_.valid())
     )
     {
-        correctFixedFluxBCs(x());
+        correctFixedFluxBCs(x.ref());
     }
 
     return x;
@@ -354,17 +354,17 @@ Foam::BlendedInterfacialModel<modelType>::Ff() const
 
     if (model_.valid())
     {
-        x() += model_->Ff()*(f1() - f2());
+        x.ref() += model_->Ff()*(f1() - f2());
     }
 
     if (model1In2_.valid())
     {
-        x() += model1In2_->Ff()*(1 - f1);
+        x.ref() += model1In2_->Ff()*(1 - f1);
     }
 
     if (model2In1_.valid())
     {
-        x() -= model2In1_->Ff()*f2; // note : subtraction
+        x.ref() -= model2In1_->Ff()*f2; // note : subtraction
     }
 
     if
@@ -373,7 +373,7 @@ Foam::BlendedInterfacialModel<modelType>::Ff() const
      && (model_.valid() || model1In2_.valid() || model2In1_.valid())
     )
     {
-        correctFixedFluxBCs(x());
+        correctFixedFluxBCs(x.ref());
     }
 
     return x;
@@ -416,17 +416,17 @@ Foam::BlendedInterfacialModel<modelType>::D() const
 
     if (model_.valid())
     {
-        x() += model_->D()*(f1() - f2());
+        x.ref() += model_->D()*(f1() - f2());
     }
 
     if (model1In2_.valid())
     {
-        x() += model1In2_->D()*(1 - f1);
+        x.ref() += model1In2_->D()*(1 - f1);
     }
 
     if (model2In1_.valid())
     {
-        x() += model2In1_->D()*f2;
+        x.ref() += model2In1_->D()*f2;
     }
 
     if
@@ -435,7 +435,7 @@ Foam::BlendedInterfacialModel<modelType>::D() const
      && (model_.valid() || model1In2_.valid() || model2In1_.valid())
     )
     {
-        correctFixedFluxBCs(x());
+        correctFixedFluxBCs(x.ref());
     }
 
     return x;
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
index 3659dcff385bea6ac9c641b95284638c1fe9e42a..ac12a39f4ef8bac81a100cb43811a564559184b2 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
@@ -1150,7 +1150,7 @@ tmp<pointField> calcOffset
     vectorField::subField fc = pp.faceCentres();
 
     tmp<pointField> toffsets(new pointField(fc.size()));
-    pointField& offsets = toffsets();
+    pointField& offsets = toffsets.ref();
 
     forAll(fc, i)
     {
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshToolsTemplates.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshToolsTemplates.C
index 04abf8a943d862bf5946a267e91081f1c378417a..143241f858371304f16a3723826b75662c47bf40 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshToolsTemplates.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshToolsTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -280,7 +280,7 @@ Foam::tmp<Foam::pointField> Foam::DelaunayMeshTools::allPoints
 )
 {
     tmp<pointField> tpts(new pointField(t.vertexCount(), point::max));
-    pointField& pts = tpts();
+    pointField& pts = tpts.ref();
 
     for
     (
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
index 301e6b8074e433549eb45b5d3bf1258c9a3e0cf9..24916fa50b34e20fdea9f19f96331da66f496bfd 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -262,7 +262,7 @@ Foam::label Foam::cellShapeControlMesh::removePoints()
 Foam::tmp<Foam::pointField> Foam::cellShapeControlMesh::cellCentres() const
 {
     tmp<pointField> tcellCentres(new pointField(number_of_finite_cells()));
-    pointField& cellCentres = tcellCentres();
+    pointField& cellCentres = tcellCentres.ref();
 
     label count = 0;
     for
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/smoothAlignmentSolver/smoothAlignmentSolver.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/smoothAlignmentSolver/smoothAlignmentSolver.C
index 456f6be9c5fa40a1650ead429a06fff10ecd5f60..af208c587b7bf98eb4b6e4d586ed2c8762b3a24e 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/smoothAlignmentSolver/smoothAlignmentSolver.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/smoothAlignmentSolver/smoothAlignmentSolver.C
@@ -179,7 +179,7 @@ Foam::tmp<Foam::triadField> Foam::smoothAlignmentSolver::buildAlignmentField
     (
         new triadField(mesh.vertexCount(), triad::unset)
     );
-    triadField& alignments = tAlignments();
+    triadField& alignments = tAlignments.ref();
 
     for
     (
@@ -211,7 +211,7 @@ Foam::tmp<Foam::pointField> Foam::smoothAlignmentSolver::buildPointField
     (
         new pointField(mesh.vertexCount(), point(GREAT, GREAT, GREAT))
     );
-    pointField& points = tPoints();
+    pointField& points = tPoints.ref();
 
     for
     (
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.C
index 5a016d88fb77779de298bece76b4684fa443f5e7..025ba5885d9cc0cc051b650996ee57b2ddc3b6ad 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -152,7 +152,7 @@ Foam::tmp<Foam::triSurfacePointScalarField> Foam::automatic::load()
         )
     );
 
-    triSurfacePointScalarField& pointCellSize = tPointCellSize();
+    triSurfacePointScalarField& pointCellSize = tPointCellSize.ref();
 
     if (readCurvature_)
     {
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C
index 805ad1196aeb311705a857312de6ac3c0ee44797..6d6b56e7b6c43de620cbeebba30c18e1002ef245 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,7 +95,7 @@ Foam::tmp<Foam::triSurfacePointScalarField> Foam::fieldFromFile::load()
         )
     );
 
-    pointCellSize() *= cellSizeMultipleCoeff_;
+    pointCellSize.ref() *= cellSizeMultipleCoeff_;
 
     return pointCellSize;
 }
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchableBoxFeatures.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchableBoxFeatures.C
index 5c3bb699e2dfa221f835a36be8bda27751eb0fff..16a18de191d0cb19f1d91eedabc4ee715a97ad07 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchableBoxFeatures.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchableBoxFeatures.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,7 +105,7 @@ Foam::searchableBoxFeatures::features() const
     edgeNormals[11][0] = 3; edgeNormals[11][1] = 0;
 
     tmp<pointField> surfacePointsTmp(surface().points());
-    pointField& surfacePoints = surfacePointsTmp();
+    pointField& surfacePoints = surfacePointsTmp.ref();
 
     forAll(edgeDirections, eI)
     {
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index 65f1f25b45792f9b52997da7434e54ff033706f5..b8e25e4f6c5b686ff8c8eb604d4d81903896973c 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,7 +84,7 @@ tmp<volScalarField> createScalarField
             zeroGradientFvPatchScalarField::typeName
         )
     );
-    volScalarField& fld = tfld();
+    volScalarField& fld = tfld.ref();
 
     forAll(fld, cellI)
     {
diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index 9494874c89598682b363e97fe5f1260a13e6fb68..6697fbecfbacea24b75fea009057ff1b6297d36f 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -181,7 +181,7 @@ void subsetVolFields
         {
             if (addedPatches.found(patchI))
             {
-                tSubFld().boundaryField()[patchI] ==
+                tSubFld.ref().boundaryField()[patchI] ==
                     pTraits<typename GeoField::value_type>::zero;
             }
         }
@@ -233,7 +233,7 @@ void subsetSurfaceFields
         {
             if (addedPatches.found(patchI))
             {
-                tSubFld().boundaryField()[patchI] ==
+                tSubFld.ref().boundaryField()[patchI] ==
                     pTraits<typename GeoField::value_type>::zero;
             }
         }
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
index e6a38a6cf6ed8fa24a5673bfa71453a8e2d1e7b0..612afc260adc59ea5510421ab810f1729424577c 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
@@ -53,8 +53,8 @@ volField
         (
             meshSubsetter.interpolate(vf)
         );
-        tfld().checkOut();
-        tfld().rename(vf.name());
+        tfld.ref().checkOut();
+        tfld.ref().rename(vf.name());
         return tfld;
     }
     else
@@ -729,7 +729,7 @@ void ensightField
         (
             volPointInterpolation::New(vf.mesh()).interpolate(vf)
         );
-        pfld().rename(vf.name());
+        pfld.ref().rename(vf.name());
 
         ensightPointField<Type>
         (
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/surfaceMeshWriterTemplates.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/surfaceMeshWriterTemplates.C
index 86ccae28675618609e9e5ebaf9fe21b9e30b32ad..c3412a1faa4d6a6738f803b07398abf546379a5d 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/surfaceMeshWriterTemplates.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/surfaceMeshWriterTemplates.C
@@ -37,7 +37,7 @@ Foam::tmp<Field<Type>> Foam::surfaceMeshWriter::getFaceField
     const polyBoundaryMesh& patches = sfld.mesh().boundaryMesh();
 
     tmp<Field<Type>> tfld(new Field<Type>(pp_.size()));
-    Field<Type>& fld = tfld();
+    Field<Type>& fld = tfld.ref();
 
     forAll(pp_.addressing(), i)
     {
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.H
index d9d65a1201219ce3afa5c3169a50642ae9858f28..f536ae04fa2f7e27faa7baf121366ba63e988fc4 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -157,7 +157,7 @@ public:
                 if (useSubMesh())
                 {
                     tmp<GeoField> subFld = subsetter_.interpolate(fld);
-                    subFld().rename(fld.name());
+                    subFld.ref().rename(fld.name());
                     return subFld;
                 }
                 else
diff --git a/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C b/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C
index 146878e272d4e8ddf44956f2f96b99ae44284e5d..2a6808776b5735bec6f5c9d64dd4fcafc006ea7c 100644
--- a/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C
+++ b/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -133,7 +133,7 @@ int main(int argc, char *argv[])
         // Calculate nut - reference nut is calculated by the turbulence model
         // on its construction
         tmp<volScalarField> tnut = turbulence->nut();
-        volScalarField& nut = tnut();
+        volScalarField& nut = tnut.ref();
         volScalarField S(mag(dev(symm(fvc::grad(U)))));
         nut = (1 - mask)*nut + mask*sqr(kappa*min(y, ybl))*::sqrt(2)*S;
 
@@ -151,7 +151,7 @@ int main(int argc, char *argv[])
 
         // Turbulence k
         tmp<volScalarField> tk = turbulence->k();
-        volScalarField& k = tk();
+        volScalarField& k = tk.ref();
         scalar ck0 = pow025(Cmu)*kappa;
         k = (1 - mask)*k + mask*sqr(nut/(ck0*min(y, ybl)));
 
@@ -165,7 +165,7 @@ int main(int argc, char *argv[])
 
         // Turbulence epsilon
         tmp<volScalarField> tepsilon = turbulence->epsilon();
-        volScalarField& epsilon = tepsilon();
+        volScalarField& epsilon = tepsilon.ref();
         scalar ce0 = ::pow(Cmu, 0.75)/kappa;
         epsilon = (1 - mask)*epsilon + mask*ce0*k*sqrt(k)/min(y, ybl);
 
diff --git a/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C b/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
index 6c1cf3cf7b44cab07ea11a375592d908086150c4..46bbfdf3952ca562639363575592bfbea6c25302 100644
--- a/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
+++ b/applications/utilities/surface/surfaceLambdaMuSmooth/surfaceLambdaMuSmooth.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,7 +57,7 @@ tmp<pointField> avg
     const labelListList& pointEdges = s.pointEdges();
 
     tmp<pointField> tavg(new pointField(s.nPoints(), vector::zero));
-    pointField& avg = tavg();
+    pointField& avg = tavg.ref();
 
     forAll(pointEdges, vertI)
     {
diff --git a/src/OpenFOAM/Make/options b/src/OpenFOAM/Make/options
index 9543b9ae27ba6fec544e4a672532d4fc5f7bc57b..b1cc1a2d10a16d24ecee7fd9ac347019fa4ed786 100644
--- a/src/OpenFOAM/Make/options
+++ b/src/OpenFOAM/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP -I$(OBJECTS_DIR)
+EXE_INC = -I$(OBJECTS_DIR)
 
 LIB_LIBS = \
     $(FOAM_LIBBIN)/libOSspecific.o \
diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.C
index c282ecb900a4fb8906d26f78ab2dd05d72f4288c..dc2a68de3897293fc0e4c78fe8f1983c30f5bcfd 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.C
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.C
@@ -60,7 +60,7 @@ pow
         )
     );
 
-    pow<Type, r, GeoMesh>(tPow().field(), df.field());
+    pow<Type, r, GeoMesh>(tPow.ref().field(), df.field());
 
     return tPow;
 }
@@ -86,7 +86,7 @@ pow
             pow(df.dimensions(), r)
         );
 
-    pow<Type, r, GeoMesh>(tPow().field(), df.field());
+    pow<Type, r, GeoMesh>(tPow.ref().field(), df.field());
 
     tdf.clear();
 
@@ -114,7 +114,7 @@ sqr(const DimensionedField<Type, GeoMesh>& df)
         )
     );
 
-    sqr(tSqr().field(), df.field());
+    sqr(tSqr.ref().field(), df.field());
 
     return tSqr;
 }
@@ -135,7 +135,7 @@ sqr(const tmp<DimensionedField<Type, GeoMesh>>& tdf)
             sqr(df.dimensions())
         );
 
-    sqr(tSqr().field(), df.field());
+    sqr(tSqr.ref().field(), df.field());
 
     tdf.clear();
 
@@ -164,7 +164,7 @@ tmp<DimensionedField<scalar, GeoMesh>> magSqr
         )
     );
 
-    magSqr(tMagSqr().field(), df.field());
+    magSqr(tMagSqr.ref().field(), df.field());
 
     return tMagSqr;
 }
@@ -185,7 +185,7 @@ tmp<DimensionedField<scalar, GeoMesh>> magSqr
             sqr(df.dimensions())
         );
 
-    magSqr(tMagSqr().field(), df.field());
+    magSqr(tMagSqr.ref().field(), df.field());
 
     tdf.clear();
 
@@ -214,7 +214,7 @@ tmp<DimensionedField<scalar, GeoMesh>> mag
         )
     );
 
-    mag(tMag().field(), df.field());
+    mag(tMag.ref().field(), df.field());
 
     return tMag;
 }
@@ -235,7 +235,7 @@ tmp<DimensionedField<scalar, GeoMesh>> mag
             df.dimensions()
         );
 
-    mag(tMag().field(), df.field());
+    mag(tMag.ref().field(), df.field());
 
     tdf.clear();
 
@@ -268,7 +268,7 @@ cmptAv(const DimensionedField<Type, GeoMesh>& df)
         )
     );
 
-    cmptAv(CmptAv().field(), df.field());
+    cmptAv(CmptAv.ref().field(), df.field());
 
     return CmptAv;
 }
@@ -294,7 +294,7 @@ cmptAv(const tmp<DimensionedField<Type, GeoMesh>>& tdf)
             df.dimensions()
         );
 
-    cmptAv(CmptAv().field(), df.field());
+    cmptAv(CmptAv.ref().field(), df.field());
 
     tdf.clear();
 
diff --git a/src/OpenFOAM/memory/tmp/tmp.H b/src/OpenFOAM/memory/tmp/tmp.H
index 4b320d8e05db0898cfecaab913cddf819aa34dea..9afb78bcfe8c468481249eee61975932806d141d 100644
--- a/src/OpenFOAM/memory/tmp/tmp.H
+++ b/src/OpenFOAM/memory/tmp/tmp.H
@@ -128,7 +128,7 @@ public:
 
     // Member operators
 
-        #ifndef CONST_TMP
+        #ifdef NON_CONST_TMP
         //- Deprecated non-const dereference operator.
         //  Use ref() where non-const access is required
         inline T& operator()();
diff --git a/src/OpenFOAM/memory/tmp/tmpI.H b/src/OpenFOAM/memory/tmp/tmpI.H
index 5c962200487e1172d4f7aaab09d41c278daedf9c..819c4714284fe64f5e3028d176187751282b37e4 100644
--- a/src/OpenFOAM/memory/tmp/tmpI.H
+++ b/src/OpenFOAM/memory/tmp/tmpI.H
@@ -217,7 +217,7 @@ inline void Foam::tmp<T>::clear() const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-#ifndef CONST_TMP
+#ifdef NON_CONST_TMP
 template<class T>
 inline T& Foam::tmp<T>::operator()()
 {
diff --git a/src/TurbulenceModels/compressible/Make/options b/src/TurbulenceModels/compressible/Make/options
index a6401963fd1764080a335633f0d5e3bdb93f298c..4a6578d628571d492b84de2023e9bc1995fae803 100644
--- a/src/TurbulenceModels/compressible/Make/options
+++ b/src/TurbulenceModels/compressible/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I../turbulenceModels/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
diff --git a/src/TurbulenceModels/incompressible/Make/options b/src/TurbulenceModels/incompressible/Make/options
index cf54ffd963bbcf322562cf402445e4f0621665fd..c44e9ae61c9cc89f7984bcebd9cef0318f44de8e 100644
--- a/src/TurbulenceModels/incompressible/Make/options
+++ b/src/TurbulenceModels/incompressible/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I../turbulenceModels/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/src/TurbulenceModels/turbulenceModels/Make/options b/src/TurbulenceModels/turbulenceModels/Make/options
index bcd627471b97e93547ed691c52bd70b90be9030e..a3ae8da833177387e9eecf75b5e2675fc7b481f5 100644
--- a/src/TurbulenceModels/turbulenceModels/Make/options
+++ b/src/TurbulenceModels/turbulenceModels/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
diff --git a/src/combustionModels/Make/options b/src/combustionModels/Make/options
index d1ee86bab4b2ce7cdf645b1922decd7c77821eb0..e660ccd9f6163b86cd418aada44cdc8dd9d3dde1 100644
--- a/src/combustionModels/Make/options
+++ b/src/combustionModels/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
diff --git a/src/conversion/Make/options b/src/conversion/Make/options
index ace5a3133b72a7e5eb692293380a51b32212b24a..4ff461186ccfa78b8effcbd50daebaf0f1007595 100644
--- a/src/conversion/Make/options
+++ b/src/conversion/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
diff --git a/src/dynamicFvMesh/Make/options b/src/dynamicFvMesh/Make/options
index 3cee7d2a3c498d13e88f730fee2413da5e075f08..966b56964d720cb3acf3a2f006a63d6b14284ef5 100644
--- a/src/dynamicFvMesh/Make/options
+++ b/src/dynamicFvMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
diff --git a/src/dynamicMesh/Make/options b/src/dynamicMesh/Make/options
index 5725e0f867119f09546876956fd17e6ee3b42b6d..0757c296783738fafdd35fc63fafb035a5fc1500 100644
--- a/src/dynamicMesh/Make/options
+++ b/src/dynamicMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
diff --git a/src/engine/Make/options b/src/engine/Make/options
index f16d25abf155c0b4f3cbd83b8219d1f51c506a8f..0929e2bc8f16444baec5e7d3b051594b0a20e031 100644
--- a/src/engine/Make/options
+++ b/src/engine/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
diff --git a/src/finiteVolume/Make/options b/src/finiteVolume/Make/options
index 25eb0cc6f8f5a02ce5c78579b32ebe96d682ab66..1c1990e5f3e5f624194a3d86a9cd7f3f30ce2aca 100644
--- a/src/finiteVolume/Make/options
+++ b/src/finiteVolume/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C b/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C
index b02f92354f7696b24d22cd843071a03bcee8c825..46a3f5db9d91b11e84579adebab4b7a5ea795b5b 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C
@@ -51,7 +51,7 @@ d2dt2
     (
         vf.mesh(),
         vf.mesh().d2dt2Scheme("d2dt2(" + vf.name() + ')')
-    )().fvmD2dt2(vf);
+    ).ref().fvmD2dt2(vf);
 }
 
 
@@ -67,7 +67,7 @@ d2dt2
     (
         vf.mesh(),
         vf.mesh().d2dt2Scheme("d2dt2(" + rho.name() + ',' + vf.name() + ')')
-    )().fvmD2dt2(rho, vf);
+    ).ref().fvmD2dt2(rho, vf);
 }
 
 
@@ -83,7 +83,7 @@ d2dt2
     (
         vf.mesh(),
         vf.mesh().d2dt2Scheme("d2dt2(" + rho.name() + ',' + vf.name() + ')')
-    )().fvmD2dt2(rho, vf);
+    ).ref().fvmD2dt2(rho, vf);
 }
 
 
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 91d7384f552e56e1bca8341053e363bf1339060c..799b5d2198863b7bfab9a0a72d3fdb19522d777c 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -1401,7 +1401,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::correction
      && A.psi().mesh().fluxRequired(A.psi().name())
     )
     {
-        tAcorr().faceFluxCorrectionPtr() = (-A.flux()).ptr();
+        tAcorr.ref().faceFluxCorrectionPtr() = (-A.flux()).ptr();
     }
 
     return tAcorr;
diff --git a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
index 33255e5c14df80817e903575f4f537ad42fbd303..7f3923aef8fbede3ac390f94bc696c48db396df2 100644
--- a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
+++ b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
@@ -406,7 +406,7 @@ fvMeshSubset::interpolate
             patchFields
         )
     );
-    GeometricField<Type, pointPatchField, pointMesh>& resF = tresF();
+    GeometricField<Type, pointPatchField, pointMesh>& resF = tresF.ref();
 
 
     // 2. Change the pointPatchFields to the correct type using a mapper
diff --git a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMeshInterpolate.C b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMeshInterpolate.C
index 186fe7d500ee5da0e6db7f591a9d35ede642fc37..171e47ecd4409fcac99ce989e0a6677f941ec7fd 100644
--- a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMeshInterpolate.C
+++ b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMeshInterpolate.C
@@ -77,7 +77,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh>> singleCellFvMesh::interpolate
             patchFields
         )
     );
-    GeometricField<Type, fvPatchField, volMesh>& resF = tresF();
+    GeometricField<Type, fvPatchField, volMesh>& resF = tresF.ref();
 
 
     // 2. Change the fvPatchFields to the correct type using a mapper
diff --git a/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options b/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options
index 63d132ca5072f60e4e205fc7099de406bb65babe..93a7e01542e3a516c7ed5370f0f7b4b77b8e9544 100644
--- a/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options
+++ b/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/OpenFOAM/lnInclude
diff --git a/src/fvMotionSolver/Make/options b/src/fvMotionSolver/Make/options
index e0b34247ad2c52cba20f58a8125f4a2e1e5181bf..29e4469610118cb4f0468cba7f80f5aacbf7a1fd 100644
--- a/src/fvMotionSolver/Make/options
+++ b/src/fvMotionSolver/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
diff --git a/src/fvOptions/Make/options b/src/fvOptions/Make/options
index 5598b4ef20f3cd991f42df32c26b8443fe65d879..124084c0575dfd4a0fb89d7c272c891ca1ea1fd6 100644
--- a/src/fvOptions/Make/options
+++ b/src/fvOptions/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
diff --git a/src/genericPatchFields/Make/options b/src/genericPatchFields/Make/options
index 287318da1ff31d9067bc255cafbdd1fa8dce4787..71b7873964d544eddf96d22aa40f4c3372c23c9c 100644
--- a/src/genericPatchFields/Make/options
+++ b/src/genericPatchFields/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
diff --git a/src/lagrangian/DSMC/Make/options b/src/lagrangian/DSMC/Make/options
index 4f10470b16005977acb0493671334966914c4446..55865dfabcd9146aa4fe204e71d664b62e535b95 100644
--- a/src/lagrangian/DSMC/Make/options
+++ b/src/lagrangian/DSMC/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
diff --git a/src/lagrangian/basic/Make/options b/src/lagrangian/basic/Make/options
index 4670e80487454a6007985424b1dc7fb49fbc2c93..f7d9a36f3b6448a99c7aaf6aecdcaefd1a967f5b 100644
--- a/src/lagrangian/basic/Make/options
+++ b/src/lagrangian/basic/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
diff --git a/src/lagrangian/coalCombustion/Make/options b/src/lagrangian/coalCombustion/Make/options
index 9c78fa29f774e03a17324a318e9b3142072b444d..04b133ba05321d3e8e29079eba7602c5e8744ccc 100644
--- a/src/lagrangian/coalCombustion/Make/options
+++ b/src/lagrangian/coalCombustion/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
diff --git a/src/lagrangian/distributionModels/Make/options b/src/lagrangian/distributionModels/Make/options
index ccd5450e02b21a31066f39bed97f5c517c629a49..41306609f208806f0c6f42a2426867d3e10d4897 100644
--- a/src/lagrangian/distributionModels/Make/options
+++ b/src/lagrangian/distributionModels/Make/options
@@ -1 +1 @@
-EXE_INC = -DCONST_TMP
+EXE_INC =
diff --git a/src/lagrangian/intermediate/Make/options b/src/lagrangian/intermediate/Make/options
index 18eeb8197017cc190adfe6a5692a0b3959df2f23..a5a3b3c8ebdb14fa2fd410aabc9fe451f0e84b4c 100644
--- a/src/lagrangian/intermediate/Make/options
+++ b/src/lagrangian/intermediate/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
diff --git a/src/lagrangian/molecularDynamics/molecule/Make/options b/src/lagrangian/molecularDynamics/molecule/Make/options
index 6ccac8015b0371444752e141f4421b4a4d0d0d3d..de95d1b75b80f1bfb1d186901e8eee45e6e29951 100644
--- a/src/lagrangian/molecularDynamics/molecule/Make/options
+++ b/src/lagrangian/molecularDynamics/molecule/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
diff --git a/src/lagrangian/solidParticle/Make/options b/src/lagrangian/solidParticle/Make/options
index 5cf0265987001e1a59307e01b82b1d13dc7cbac2..9fe79c9220b60af2b5c87cc7e9f6974096372b1a 100644
--- a/src/lagrangian/solidParticle/Make/options
+++ b/src/lagrangian/solidParticle/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude
diff --git a/src/lagrangian/spray/Make/options b/src/lagrangian/spray/Make/options
index d5eb8a74c56200b6c8503c4a042d1436d312b2f3..3b5a4e72e9565bfce6a407b8f8641dcffa5201b5 100644
--- a/src/lagrangian/spray/Make/options
+++ b/src/lagrangian/spray/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
diff --git a/src/lagrangian/turbulence/Make/options b/src/lagrangian/turbulence/Make/options
index 67f1b6608c5efb4516d60cec9fa443ce08355e40..7b2fc4bea88fb6df9d69954051ba03a7af2872e6 100644
--- a/src/lagrangian/turbulence/Make/options
+++ b/src/lagrangian/turbulence/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
diff --git a/src/mesh/autoMesh/Make/options b/src/mesh/autoMesh/Make/options
index 78bd91e5936551e82b92425fbae0231c4d32ceba..d37a26a9dd460742d9203623e4ddf6f69a4c99de 100644
--- a/src/mesh/autoMesh/Make/options
+++ b/src/mesh/autoMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/src/mesh/blockMesh/Make/options b/src/mesh/blockMesh/Make/options
index d398d1687fc4a2bec5520e06d00c7206b573f316..1618ab57ec2d869d439de31fd8ead10e6f2b8ae2 100644
--- a/src/mesh/blockMesh/Make/options
+++ b/src/mesh/blockMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
 
diff --git a/src/mesh/extrudeModel/Make/options b/src/mesh/extrudeModel/Make/options
index 92abbacf67a1b6c6ed2e20d4a294673f357006dd..eabd0b53a8ffd8a16a228d71ae038fac7e2ea6a2 100644
--- a/src/mesh/extrudeModel/Make/options
+++ b/src/mesh/extrudeModel/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
 
diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options
index 2f8d1a7b478857c6e18402dff04058f3a3fdb892..1713152e9e502c6789d21eefb44ffd4904e08f91 100644
--- a/src/meshTools/Make/options
+++ b/src/meshTools/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude
diff --git a/src/parallel/decompose/decompose/fvFieldDecomposerDecomposeFields.C b/src/parallel/decompose/decompose/fvFieldDecomposerDecomposeFields.C
index 24a3b377dae7cf1e6b6e177e1ced2d67047abd32..65875e91fbef82627f3f5311e114416d21dc2ec0 100644
--- a/src/parallel/decompose/decompose/fvFieldDecomposerDecomposeFields.C
+++ b/src/parallel/decompose/decompose/fvFieldDecomposerDecomposeFields.C
@@ -76,7 +76,7 @@ Foam::fvFieldDecomposer::decomposeField
             patchFields
         )
     );
-    GeometricField<Type, fvPatchField, volMesh>& resF = tresF();
+    GeometricField<Type, fvPatchField, volMesh>& resF = tresF.ref();
 
 
     // 2. Change the fvPatchFields to the correct type using a mapper
@@ -246,7 +246,7 @@ Foam::fvFieldDecomposer::decomposeField
             patchFields
         )
     );
-    GeometricField<Type, fvsPatchField, surfaceMesh>& resF = tresF();
+    GeometricField<Type, fvsPatchField, surfaceMesh>& resF = tresF.ref();
 
 
     // 2. Change the fvsPatchFields to the correct type using a mapper
diff --git a/src/parallel/reconstruct/reconstruct/reconstructLagrangianFields.C b/src/parallel/reconstruct/reconstruct/reconstructLagrangianFields.C
index 1fb0a26ac06d772d72e50b14360f6cc721b45c8b..960a65d76a8bf0b64e90b1e8f0468fc2a9d80b30 100644
--- a/src/parallel/reconstruct/reconstruct/reconstructLagrangianFields.C
+++ b/src/parallel/reconstruct/reconstruct/reconstructLagrangianFields.C
@@ -55,7 +55,7 @@ Foam::tmp<Foam::IOField<Type>> Foam::reconstructLagrangianField
             Field<Type>(0)
         )
     );
-    Field<Type>& field = tfield();
+    Field<Type>& field = tfield.ref();
 
     forAll(meshes, i)
     {
@@ -115,7 +115,7 @@ Foam::reconstructLagrangianFieldField
             Field<Field<Type>>(0)
         )
     );
-    Field<Field<Type>>& field = tfield();
+    Field<Field<Type>>& field = tfield.ref();
 
     forAll(meshes, i)
     {
diff --git a/src/postProcessing/functionObjects/field/Make/options b/src/postProcessing/functionObjects/field/Make/options
index 583c538b9e62e6355c5925432239a9c0da048f38..0732822dba6309b4f07c986b07406c8eece747b1 100644
--- a/src/postProcessing/functionObjects/field/Make/options
+++ b/src/postProcessing/functionObjects/field/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
diff --git a/src/postProcessing/functionObjects/utilities/Make/options b/src/postProcessing/functionObjects/utilities/Make/options
index f50a857be26bbfbdd36450ae5c271b74a7c279b2..ad4e27c1967426d1d3fdfdeb9b45fa67744396d6 100644
--- a/src/postProcessing/functionObjects/utilities/Make/options
+++ b/src/postProcessing/functionObjects/utilities/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/DSMC/lnInclude \
     -I$(LIB_SRC)/transportModels \
diff --git a/src/randomProcesses/Make/options b/src/randomProcesses/Make/options
index 287318da1ff31d9067bc255cafbdd1fa8dce4787..71b7873964d544eddf96d22aa40f4c3372c23c9c 100644
--- a/src/randomProcesses/Make/options
+++ b/src/randomProcesses/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
diff --git a/src/regionCoupled/Make/options b/src/regionCoupled/Make/options
index 21a2486c48c15020941cff955e95ea0ec38f35b4..a9753c5115e98af376181d88ac517ec012762443 100644
--- a/src/regionCoupled/Make/options
+++ b/src/regionCoupled/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
diff --git a/src/regionModels/regionModel/Make/options b/src/regionModels/regionModel/Make/options
index bcd627471b97e93547ed691c52bd70b90be9030e..a3ae8da833177387e9eecf75b5e2675fc7b481f5 100644
--- a/src/regionModels/regionModel/Make/options
+++ b/src/regionModels/regionModel/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
diff --git a/src/regionModels/surfaceFilmModels/Make/options b/src/regionModels/surfaceFilmModels/Make/options
index 0cd68b47e811caae0f788214546a2d504ba7db3a..a0eb4828eaa8393b4c2a6c690113957d7bd4359b 100644
--- a/src/regionModels/surfaceFilmModels/Make/options
+++ b/src/regionModels/surfaceFilmModels/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
diff --git a/src/sampling/Make/options b/src/sampling/Make/options
index 15ca6b932f8ba3574029ca7ec291cad669d9b7a3..23e6bc81a7e708d69b5fd14d27073827646e14ef 100644
--- a/src/sampling/Make/options
+++ b/src/sampling/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
diff --git a/src/sixDoFRigidBodyMotion/Make/options b/src/sixDoFRigidBodyMotion/Make/options
index fd1f401690d3eb816fb956ec3114c8b18d0d4e00..f2367d2e1d47e8b701978c06a03ad5696c8c8668 100644
--- a/src/sixDoFRigidBodyMotion/Make/options
+++ b/src/sixDoFRigidBodyMotion/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/postProcessing/functionObjects/forces/lnInclude \
diff --git a/src/surfMesh/Make/options b/src/surfMesh/Make/options
index 079607db058a7cc9acdcca3a7aa9b0f20a0a36a0..7e207d0dbeaca258e5a72af8b4eb7bacefc0dee8 100644
--- a/src/surfMesh/Make/options
+++ b/src/surfMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/fileFormats/lnInclude
 
 LIB_LIBS = \
diff --git a/src/thermophysicalModels/SLGThermo/Make/options b/src/thermophysicalModels/SLGThermo/Make/options
index 3dac244d79650919f883a3d0123751ab278f633a..c4f0b2a8a95fa2420cb555ebad3e0f96387baca4 100644
--- a/src/thermophysicalModels/SLGThermo/Make/options
+++ b/src/thermophysicalModels/SLGThermo/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
diff --git a/src/thermophysicalModels/barotropicCompressibilityModel/Make/options b/src/thermophysicalModels/barotropicCompressibilityModel/Make/options
index 287318da1ff31d9067bc255cafbdd1fa8dce4787..71b7873964d544eddf96d22aa40f4c3372c23c9c 100644
--- a/src/thermophysicalModels/barotropicCompressibilityModel/Make/options
+++ b/src/thermophysicalModels/barotropicCompressibilityModel/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
diff --git a/src/thermophysicalModels/basic/Make/options b/src/thermophysicalModels/basic/Make/options
index 95f178b455134e96eb3b0c998c3bab5429182cdc..b5c859baf1457da06327d60192105d4b8fbaad53 100644
--- a/src/thermophysicalModels/basic/Make/options
+++ b/src/thermophysicalModels/basic/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/src/thermophysicalModels/chemistryModel/Make/options b/src/thermophysicalModels/chemistryModel/Make/options
index 7a1e7f43d2b62cfdf78057946e1634fa19f1b922..de52d7e6ae8750b7e48f9e7dc740e3c2da538042 100644
--- a/src/thermophysicalModels/chemistryModel/Make/options
+++ b/src/thermophysicalModels/chemistryModel/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
diff --git a/src/thermophysicalModels/laminarFlameSpeed/Make/options b/src/thermophysicalModels/laminarFlameSpeed/Make/options
index 925d4a58aab37f576d449679c9073303acf44228..f97edb8458c3d408db331fb5aa10930eadd2fc73 100644
--- a/src/thermophysicalModels/laminarFlameSpeed/Make/options
+++ b/src/thermophysicalModels/laminarFlameSpeed/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
diff --git a/src/thermophysicalModels/properties/liquidMixtureProperties/Make/options b/src/thermophysicalModels/properties/liquidMixtureProperties/Make/options
index 690761730fada84af3821927c633c9bdd8f77e9b..6e1d19dbb60b095df0f8c87b45e1d7cdf4449c9e 100644
--- a/src/thermophysicalModels/properties/liquidMixtureProperties/Make/options
+++ b/src/thermophysicalModels/properties/liquidMixtureProperties/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
diff --git a/src/thermophysicalModels/properties/liquidProperties/Make/options b/src/thermophysicalModels/properties/liquidProperties/Make/options
index e1fbd3ed45354782b62c37ba685a67beea951632..b964b61294c787fe650b03faee1b09bbb9b48256 100644
--- a/src/thermophysicalModels/properties/liquidProperties/Make/options
+++ b/src/thermophysicalModels/properties/liquidProperties/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude
 
 LIB_LIBS = \
diff --git a/src/thermophysicalModels/properties/solidMixtureProperties/Make/options b/src/thermophysicalModels/properties/solidMixtureProperties/Make/options
index 97cc50921233c7a9d0cee3aab90633b64e1b17e6..f2d8f809d5d12dac81a95c25e5a85a21cc9adac1 100644
--- a/src/thermophysicalModels/properties/solidMixtureProperties/Make/options
+++ b/src/thermophysicalModels/properties/solidMixtureProperties/Make/options
@@ -1,3 +1,3 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I${LIB_SRC}/thermophysicalModels/properties/solidProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
diff --git a/src/thermophysicalModels/properties/solidProperties/Make/options b/src/thermophysicalModels/properties/solidProperties/Make/options
index f8779f538f675e83239fbc61890c5653da29ec84..848cce789f2aad33d9ed6eda979ea31c32377425 100644
--- a/src/thermophysicalModels/properties/solidProperties/Make/options
+++ b/src/thermophysicalModels/properties/solidProperties/Make/options
@@ -1,2 +1,2 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
diff --git a/src/thermophysicalModels/radiation/Make/options b/src/thermophysicalModels/radiation/Make/options
index ea189b0b63517b6784f7380df21d424531b5b4b9..726b76e76442077c0a0c4f77334e54d9bfe6a6ad 100644
--- a/src/thermophysicalModels/radiation/Make/options
+++ b/src/thermophysicalModels/radiation/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
diff --git a/src/thermophysicalModels/reactionThermo/Make/options b/src/thermophysicalModels/reactionThermo/Make/options
index bded2d7118f6d62b9e7b73c14ac287c60c1b4b9b..f59f44fc8d3a1d08df9534646ba44212a3034ba2 100644
--- a/src/thermophysicalModels/reactionThermo/Make/options
+++ b/src/thermophysicalModels/reactionThermo/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
diff --git a/src/thermophysicalModels/solidChemistryModel/Make/options b/src/thermophysicalModels/solidChemistryModel/Make/options
index e922aa724f5f1fa2eb49dd3b39eaebacfa6f38e3..e871169b17f52e46e9fdf1ab5878f7c4f4532acc 100644
--- a/src/thermophysicalModels/solidChemistryModel/Make/options
+++ b/src/thermophysicalModels/solidChemistryModel/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
diff --git a/src/thermophysicalModels/solidSpecie/Make/options b/src/thermophysicalModels/solidSpecie/Make/options
index 85a93d1d170df8cdd0fe01f707948075b80fe868..bdc8b074c03a509ce2f001f40aff6c0c189b0222 100644
--- a/src/thermophysicalModels/solidSpecie/Make/options
+++ b/src/thermophysicalModels/solidSpecie/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
 
 LIB_LIBS = \
diff --git a/src/thermophysicalModels/solidThermo/Make/options b/src/thermophysicalModels/solidThermo/Make/options
index c84b6c55464449cf1fd5e8853db822c7a96f41c2..d3e1cb87ce5eac43f9599cf03d7d2c861d24bfab 100644
--- a/src/thermophysicalModels/solidThermo/Make/options
+++ b/src/thermophysicalModels/solidThermo/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
diff --git a/src/topoChangerFvMesh/Make/options b/src/topoChangerFvMesh/Make/options
index 82165ae6e5ec39ca860ac3ef18b61976d679f3ab..44753e64c96546ae740c2d54605959f8ba6890fd 100644
--- a/src/topoChangerFvMesh/Make/options
+++ b/src/topoChangerFvMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
diff --git a/src/transportModels/compressible/Make/options b/src/transportModels/compressible/Make/options
index 287318da1ff31d9067bc255cafbdd1fa8dce4787..71b7873964d544eddf96d22aa40f4c3372c23c9c 100644
--- a/src/transportModels/compressible/Make/options
+++ b/src/transportModels/compressible/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
diff --git a/src/transportModels/immiscibleIncompressibleTwoPhaseMixture/Make/options b/src/transportModels/immiscibleIncompressibleTwoPhaseMixture/Make/options
index 138dcbf7f55763cb937eb1952e92913782558664..4ae2fbf96644b3ac1443ae0bea87926356d158fd 100644
--- a/src/transportModels/immiscibleIncompressibleTwoPhaseMixture/Make/options
+++ b/src/transportModels/immiscibleIncompressibleTwoPhaseMixture/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I.. \
     -I../incompressible/lnInclude \
     -I../interfaceProperties/lnInclude \
diff --git a/src/transportModels/incompressible/Make/options b/src/transportModels/incompressible/Make/options
index c55a98a0b64a5facfe106252a638c7a3772b4271..46aa26d4ad8980819f6d6b939101f3f3fef92f64 100644
--- a/src/transportModels/incompressible/Make/options
+++ b/src/transportModels/incompressible/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I.. \
     -I../twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/src/transportModels/interfaceProperties/Make/options b/src/transportModels/interfaceProperties/Make/options
index efbf1b41ca3431db1b4527a8b8f188a73ad239f7..cc63db7dcc45c854ad4a57862b1b9ad635ec5ce7 100644
--- a/src/transportModels/interfaceProperties/Make/options
+++ b/src/transportModels/interfaceProperties/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels/twoPhaseProperties/alphaContactAngle/alphaContactAngle \
     -I$(LIB_SRC)/finiteVolume/lnInclude
diff --git a/src/transportModels/twoPhaseMixture/Make/options b/src/transportModels/twoPhaseMixture/Make/options
index 287318da1ff31d9067bc255cafbdd1fa8dce4787..71b7873964d544eddf96d22aa40f4c3372c23c9c 100644
--- a/src/transportModels/twoPhaseMixture/Make/options
+++ b/src/transportModels/twoPhaseMixture/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
diff --git a/src/transportModels/twoPhaseProperties/Make/options b/src/transportModels/twoPhaseProperties/Make/options
index 547d45694c84a3d980b5c1b88dc36aaea3d2e224..32ce36b9594d8e970ba45526b49fce90d0d9fef0 100644
--- a/src/transportModels/twoPhaseProperties/Make/options
+++ b/src/transportModels/twoPhaseProperties/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
diff --git a/src/triSurface/Make/options b/src/triSurface/Make/options
index 948c54f7560fdc0d13aef5830827846f1d0dd459..9ee5884e5908ccdb6b1ba3acf323b8c524ec996a 100644
--- a/src/triSurface/Make/options
+++ b/src/triSurface/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = -DCONST_TMP \
+EXE_INC = \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude