diff --git a/applications/test/fieldMapping/Test-fieldMapping.C b/applications/test/fieldMapping/Test-fieldMapping.C
index 38a4b110f96446249cc9d041aeb9b088bcfb4f33..41af9c79231cafcf614d708e0a547182ebe61105 100644
--- a/applications/test/fieldMapping/Test-fieldMapping.C
+++ b/applications/test/fieldMapping/Test-fieldMapping.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
@@ -39,6 +39,7 @@ Description
 #include "mapPolyMesh.H"
 #include "polyTopoChange.H"
 #include "fvcDiv.H"
+#include "zeroGradientFvPatchFields.H"
 #include "Random.H"
 
 using namespace Foam;
diff --git a/src/OpenFOAM/fields/Fields/Field/Field.C b/src/OpenFOAM/fields/Fields/Field/Field.C
index 596f74f1c30c33138275bde95debf41064b577de..e7e7479c1528326908307e22c4b85d676a188e59 100644
--- a/src/OpenFOAM/fields/Fields/Field/Field.C
+++ b/src/OpenFOAM/fields/Fields/Field/Field.C
@@ -609,7 +609,7 @@ template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::Field<Type>::T() const
 {
     tmp<Field<Type>> transpose(new Field<Type>(this->size()));
-    ::Foam::T(transpose(), *this);
+    ::Foam::T(transpose.ref(), *this);
     return transpose;
 }
 
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C
index e53fdf778a6ca601aa5daae3d6385969c18aa882..a70abece0ffd6fea1d21a481989c36a91527f8f3 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/uniformInterpolate.C
@@ -92,7 +92,7 @@ Foam::tmp<GeoField> Foam::uniformInterpolate
 
     // Interpolate
     tmp<GeoField> tfld(new GeoField(fieldIO, weights[0]*field0));
-    GeoField& fld = tfld();
+    GeoField& fld = tfld.ref();
 
     for (label i = 1; i < times.size(); ++i)
     {
diff --git a/src/OpenFOAM/interpolations/primitivePatchInterpolation/PrimitivePatchInterpolation.C b/src/OpenFOAM/interpolations/primitivePatchInterpolation/PrimitivePatchInterpolation.C
index fd77c9a677892d85101bc15a746dda5eaf780d4c..e9253134ef47e7e9796565cdc6c38cb7dee032b9 100644
--- a/src/OpenFOAM/interpolations/primitivePatchInterpolation/PrimitivePatchInterpolation.C
+++ b/src/OpenFOAM/interpolations/primitivePatchInterpolation/PrimitivePatchInterpolation.C
@@ -192,7 +192,7 @@ tmp<Field<Type>> PrimitivePatchInterpolation<Patch>::faceToPointInterpolate
         )
     );
 
-    Field<Type>& result = tresult();
+    Field<Type>& result = tresult.ref();
 
     const labelListList& pointFaces = patch_.pointFaces();
     const scalarListList& weights = faceToPointWeights();
@@ -249,7 +249,7 @@ tmp<Field<Type>> PrimitivePatchInterpolation<Patch>::pointToFaceInterpolate
         )
     );
 
-    Field<Type>& result = tresult();
+    Field<Type>& result = tresult.ref();
 
     const List<typename Patch::FaceType>& localFaces = patch_.localFaces();
 
@@ -303,7 +303,7 @@ tmp<Field<Type>> PrimitivePatchInterpolation<Patch>::faceToEdgeInterpolate
         new Field<Type>(patch_.nEdges(), pTraits<Type>::zero)
     );
 
-    Field<Type>& result = tresult();
+    Field<Type>& result = tresult.ref();
 
     const edgeList& edges = patch_.edges();
     const labelListList& edgeFaces = patch_.edgeFaces();
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C
index 82f5851f3308950d18b3f8063eb82ee8aa66be80..34be2afdcae74417b98e479dae4571aa5d3d3980 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C
@@ -159,7 +159,7 @@ Foam::PatchTools::pointNormals
     //    to avoid them being stored)
 
     tmp<pointField> textrudeN(new pointField(p.nPoints(), vector::zero));
-    pointField& extrudeN = textrudeN();
+    pointField& extrudeN = textrudeN.ref();
     {
         const faceList& localFaces = p.localFaces();
         const vectorField& faceNormals = p.faceNormals();
@@ -213,7 +213,7 @@ Foam::PatchTools::edgeNormals
     // 1. Start off with local normals
 
     tmp<pointField> tedgeNormals(new pointField(p.nEdges(), vector::zero));
-    pointField& edgeNormals = tedgeNormals();
+    pointField& edgeNormals = tedgeNormals.ref();
     {
         const labelListList& edgeFaces = p.edgeFaces();
         const vectorField& faceNormals = p.faceNormals();
diff --git a/src/combustionModels/FSD/FSD.C b/src/combustionModels/FSD/FSD.C
index 6d846b3d95e07d811ed4d1e539bc6a44f8749e94..8f10a8940ce5c7c39282cdc117919c64f959687d 100644
--- a/src/combustionModels/FSD/FSD.C
+++ b/src/combustionModels/FSD/FSD.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
@@ -161,7 +161,7 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
         )
     );
 
-    volScalarField& pc = tPc();
+    volScalarField& pc = tPc.ref();
 
     tmp<volScalarField> tomegaFuel
     (
@@ -185,7 +185,7 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
         )
     );
 
-    volScalarField& omegaFuelBar = tomegaFuel();
+    volScalarField& omegaFuelBar = tomegaFuel.ref();
 
     // Calculation of the mixture fraction variance (ftVar)
     const compressible::LESModel& lesModel =
@@ -311,7 +311,7 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
         )
     );
 
-    volScalarField& products = tproducts();
+    volScalarField& products = tproducts.ref();
 
     forAll(productsIndex, j)
     {
diff --git a/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C b/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
index 2a5cbbd7dafb8ed25b8379fc190d107be5f5a77c..1aa0e0ae83ed4c38282ba5e4af80d1ee1270155c 100644
--- a/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
+++ b/src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.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
@@ -102,7 +102,7 @@ Foam::tmp<Foam::volScalarField> Foam::consumptionSpeed::omega0Sigma
         )
     );
 
-    volScalarField& omega0 = tomega0();
+    volScalarField& omega0 = tomega0.ref();
 
     volScalarField::InternalField& iomega0 = omega0.internalField();
 
diff --git a/src/combustionModels/Make/options b/src/combustionModels/Make/options
index e660ccd9f6163b86cd418aada44cdc8dd9d3dde1..d1ee86bab4b2ce7cdf645b1922decd7c77821eb0 100644
--- a/src/combustionModels/Make/options
+++ b/src/combustionModels/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
diff --git a/src/combustionModels/laminar/laminar.C b/src/combustionModels/laminar/laminar.C
index d53a859a733d07a8a944210ac09fd0f8ef2965ba..155c8bf01e3d43d0bf905b84c30687ceef089bdb 100644
--- a/src/combustionModels/laminar/laminar.C
+++ b/src/combustionModels/laminar/laminar.C
@@ -119,7 +119,7 @@ Foam::combustionModels::laminar<Type>::R(volScalarField& Y) const
 {
     tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
 
-    fvScalarMatrix& Su = tSu();
+    fvScalarMatrix& Su = tSu.ref();
 
     if (this->active())
     {
@@ -157,7 +157,7 @@ Foam::combustionModels::laminar<Type>::dQ() const
 
     if (this->active())
     {
-        tdQ() = this->chemistryPtr_->dQ();
+        tdQ.ref() = this->chemistryPtr_->dQ();
     }
 
     return tdQ;
@@ -188,7 +188,7 @@ Foam::combustionModels::laminar<Type>::Sh() const
 
     if (this->active())
     {
-        tSh() = this->chemistryPtr_->Sh();
+        tSh.ref() = this->chemistryPtr_->Sh();
     }
 
     return tSh;
diff --git a/src/combustionModels/singleStepCombustion/singleStepCombustion.C b/src/combustionModels/singleStepCombustion/singleStepCombustion.C
index ba586bc1cd42ac5abb4e01e586c36db96863f500..d8bc1cfac30a0bcc59404ed732d7d0cb6cc571f3 100644
--- a/src/combustionModels/singleStepCombustion/singleStepCombustion.C
+++ b/src/combustionModels/singleStepCombustion/singleStepCombustion.C
@@ -160,7 +160,7 @@ singleStepCombustion<CombThermoType, ThermoType>::dQ() const
 
     if (this->active())
     {
-        volScalarField& dQ = tdQ();
+        volScalarField& dQ = tdQ.ref();
         dQ.dimensionedInternalField() = this->mesh().V()*Sh()();
     }
     return tdQ;
diff --git a/src/conversion/Make/options b/src/conversion/Make/options
index 4ff461186ccfa78b8effcbd50daebaf0f1007595..ace5a3133b72a7e5eb692293380a51b32212b24a 100644
--- a/src/conversion/Make/options
+++ b/src/conversion/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 966b56964d720cb3acf3a2f006a63d6b14284ef5..3cee7d2a3c498d13e88f730fee2413da5e075f08 100644
--- a/src/dynamicFvMesh/Make/options
+++ b/src/dynamicFvMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 0757c296783738fafdd35fc63fafb035a5fc1500..5725e0f867119f09546876956fd17e6ee3b42b6d 100644
--- a/src/dynamicMesh/Make/options
+++ b/src/dynamicMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
diff --git a/src/dynamicMesh/layerAdditionRemoval/addCellLayer.C b/src/dynamicMesh/layerAdditionRemoval/addCellLayer.C
index 572b8b7d790a2bf093fb23111bfef2bbbc74fed6..085cc99079e6b008dcfbed37b8356ba2ffac28ab 100644
--- a/src/dynamicMesh/layerAdditionRemoval/addCellLayer.C
+++ b/src/dynamicMesh/layerAdditionRemoval/addCellLayer.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
@@ -45,7 +45,7 @@ Foam::tmp<Foam::vectorField> Foam::layerAdditionRemoval::extrusionDir() const
     const labelList& mp = masterFaceLayer.meshPoints();
 
     tmp<vectorField> textrusionDir(new vectorField(mp.size()));
-    vectorField& extrusionDir = textrusionDir();
+    vectorField& extrusionDir = textrusionDir.ref();
 
     if (setLayerPairing())
     {
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
index bebb9ea3aa8ab229a465a4ba592c084ec8e41156..384663939e8b3eb86da321035b513cfcad679951 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
@@ -119,7 +119,7 @@ Foam::tmp<Foam::scalarField> Foam::motionSmootherAlgo::calcEdgeWeights
     const edgeList& edges = mesh_.edges();
 
     tmp<scalarField> twght(new scalarField(edges.size()));
-    scalarField& wght = twght();
+    scalarField& wght = twght.ref();
 
     forAll(edges, edgeI)
     {
@@ -815,7 +815,7 @@ Foam::tmp<Foam::pointField> Foam::motionSmootherAlgo::curPoints() const
     tmp<pointField> tnewPoints(oldPoints_ + totalDisplacement.internalField());
 
     // Correct for 2-D motion
-    modifyMotionPoints(tnewPoints());
+    modifyMotionPoints(tnewPoints.ref());
 
     return tnewPoints;
 }
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
index f59703c56f3bb62593333dfb06346adc2c017a6a..f7a9174d8e9663bd88f267f31d6c8b4a39f35b1d 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
@@ -159,7 +159,7 @@ Foam::motionSmootherAlgo::avg
             dimensioned<Type>("zero", fld.dimensions(), pTraits<Type>::zero)
         )
     );
-    GeometricField<Type, pointPatchField, pointMesh>& res = tres();
+    GeometricField<Type, pointPatchField, pointMesh>& res = tres.ref();
 
     const polyMesh& mesh = fld.mesh()();
 
diff --git a/src/engine/Make/options b/src/engine/Make/options
index 0929e2bc8f16444baec5e7d3b051594b0a20e031..f16d25abf155c0b4f3cbd83b8219d1f51c506a8f 100644
--- a/src/engine/Make/options
+++ b/src/engine/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcCellReduce.C b/src/finiteVolume/finiteVolume/fvc/fvcCellReduce.C
index 3ec972fc2abc97b9cbb5cd8676f12cc923716181..c078bd365d2c96eb0aab4cf3eecce4433f2f2e7b 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcCellReduce.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcCellReduce.C
@@ -70,7 +70,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh>> cellReduce
         )
     );
 
-    volFieldType& result = tresult();
+    volFieldType& result = tresult.ref();
 
     const labelUList& own = mesh.owner();
     const labelUList& nbr = mesh.neighbour();
@@ -102,7 +102,8 @@ tmp<GeometricField<Type, fvPatchField, volMesh>> cellReduce
     tmp<GeometricField<Type, fvPatchField, volMesh>>
         tvf(cellReduce(cop, tssf));
 
-   tssf.clear();
+    tssf.clear();
+
     return tvf;
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcCurl.C b/src/finiteVolume/finiteVolume/fvc/fvcCurl.C
index 38dee36681d650674d132cc0ae1e76cb97ea2b9a..1ec6024fb9b41db8fdd1eab0be532ab787edbdb3 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcCurl.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcCurl.C
@@ -56,7 +56,7 @@ curl
     tmp<GeometricField<Type, fvPatchField, volMesh>> tcurlVf =
         2.0*(*skew(fvc::grad(vf, nameCurlVf)));
 
-    tcurlVf().rename(nameCurlVf);
+    tcurlVf.ref().rename(nameCurlVf);
 
     return tcurlVf;
 }
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
index 9b2f375b0865751bcfc547162a498eb69fb7d9a7..32c8a51dda769fbe68617068856e7206ed239d64 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
@@ -78,7 +78,7 @@ reconstruct
         )
     );
 
-    treconField().correctBoundaryConditions();
+    treconField.ref().correctBoundaryConditions();
 
     return treconField;
 }
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 6d4b18f3714a159dac7bf45db556254ff72dfc8b..91d7384f552e56e1bca8341053e363bf1339060c 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -693,7 +693,7 @@ template<class Type>
 Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
 {
     tmp<scalarField> tdiag(new scalarField(diag()));
-    addCmptAvBoundaryDiag(tdiag());
+    addCmptAvBoundaryDiag(tdiag.ref());
     return tdiag;
 }
 
@@ -713,7 +713,7 @@ Foam::tmp<Foam::Field<Type>> Foam::fvMatrix<Type>::DD() const
             (
                 lduAddr().patchAddr(patchI),
                 internalCoeffs_[patchI],
-                tdiag()
+                tdiag.ref()
             );
         }
     }
diff --git a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
index dbf2c09d7753064876990e4999be16b0659ce281..33255e5c14df80817e903575f4f537ad42fbd303 100644
--- a/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
+++ b/src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C
@@ -100,7 +100,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh>> fvMeshSubset::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
@@ -243,7 +243,7 @@ tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> fvMeshSubset::interpolate
             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/fvAgglomerationMethods/pairPatchAgglomeration/Make/options b/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options
index 93a7e01542e3a516c7ed5370f0f7b4b77b8e9544..63d132ca5072f60e4e205fc7099de406bb65babe 100644
--- a/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options
+++ b/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/OpenFOAM/lnInclude
diff --git a/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C
index ab49334a5d942c5b029f38034a19e0021bfa90be..32bf9efa12bb85432342565777c285da39328056 100644
--- a/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C
+++ b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C
@@ -440,7 +440,7 @@ Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
     const label nFineFaces = patch.size();
 
     tmp<labelField> tcoarseCellMap(new labelField(nFineFaces, -1));
-    labelField& coarseCellMap = tcoarseCellMap();
+    labelField& coarseCellMap = tcoarseCellMap.ref();
 
     const labelListList& faceFaces = patch.faceFaces();
 
diff --git a/src/fvMotionSolver/Make/options b/src/fvMotionSolver/Make/options
index 29e4469610118cb4f0468cba7f80f5aacbf7a1fd..e0b34247ad2c52cba20f58a8125f4a2e1e5181bf 100644
--- a/src/fvMotionSolver/Make/options
+++ b/src/fvMotionSolver/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
diff --git a/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
index 8f90cb8444538f561b40325c94cb858b67876470..fe9ac195c73655b0cee59a9308fea0f4b0285c13 100644
--- a/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -188,8 +188,9 @@ Foam::displacementComponentLaplacianFvMotionSolver::curPoints() const
     else
     {
         tmp<pointField> tcurPoints(new pointField(fvMesh_.points()));
+        pointField& curPoints = tcurPoints.ref();
 
-        tcurPoints().replace
+        curPoints.replace
         (
             cmpt_,
             points0_ + pointDisplacement_.internalField()
@@ -204,11 +205,11 @@ Foam::displacementComponentLaplacianFvMotionSolver::curPoints() const
             {
                 label pointI = pz[i];
 
-                tcurPoints()[pointI][cmpt_] = points0_[pointI];
+                curPoints[pointI][cmpt_] = points0_[pointI];
             }
         }
 
-        twoDCorrectPoints(tcurPoints());
+        twoDCorrectPoints(curPoints);
 
         return tcurPoints;
     }
diff --git a/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C
index 124a679e1543db698307c02efa6d29fad12b6012..bb40c74881edc4454f763c465f818e473057be48 100644
--- a/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -101,14 +101,14 @@ Foam::velocityComponentLaplacianFvMotionSolver::curPoints() const
 
     tmp<pointField> tcurPoints(new pointField(fvMesh_.points()));
 
-    tcurPoints().replace
+    tcurPoints.ref().replace
     (
         cmpt_,
         tcurPoints().component(cmpt_)
       + fvMesh_.time().deltaTValue()*pointMotionU_.internalField()
     );
 
-    twoDCorrectPoints(tcurPoints());
+    twoDCorrectPoints(tcurPoints.ref());
 
     return tcurPoints;
 }
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
index 7aafefc0c52dba4c5e0bfee22bf5c812dd2fca3b..ba28c423c9e9188e9fa61f7253762eaea71c711c 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.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::displacementSBRStressFvMotionSolver::curPoints() const
         points0() + pointDisplacement().internalField()
     );
 
-    twoDCorrectPoints(tcurPoints());
+    twoDCorrectPoints(tcurPoints.ref());
 
     return tcurPoints;
 }
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C
index f3286af26a4b6cb636e9212a1ef000ef3b7b4df6..b1a72eb6f80902e3aeb34a811d41c91dc4cbfa8f 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C
@@ -313,7 +313,7 @@ Foam::displacementInterpolationMotionSolver::curPoints() const
     }
 
     tmp<pointField> tcurPoints(new pointField(points0()));
-    pointField& curPoints = tcurPoints();
+    pointField& curPoints = tcurPoints.ref();
 
     // Interpolate the displacement of the face zones.
     vectorField zoneDisp(displacements_.size(), vector::zero);
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
index d4a3ca5a0a7c1f7b3dfdd2fe55eb3879d26a68b9..863f800d190e17384609c1988eb637ff11d73326 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -200,6 +200,7 @@ Foam::displacementLaplacianFvMotionSolver::curPoints() const
         (
             points0() + pointDisplacement_.internalField()
         );
+        pointField& curPoints = tcurPoints.ref();
 
         // Implement frozen points
         if (frozenPointsZone_ != -1)
@@ -208,11 +209,11 @@ Foam::displacementLaplacianFvMotionSolver::curPoints() const
 
             forAll(pz, i)
             {
-                tcurPoints()[pz[i]] = points0()[pz[i]];
+                curPoints[pz[i]] = points0()[pz[i]];
             }
         }
 
-        twoDCorrectPoints(tcurPoints());
+        twoDCorrectPoints(curPoints);
 
         return tcurPoints;
     }
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
index 832f249e20265666b709423e16e611f1a4f0dc66..0900344564ac8731879876ce9969e70d71a6e294 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.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
@@ -223,7 +223,7 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
 ) const
 {
     tmp<vectorField> tfld(new vectorField(meshPoints.size()));
-    vectorField& fld = tfld();
+    vectorField& fld = tfld.ref();
 
     const word type(dict.lookup("type"));
 
diff --git a/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
index ce63274ed9fc8a25dfab738bde02878932a521ba..ba05811ba6dcbb0092eeccef56897d7939fb135a 100644
--- a/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,7 +103,7 @@ Foam::velocityLaplacianFvMotionSolver::curPoints() const
       + fvMesh_.time().deltaTValue()*pointMotionU_.internalField()
     );
 
-    twoDCorrectPoints(tcurPoints());
+    twoDCorrectPoints(tcurPoints.ref());
 
     return tcurPoints;
 }
diff --git a/src/fvOptions/Make/options b/src/fvOptions/Make/options
index 124084c0575dfd4a0fb89d7c272c891ca1ea1fd6..5598b4ef20f3cd991f42df32c26b8443fe65d879 100644
--- a/src/fvOptions/Make/options
+++ b/src/fvOptions/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
diff --git a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
index 181db2575e4549a5662000c52966153f15735269..3ca231a5c1ddde0bbad0b72c5bbbc3817e546ecc 100644
--- a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
+++ b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.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
@@ -186,7 +186,7 @@ void Foam::fv::rotorDiskSource::writeField
             )
         );
 
-        Field<Type>& field = tfield().internalField();
+        Field<Type>& field = tfield.ref().internalField();
 
         if (cells_.size() != values.size())
         {
diff --git a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
index 7ac47f1b2b2e6626f3aee28fadbfc2e2f407092d..1f9905b94bd6c1c298d8fac7436f99db3d606a7b 100644
--- a/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
+++ b/src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -260,7 +260,7 @@ Foam::tmp<Foam::scalarField> Foam::targetCoeffTrim::thetag() const
     const List<vector>& x = rotor_.x();
 
     tmp<scalarField> ttheta(new scalarField(x.size()));
-    scalarField& t = ttheta();
+    scalarField& t = ttheta.ref();
 
     forAll(t, i)
     {
diff --git a/src/fvOptions/sources/interRegion/interRegionHeatTransfer/interRegionHeatTransferModel/interRegionHeatTransferModel.C b/src/fvOptions/sources/interRegion/interRegionHeatTransfer/interRegionHeatTransferModel/interRegionHeatTransferModel.C
index 170f72517ba967766ed22287a185f2ec177b21cc..ee213ae246f4718ee07a742e656c43b32ff86388 100644
--- a/src/fvOptions/sources/interRegion/interRegionHeatTransfer/interRegionHeatTransferModel/interRegionHeatTransferModel.C
+++ b/src/fvOptions/sources/interRegion/interRegionHeatTransfer/interRegionHeatTransferModel/interRegionHeatTransferModel.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
@@ -194,7 +194,7 @@ void Foam::fv::interRegionHeatTransferModel::addSup
         )
     );
 
-    volScalarField& Tmapped = tTmapped();
+    volScalarField& Tmapped = tTmapped.ref();
 
     const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);
 
diff --git a/src/genericPatchFields/Make/options b/src/genericPatchFields/Make/options
index 71b7873964d544eddf96d22aa40f4c3372c23c9c..287318da1ff31d9067bc255cafbdd1fa8dce4787 100644
--- a/src/genericPatchFields/Make/options
+++ b/src/genericPatchFields/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
diff --git a/src/lagrangian/DSMC/Make/options b/src/lagrangian/DSMC/Make/options
index 55865dfabcd9146aa4fe204e71d664b62e535b95..4f10470b16005977acb0493671334966914c4446 100644
--- a/src/lagrangian/DSMC/Make/options
+++ b/src/lagrangian/DSMC/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 f7d9a36f3b6448a99c7aaf6aecdcaefd1a967f5b..4670e80487454a6007985424b1dc7fb49fbc2c93 100644
--- a/src/lagrangian/basic/Make/options
+++ b/src/lagrangian/basic/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 04b133ba05321d3e8e29079eba7602c5e8744ccc..9c78fa29f774e03a17324a318e9b3142072b444d 100644
--- a/src/lagrangian/coalCombustion/Make/options
+++ b/src/lagrangian/coalCombustion/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
diff --git a/src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H b/src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H
index 16152ff9dc464dea452e30fa54650910155a4e5e..a249a875d1b67c3f5f09fbe47407a23cc668041a 100644
--- a/src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H
+++ b/src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H
@@ -47,7 +47,7 @@ Foam::coalCloudList::UTrans() const
         )
     );
 
-    DimensionedField<vector, volMesh>& fld = tfld();
+    DimensionedField<vector, volMesh>& fld = tfld.ref();
 
     forAll(*this, i)
     {
@@ -64,7 +64,7 @@ Foam::tmp<Foam::fvVectorMatrix> Foam::coalCloudList::SU
 ) const
 {
     tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
-    fvVectorMatrix& fvm = tfvm();
+    fvVectorMatrix& fvm = tfvm.ref();
 
     forAll(*this, i)
     {
@@ -95,7 +95,7 @@ Foam::coalCloudList::hsTrans() const
         )
     );
 
-    DimensionedField<scalar, volMesh>& fld = tfld();
+    DimensionedField<scalar, volMesh>& fld = tfld.ref();
 
     forAll(*this, i)
     {
@@ -112,7 +112,7 @@ Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Sh
 ) const
 {
     tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(hs, dimEnergy/dimTime));
-    fvScalarMatrix& fvm = tfvm();
+    fvScalarMatrix& fvm = tfvm.ref();
 
     forAll(*this, i)
     {
@@ -130,7 +130,7 @@ Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::SYi
 ) const
 {
     tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(Yi, dimMass/dimTime));
-    fvScalarMatrix& fvm = tfvm();
+    fvScalarMatrix& fvm = tfvm.ref();
 
     forAll(*this, i)
     {
@@ -161,7 +161,7 @@ Foam::coalCloudList::rhoTrans() const
         )
     );
 
-    DimensionedField<scalar, volMesh>& fld = tfld();
+    DimensionedField<scalar, volMesh>& fld = tfld.ref();
 
     forAll(*this, i)
     {
@@ -197,7 +197,7 @@ Foam::coalCloudList::Srho() const
         )
     );
 
-    DimensionedField<scalar, volMesh>& fld = tfld();
+    DimensionedField<scalar, volMesh>& fld = tfld.ref();
 
     forAll(*this, i)
     {
@@ -231,7 +231,7 @@ Foam::coalCloudList::Srho
         )
     );
 
-    DimensionedField<scalar, volMesh>& fld = tfld();
+    DimensionedField<scalar, volMesh>& fld = tfld.ref();
 
     forAll(*this, j)
     {
@@ -248,7 +248,7 @@ Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Srho
 ) const
 {
     tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(rho, dimMass/dimTime));
-    fvScalarMatrix& fvm = tfvm();
+    fvScalarMatrix& fvm = tfvm.ref();
 
     forAll(*this, i)
     {
diff --git a/src/lagrangian/distributionModels/Make/options b/src/lagrangian/distributionModels/Make/options
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..ccd5450e02b21a31066f39bed97f5c517c629a49 100644
--- a/src/lagrangian/distributionModels/Make/options
+++ b/src/lagrangian/distributionModels/Make/options
@@ -0,0 +1 @@
+EXE_INC = -DCONST_TMP
diff --git a/src/lagrangian/intermediate/Make/options b/src/lagrangian/intermediate/Make/options
index a5a3b3c8ebdb14fa2fd410aabc9fe451f0e84b4c..18eeb8197017cc190adfe6a5692a0b3959df2f23 100644
--- a/src/lagrangian/intermediate/Make/options
+++ b/src/lagrangian/intermediate/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index 3ed90408757bd21ad4c78f8747d186aa571de23a..69b0eb60af82104de726298640db1e3b85f0ff40 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -439,7 +439,7 @@ Foam::KinematicCloud<CloudType>::SU(volVectorField& U) const
         else
         {
             tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
-            fvVectorMatrix& fvm = tfvm();
+            fvVectorMatrix& fvm = tfvm.ref();
 
             fvm.source() = -UTrans()/(this->db().time().deltaT());
 
@@ -474,7 +474,7 @@ Foam::KinematicCloud<CloudType>::vDotSweep() const
         )
     );
 
-    volScalarField& vDotSweep = tvDotSweep();
+    volScalarField& vDotSweep = tvDotSweep.ref();
     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
@@ -513,7 +513,7 @@ Foam::KinematicCloud<CloudType>::theta() const
         )
     );
 
-    volScalarField& theta = ttheta();
+    volScalarField& theta = ttheta.ref();
     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
@@ -551,7 +551,7 @@ Foam::KinematicCloud<CloudType>::alpha() const
         )
     );
 
-    scalarField& alpha = talpha().internalField();
+    scalarField& alpha = talpha.ref().internalField();
     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
@@ -588,7 +588,7 @@ Foam::KinematicCloud<CloudType>::rhoEff() const
         )
     );
 
-    scalarField& rhoEff = trhoEff().internalField();
+    scalarField& rhoEff = trhoEff.ref().internalField();
     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
index d15938121505a72e0a46c92d732ece78531170bc..1c9907ecdb3702c48ea0982143789a96d802eba8 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
@@ -127,7 +127,7 @@ inline Foam::tmp<Foam::fvScalarMatrix> Foam::ReactingCloud<CloudType>::SYi
                 )
             );
 
-            volScalarField& sourceField = trhoTrans();
+            volScalarField& sourceField = trhoTrans.ref();
 
             sourceField.internalField() =
                 rhoTrans_[i]/(this->db().time().deltaTValue()*this->mesh().V());
@@ -141,7 +141,7 @@ inline Foam::tmp<Foam::fvScalarMatrix> Foam::ReactingCloud<CloudType>::SYi
         else
         {
             tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(Yi, dimMass/dimTime));
-            fvScalarMatrix& fvm = tfvm();
+            fvScalarMatrix& fvm = tfvm.ref();
 
             fvm.source() = -rhoTrans_[i]/this->db().time().deltaTValue();
 
@@ -182,7 +182,7 @@ Foam::ReactingCloud<CloudType>::Srho(const label i) const
 
     if (this->solution().coupled())
     {
-        scalarField& rhoi = tRhoi();
+        scalarField& rhoi = tRhoi.ref();
         rhoi = rhoTrans_[i]/(this->db().time().deltaTValue()*this->mesh().V());
     }
 
@@ -219,7 +219,7 @@ Foam::ReactingCloud<CloudType>::Srho() const
 
     if (this->solution().coupled())
     {
-        scalarField& sourceField = trhoTrans();
+        scalarField& sourceField = trhoTrans.ref();
         forAll(rhoTrans_, i)
         {
             sourceField += rhoTrans_[i];
@@ -256,7 +256,7 @@ Foam::ReactingCloud<CloudType>::Srho(volScalarField& rho) const
             )
         );
 
-        scalarField& sourceField = trhoTrans();
+        scalarField& sourceField = trhoTrans.ref();
 
         if (this->solution().semiImplicit("rho"))
         {
@@ -272,7 +272,7 @@ Foam::ReactingCloud<CloudType>::Srho(volScalarField& rho) const
         else
         {
             tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(rho, dimMass/dimTime));
-            fvScalarMatrix& fvm = tfvm();
+            fvScalarMatrix& fvm = tfvm.ref();
 
             forAll(rhoTrans_, i)
             {
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
index 9655ef95f36ef7c45da79a2ee6c6f2e294e3e558..d23371cb8213b3b9d03e3bed796e9ad560b24920 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
@@ -247,7 +247,7 @@ Foam::ThermoCloud<CloudType>::Sh(volScalarField& hs) const
         else
         {
             tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(hs, dimEnergy/dimTime));
-            fvScalarMatrix& fvm = tfvm();
+            fvScalarMatrix& fvm = tfvm.ref();
 
             fvm.source() = -hsTrans()/(this->db().time().deltaT());
 
@@ -282,7 +282,7 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::Ep() const
 
     if (radiation_)
     {
-        scalarField& Ep = tEp().internalField();
+        scalarField& Ep = tEp.ref().internalField();
         const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
@@ -318,7 +318,7 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::ap() const
 
     if (radiation_)
     {
-        scalarField& ap = tap().internalField();
+        scalarField& ap = tap.ref().internalField();
         const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
@@ -355,7 +355,7 @@ Foam::ThermoCloud<CloudType>::sigmap() const
 
     if (radiation_)
     {
-        scalarField& sigmap = tsigmap().internalField();
+        scalarField& sigmap = tsigmap.ref().internalField();
         const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
index 73a72f5540c10e82a741e12227ea8c53182dfddc..ad22dc25fb69b3cd52dd921856f3aac285dcfd12 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
@@ -267,14 +267,14 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
 
             if (applyGravity_)
             {
-                phiCorrect_() -= phiGByA();
+                phiCorrect_.ref() -= phiGByA();
             }
 
             forAll(phiCorrect_(), faceI)
             {
                 // Current and correction fluxes
                 const scalar phiCurr = phi[faceI];
-                scalar& phiCorr = phiCorrect_()[faceI];
+                scalar& phiCorr = phiCorrect_.ref()[faceI];
 
                 // Don't limit if the correction is in the opposite direction to
                 // the flux. We need all the help we can get in this state.
@@ -296,7 +296,7 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
 
             if (applyGravity_)
             {
-                phiCorrect_() += phiGByA();
+                phiCorrect_.ref() += phiGByA();
             }
         }
 
diff --git a/src/lagrangian/intermediate/submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C b/src/lagrangian/intermediate/submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
index 56153ff88bcd4f3c19cddef902815a0e28726370..bf3a35272bd489105f436f38b01e1661bee706f0 100644
--- a/src/lagrangian/intermediate/submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
+++ b/src/lagrangian/intermediate/submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.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
@@ -95,7 +95,7 @@ Foam::radiation::cloudAbsorptionEmission::aDisp(const label) const
             mesh_.objectRegistry::lookupObject<thermoCloud>(cloudNames_[i])
         );
 
-        ta() += tc.ap();
+        ta.ref() += tc.ap();
     }
 
     return ta;
@@ -155,7 +155,7 @@ Foam::radiation::cloudAbsorptionEmission::EDisp(const label bandI) const
             mesh_.objectRegistry::lookupObject<thermoCloud>(cloudNames_[i])
         );
 
-        tE() += tc.Ep();
+        tE.ref() += tc.Ep();
     }
 
     // Total emission is 4 times the projected emission
diff --git a/src/lagrangian/intermediate/submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C b/src/lagrangian/intermediate/submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
index aa38692a769249afbc881bf6310787b2daf60027..0f88407b4e7b5cbd332e40e0426978a2200f9046 100644
--- a/src/lagrangian/intermediate/submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
+++ b/src/lagrangian/intermediate/submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,7 +95,7 @@ Foam::radiation::cloudScatter::sigmaEff() const
             mesh_.objectRegistry::lookupObject<thermoCloud>(cloudNames_[i])
         );
 
-        tsigma() += tc.sigmap();
+        tsigma.ref() += tc.sigmap();
     }
 
     return 3.0*tsigma;
diff --git a/src/lagrangian/molecularDynamics/molecule/Make/options b/src/lagrangian/molecularDynamics/molecule/Make/options
index 5698f81a651a3be1d4ee7a32fe2c7867f3a343af..6ccac8015b0371444752e141f4421b4a4d0d0d3d 100644
--- a/src/lagrangian/molecularDynamics/molecule/Make/options
+++ b/src/lagrangian/molecularDynamics/molecule/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
@@ -11,4 +11,3 @@ LIB_LIBS = \
     -llagrangian \
     -lpotential \
     -lmolecularMeasurements
-
diff --git a/src/lagrangian/solidParticle/Make/options b/src/lagrangian/solidParticle/Make/options
index 9fe79c9220b60af2b5c87cc7e9f6974096372b1a..5cf0265987001e1a59307e01b82b1d13dc7cbac2 100644
--- a/src/lagrangian/solidParticle/Make/options
+++ b/src/lagrangian/solidParticle/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 3b5a4e72e9565bfce6a407b8f8641dcffa5201b5..d5eb8a74c56200b6c8503c4a042d1436d312b2f3 100644
--- a/src/lagrangian/spray/Make/options
+++ b/src/lagrangian/spray/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 7b2fc4bea88fb6df9d69954051ba03a7af2872e6..67f1b6608c5efb4516d60cec9fa443ce08355e40 100644
--- a/src/lagrangian/turbulence/Make/options
+++ b/src/lagrangian/turbulence/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 d37a26a9dd460742d9203623e4ddf6f69a4c99de..78bd91e5936551e82b92425fbae0231c4d32ceba 100644
--- a/src/mesh/autoMesh/Make/options
+++ b/src/mesh/autoMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index f1727c44b0a84d99f66602f94fd6a5f51384165d..91ce9cce7594269f117b3c46305056cfd7157267 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -112,7 +112,7 @@ Foam::tmp<Foam::scalarField> Foam::autoLayerDriver::avgPointData
 )
 {
     tmp<scalarField> tfaceFld(new scalarField(pp.size(), 0.0));
-    scalarField& faceFld = tfaceFld();
+    scalarField& faceFld = tfaceFld.ref();
 
     forAll(pp.localFaces(), faceI)
     {
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 7e4fd3a1af58c27d8a2a04a844f454468a620488..29475a3958b2688eca91aadec30af5f61b9476d5 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.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
@@ -514,7 +514,7 @@ Foam::tmp<Foam::scalarField> Foam::autoSnapDriver::edgePatchDist
 
     // Copy edge values into scalarField
     tmp<scalarField> tedgeDist(new scalarField(mesh.nEdges()));
-    scalarField& edgeDist = tedgeDist();
+    scalarField& edgeDist = tedgeDist.ref();
 
     forAll(allEdgeInfo, edgeI)
     {
@@ -874,7 +874,7 @@ Foam::tmp<Foam::pointField> Foam::autoSnapDriver::avgCellCentres
     (
         new pointField(pointFaces.size(), vector::zero)
     );
-    pointField& avgBoundary = tavgBoundary();
+    pointField& avgBoundary = tavgBoundary.ref();
     labelList nBoundary(pointFaces.size(), 0);
 
     forAll(pointFaces, pointI)
diff --git a/src/mesh/blockMesh/Make/options b/src/mesh/blockMesh/Make/options
index 1618ab57ec2d869d439de31fd8ead10e6f2b8ae2..d398d1687fc4a2bec5520e06d00c7206b573f316 100644
--- a/src/mesh/blockMesh/Make/options
+++ b/src/mesh/blockMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 eabd0b53a8ffd8a16a228d71ae038fac7e2ea6a2..92abbacf67a1b6c6ed2e20d4a294673f357006dd 100644
--- a/src/mesh/extrudeModel/Make/options
+++ b/src/mesh/extrudeModel/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
 
diff --git a/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.C b/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.C
index d9e3a0fb77ea458d335683192945561b6de0ee82..50dad9b28b83234ec4974095ba1d098f1742693b 100644
--- a/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.C
+++ b/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -184,7 +184,7 @@ Foam::cyclicACMIGAMGInterface::internalFieldTransfer
     const labelUList& nbrFaceCells = nbr.faceCells();
 
     tmp<labelField> tpnf(new labelField(nbrFaceCells.size()));
-    labelField& pnf = tpnf();
+    labelField& pnf = tpnf.ref();
 
     forAll(pnf, facei)
     {
diff --git a/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C b/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C
index 151b021aba86b9d3ec74c188df743c767ad33e34..0270d4c645d2052aaba6a8642e2c3c6bba314120 100644
--- a/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C
+++ b/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -183,7 +183,7 @@ Foam::tmp<Foam::labelField> Foam::cyclicAMIGAMGInterface::internalFieldTransfer
     const labelUList& nbrFaceCells = nbr.faceCells();
 
     tmp<labelField> tpnf(new labelField(nbrFaceCells.size()));
-    labelField& pnf = tpnf();
+    labelField& pnf = tpnf.ref();
 
     forAll(pnf, facei)
     {
diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options
index 1713152e9e502c6789d21eefb44ffd4904e08f91..2f8d1a7b478857c6e18402dff04058f3a3fdb892 100644
--- a/src/meshTools/Make/options
+++ b/src/meshTools/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude
diff --git a/src/meshTools/cellQuality/cellQuality.C b/src/meshTools/cellQuality/cellQuality.C
index 95fe7e53399487da905660c1782dbf91793d75d9..c1165ee375b147b981c720da5821a7ff2b258020 100644
--- a/src/meshTools/cellQuality/cellQuality.C
+++ b/src/meshTools/cellQuality/cellQuality.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,7 +47,7 @@ Foam::tmp<Foam::scalarField> Foam::cellQuality::nonOrthogonality() const
         )
     );
 
-    scalarField& result = tresult();
+    scalarField& result = tresult.ref();
 
     scalarField sumArea(mesh_.nCells(), 0.0);
 
@@ -108,7 +108,7 @@ Foam::tmp<Foam::scalarField> Foam::cellQuality::skewness() const
             mesh_.nCells(), 0.0
         )
     );
-    scalarField& result = tresult();
+    scalarField& result = tresult.ref();
 
     scalarField sumArea(mesh_.nCells(), 0.0);
 
@@ -187,7 +187,7 @@ Foam::tmp<Foam::scalarField> Foam::cellQuality::faceNonOrthogonality() const
             mesh_.nFaces(), 0.0
         )
     );
-    scalarField& result = tresult();
+    scalarField& result = tresult.ref();
 
 
     const vectorField& centres = mesh_.cellCentres();
@@ -247,7 +247,7 @@ Foam::tmp<Foam::scalarField> Foam::cellQuality::faceSkewness() const
             mesh_.nFaces(), 0.0
         )
     );
-    scalarField& result = tresult();
+    scalarField& result = tresult.ref();
 
 
     const vectorField& cellCtrs = mesh_.cellCentres();
diff --git a/src/meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.C b/src/meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.C
index 0ffe75d2e657ba40a62a791e49d9c91224d6a2e1..a0bfb14b8fc907164e4c8981aa318a086e575c77 100644
--- a/src/meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.C
+++ b/src/meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.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
@@ -128,7 +128,7 @@ transformVector
 ) const
 {
     tmp<symmTensorField> tfld(new symmTensorField(st.size()));
-    symmTensorField& fld = tfld();
+    symmTensorField& fld = tfld.ref();
 
     forAll(fld, i)
     {
diff --git a/src/meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.C b/src/meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.C
index fde7eb39cc270e744967d6575e081c57082929d4..41e331c2b5d7d09141b8fa200a2096b746d54da3 100644
--- a/src/meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.C
+++ b/src/meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.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
@@ -129,7 +129,7 @@ transformVector
 ) const
 {
     tmp<symmTensorField> tfld(new symmTensorField(st.size()));
-    symmTensorField& fld = tfld();
+    symmTensorField& fld = tfld.ref();
 
     forAll(fld, i)
     {
diff --git a/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.C b/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.C
index c28a6fd3b480896e7c82a74ac977700569b147d2..ae0614f8de1063b4c84b8bb7de462eacdbab21e2 100644
--- a/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.C
+++ b/src/meshTools/coordinateSystems/coordinateRotation/axesRotation.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
@@ -232,7 +232,7 @@ Foam::tmp<Foam::symmTensorField> Foam::axesRotation::transformVector
 ) const
 {
     tmp<symmTensorField> tfld(new symmTensorField(st.size()));
-    symmTensorField& fld = tfld();
+    symmTensorField& fld = tfld.ref();
 
     forAll(fld, i)
     {
diff --git a/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C b/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C
index a5427b5fa1e3acf2dfb2706626a2e36e687afacf..572db42283598781d566c7f9d446fe0401670933 100644
--- a/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C
+++ b/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.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
@@ -304,7 +304,7 @@ Foam::tmp<Foam::tensorField> Foam::cylindrical::transformTensor
     const tensorField& R = Rptr_();
     const tensorField Rtr(R.T());
     tmp<tensorField> tt(new tensorField(cellMap.size()));
-    tensorField& t = tt();
+    tensorField& t = tt.ref();
     forAll(cellMap, i)
     {
         const label cellI = cellMap[i];
@@ -328,7 +328,7 @@ Foam::tmp<Foam::symmTensorField> Foam::cylindrical::transformVector
     }
 
     tmp<symmTensorField> tfld(new symmTensorField(Rptr_->size()));
-    symmTensorField& fld = tfld();
+    symmTensorField& fld = tfld.ref();
 
     const tensorField& R = Rptr_();
     forAll(fld, i)
diff --git a/src/meshTools/coordinateSystems/cylindricalCS.C b/src/meshTools/coordinateSystems/cylindricalCS.C
index 18be4ab3a8a8e8816b28b4d70012e71fb4b286ab..fda455545ac1e325adb458d0e66e8c9a1c82b40b 100644
--- a/src/meshTools/coordinateSystems/cylindricalCS.C
+++ b/src/meshTools/coordinateSystems/cylindricalCS.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -208,7 +208,7 @@ Foam::tmp<Foam::vectorField> Foam::cylindricalCS::globalToLocal
     );
 
     tmp<vectorField> tresult(new vectorField(lc.size()));
-    vectorField& result = tresult();
+    vectorField& result = tresult.ref();
 
     result.replace
     (
diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
index 22e92fa7c364f8fb380a09b1121e6d5dd46ecc89..0ecb78004f440c55df6a5e6f60f9f6656b5b1c1a 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
@@ -99,7 +99,7 @@ Foam::tmp<Foam::pointField> Foam::mappedPatchBase::facePoints
 
     // Initialise to face-centre
     tmp<pointField> tfacePoints(new pointField(patch_.size()));
-    pointField& facePoints = tfacePoints();
+    pointField& facePoints = tfacePoints.ref();
 
     forAll(pp, faceI)
     {
@@ -867,7 +867,7 @@ Foam::tmp<Foam::pointField> Foam::mappedPatchBase::readListOrField
 )
 {
     tmp<pointField> tfld(new pointField());
-    pointField& fld = tfld();
+    pointField& fld = tfld.ref();
 
     if (size)
     {
@@ -1267,7 +1267,7 @@ Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints
 ) const
 {
     tmp<pointField> tfld(new pointField(fc));
-    pointField& fld = tfld();
+    pointField& fld = tfld.ref();
 
     switch (offsetMode_)
     {
diff --git a/src/meshTools/momentOfInertia/momentOfInertia.C b/src/meshTools/momentOfInertia/momentOfInertia.C
index bb17f5833cc6a8b28a20ef2b1b56ba54e5d4a6c1..c28ca0325b9cbc473857d912e498380d4068f5c8 100644
--- a/src/meshTools/momentOfInertia/momentOfInertia.C
+++ b/src/meshTools/momentOfInertia/momentOfInertia.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -306,7 +306,7 @@ Foam::tmp<Foam::tensorField> Foam::momentOfInertia::meshInertia
 {
     tmp<tensorField> tTf = tmp<tensorField>(new tensorField(mesh.nCells()));
 
-    tensorField& tf = tTf();
+    tensorField& tf = tTf.ref();
 
     forAll(tf, cI)
     {
diff --git a/src/meshTools/searchableSurface/searchableBox.C b/src/meshTools/searchableSurface/searchableBox.C
index a29338543651350e39d4e519561c108d048c67ff..c42150e7efe3093b1210cf0aaade27c56c2806b7 100644
--- a/src/meshTools/searchableSurface/searchableBox.C
+++ b/src/meshTools/searchableSurface/searchableBox.C
@@ -223,7 +223,7 @@ const Foam::wordList& Foam::searchableBox::regions() const
 Foam::tmp<Foam::pointField> Foam::searchableBox::coordinates() const
 {
     tmp<pointField> tCtrs = tmp<pointField>(new pointField(6));
-    pointField& ctrs = tCtrs();
+    pointField& ctrs = tCtrs.ref();
 
     const pointField pts(treeBoundBox::points());
     const faceList& fcs = treeBoundBox::faces;
diff --git a/src/meshTools/searchableSurface/searchableCylinder.C b/src/meshTools/searchableSurface/searchableCylinder.C
index 01eee6e635eda8e17b887c51b31683af4f0be78c..69568fa18d42f9e5269bdd8a9360e05052936685 100644
--- a/src/meshTools/searchableSurface/searchableCylinder.C
+++ b/src/meshTools/searchableSurface/searchableCylinder.C
@@ -30,10 +30,8 @@ License
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(searchableCylinder, 0);
-addToRunTimeSelectionTable(searchableSurface, searchableCylinder, dict);
-
+    defineTypeNameAndDebug(searchableCylinder, 0);
+    addToRunTimeSelectionTable(searchableSurface, searchableCylinder, dict);
 }
 
 
@@ -67,7 +65,7 @@ void Foam::searchableCylinder::boundingSpheres
 Foam::tmp<Foam::pointField> Foam::searchableCylinder::points() const
 {
     tmp<pointField> tPts(new pointField(2));
-    pointField& pts = tPts();
+    pointField& pts = tPts.ref();
 
     pts[0] = point1_;
     pts[1] = point2_;
diff --git a/src/meshTools/searchableSurface/searchablePlate.C b/src/meshTools/searchableSurface/searchablePlate.C
index f03b4ac26f6899a133fd79adfec1eb6689905063..166aef7752ff76dcd3aa8fa84b6bdfe54a1e1db1 100644
--- a/src/meshTools/searchableSurface/searchablePlate.C
+++ b/src/meshTools/searchableSurface/searchablePlate.C
@@ -298,7 +298,7 @@ void Foam::searchablePlate::boundingSpheres
 Foam::tmp<Foam::pointField> Foam::searchablePlate::points() const
 {
     tmp<pointField> tPts(new pointField(4));
-    pointField& pts = tPts();
+    pointField& pts = tPts.ref();
 
     pts[0] = origin_;
     pts[2] = origin_ + span_;
diff --git a/src/meshTools/searchableSurface/searchableSurfaceCollection.C b/src/meshTools/searchableSurface/searchableSurfaceCollection.C
index 868a959941060dcd06c6aeb4d7ff131d1d026e78..17477023caeb08c31b2ae4ee3d8d9c979890f69d 100644
--- a/src/meshTools/searchableSurface/searchableSurfaceCollection.C
+++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.C
@@ -325,7 +325,7 @@ Foam::tmp<Foam::pointField>
 Foam::searchableSurfaceCollection::coordinates() const
 {
     tmp<pointField> tCtrs = tmp<pointField>(new pointField(size()));
-    pointField& ctrs = tCtrs();
+    pointField& ctrs = tCtrs.ref();
 
     // Append individual coordinates
     label coordI = 0;
@@ -400,7 +400,7 @@ Foam::searchableSurfaceCollection::points() const
     }
 
     tmp<pointField> tPts(new pointField(nPoints));
-    pointField& pts = tPts();
+    pointField& pts = tPts.ref();
 
     // Append individual coordinates
     nPoints = 0;
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C
index de185b181ce5c79f71a437b64cedf0e789e7f2b9..507e839c92520d541839e35f8db4613b44113d58 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.C
@@ -382,7 +382,7 @@ void Foam::triSurfaceMesh::clearOut()
 Foam::tmp<Foam::pointField> Foam::triSurfaceMesh::coordinates() const
 {
     tmp<pointField> tPts(new pointField(8));
-    pointField& pt = tPts();
+    pointField& pt = tPts.ref();
 
     // Use copy to calculate face centres so they don't get stored
     pt = PrimitivePatch<triSurface::FaceType, SubList, const pointField&>
diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C
index 87300a99284236ce96983e538f75b4d4feba914a..e61a40391e9c6f836c056d18a01c2176fcca41c6 100644
--- a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C
+++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C
@@ -197,7 +197,7 @@ void Foam::pointToPointPlanarInterpolation::calcWeights
         (
             referenceCS_.localPosition(sourcePoints)
         );
-        vectorField& localVertices = tlocalVertices();
+        vectorField& localVertices = tlocalVertices.ref();
 
         const boundBox bb(localVertices, true);
         const point bbMid(bb.midpoint());
diff --git a/src/postProcessing/functionObjects/field/Make/options b/src/postProcessing/functionObjects/field/Make/options
index 0732822dba6309b4f07c986b07406c8eece747b1..583c538b9e62e6355c5925432239a9c0da048f38 100644
--- a/src/postProcessing/functionObjects/field/Make/options
+++ b/src/postProcessing/functionObjects/field/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
index e633d9e72ba8f6e7109a8b27999e45a6157e5b49..ca7c9976625987e7e368e29df99edc6e20109632 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C
@@ -83,7 +83,7 @@ Foam::tmp<Foam::Field<Type>> Foam::fieldValues::faceSource::getFieldValues
                 (
                     new Field<Type>(faces.size(), pTraits<Type>::zero)
                 );
-                Field<Type>& avg = tavg();
+                Field<Type>& avg = tavg.ref();
 
                 forAll(faces, faceI)
                 {
@@ -355,7 +355,7 @@ Foam::tmp<Foam::Field<Type>> Foam::fieldValues::faceSource::filterField
 ) const
 {
     tmp<Field<Type>> tvalues(new Field<Type>(faceId_.size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
 
     forAll(values, i)
     {
@@ -396,7 +396,7 @@ Foam::tmp<Foam::Field<Type>> Foam::fieldValues::faceSource::filterField
 ) const
 {
     tmp<Field<Type>> tvalues(new Field<Type>(faceId_.size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
 
     forAll(values, i)
     {
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
index 6d02c13969915605a82895b8c3770d60b8320fe3..8e6698b326854f1c9d34a99c9570615c06cd85f6 100644
--- a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
@@ -216,7 +216,7 @@ Foam::tmp<Foam::scalarField> Foam::regionSizeDistribution::divide
 )
 {
     tmp<scalarField> tresult(new scalarField(num.size()));
-    scalarField& result = tresult();
+    scalarField& result = tresult.ref();
 
     forAll(denom, i)
     {
diff --git a/src/postProcessing/functionObjects/utilities/Make/options b/src/postProcessing/functionObjects/utilities/Make/options
index ad4e27c1967426d1d3fdfdeb9b45fa67744396d6..f50a857be26bbfbdd36450ae5c271b74a7c279b2 100644
--- a/src/postProcessing/functionObjects/utilities/Make/options
+++ b/src/postProcessing/functionObjects/utilities/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/DSMC/lnInclude \
     -I$(LIB_SRC)/transportModels \
diff --git a/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.C b/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.C
index 1d9013d96e6b4aea4ba7c9e13cafe97638c9abc6..3dd1877a87c3bca1162c10fcdb27995fa423a106 100644
--- a/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.C
+++ b/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.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
@@ -148,7 +148,7 @@ Foam::tmp<Foam::volScalarField> Foam::pressureTools::pDyn
     {
         const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
 
-        tpDyn() == rho(p)*0.5*magSqr(U);
+        tpDyn.ref() == rho(p)*0.5*magSqr(U);
     }
 
     return tpDyn;
@@ -164,13 +164,13 @@ Foam::tmp<Foam::volScalarField> Foam::pressureTools::convertToCoeff
 
     if (calcCoeff_)
     {
-        tCoeff() -= dimensionedScalar("pInf", dimPressure, pInf_);
+        tCoeff.ref() -= dimensionedScalar("pInf", dimPressure, pInf_);
 
         const dimensionedScalar p0("p0", dimPressure, SMALL);
         const dimensionedVector U("U", dimVelocity, UInf_);
         const dimensionedScalar rho("rho", dimDensity, rhoInf_);
 
-        tCoeff() /= 0.5*rho*magSqr(U) + p0;
+        tCoeff.ref() /= 0.5*rho*magSqr(U) + p0;
     }
 
     return tCoeff;
diff --git a/src/randomProcesses/Make/options b/src/randomProcesses/Make/options
index 71b7873964d544eddf96d22aa40f4c3372c23c9c..287318da1ff31d9067bc255cafbdd1fa8dce4787 100644
--- a/src/randomProcesses/Make/options
+++ b/src/randomProcesses/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
diff --git a/src/randomProcesses/fft/fft.C b/src/randomProcesses/fft/fft.C
index fdc5e856ed1ceca9bba631d762ea7a98d91e9aa4..f859b37f529212cd0fe6393eaa8338736bae37e4 100644
--- a/src/randomProcesses/fft/fft.C
+++ b/src/randomProcesses/fft/fft.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
@@ -191,7 +191,7 @@ tmp<complexField> fft::forwardTransform
 {
     tmp<complexField> tfftField(new complexField(tfield));
 
-    transform(tfftField(), nn, FORWARD_TRANSFORM);
+    transform(tfftField.ref(), nn, FORWARD_TRANSFORM);
 
     tfield.clear();
 
@@ -207,7 +207,7 @@ tmp<complexField> fft::reverseTransform
 {
     tmp<complexField> tifftField(new complexField(tfield));
 
-    transform(tifftField(), nn, REVERSE_TRANSFORM);
+    transform(tifftField.ref(), nn, REVERSE_TRANSFORM);
 
     tfield.clear();
 
@@ -231,7 +231,7 @@ tmp<complexVectorField> fft::forwardTransform
 
     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
     {
-        tfftVectorField().replace
+        tfftVectorField.ref().replace
         (
             cmpt,
             forwardTransform(tfield().component(cmpt), nn)
@@ -260,7 +260,7 @@ tmp<complexVectorField> fft::reverseTransform
 
     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
     {
-        tifftVectorField().replace
+        tifftVectorField.ref().replace
         (
             cmpt,
             reverseTransform(tfield().component(cmpt), nn)
diff --git a/src/randomProcesses/noise/noiseFFT.C b/src/randomProcesses/noise/noiseFFT.C
index 88cbe3c5ba121d17f919d0bba9ed8a229d9f586c..0f09cdf817450e18461c9eb49371a82847538d40 100644
--- a/src/randomProcesses/noise/noiseFFT.C
+++ b/src/randomProcesses/noise/noiseFFT.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
@@ -140,7 +140,7 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::window
     }
 
     tmp<scalarField> tpw(new scalarField(N));
-    scalarField& pw = tpw();
+    scalarField& pw = tpw.ref();
 
     label offset = ni*windowOffset;
 
@@ -193,7 +193,7 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
             scalarField::subField(tPn2(), tPn2().size()/2)
         )
     );
-    scalarField& Pn = tPn();
+    scalarField& Pn = tPn.ref();
 
     Pn *= 2.0/sqrt(scalar(tPn2().size()));
     Pn[0] /= 2.0;
diff --git a/src/regionCoupled/Make/options b/src/regionCoupled/Make/options
index a9753c5115e98af376181d88ac517ec012762443..21a2486c48c15020941cff955e95ea0ec38f35b4 100644
--- a/src/regionCoupled/Make/options
+++ b/src/regionCoupled/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
diff --git a/src/regionCoupled/derivedFvPatchFields/energyRegionCoupled/energyRegionCoupledFvPatchScalarField.C b/src/regionCoupled/derivedFvPatchFields/energyRegionCoupled/energyRegionCoupledFvPatchScalarField.C
index 697405daf7f351e11d458922ba1d6fb408f6f66b..004fecb0e20e30d60b803ad3c38643f1b5795c4b 100644
--- a/src/regionCoupled/derivedFvPatchFields/energyRegionCoupled/energyRegionCoupledFvPatchScalarField.C
+++ b/src/regionCoupled/derivedFvPatchFields/energyRegionCoupled/energyRegionCoupledFvPatchScalarField.C
@@ -190,7 +190,7 @@ weights() const
     const scalarField nbrAlphaDelta(nbrAlpha/nbrDeltas);
 
     tmp<scalarField> tw(new scalarField(deltas.size()));
-    scalarField& w = tw();
+    scalarField& w = tw.ref();
 
     forAll(alphaDelta, faceI)
     {
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
index 501f4aa64c3048323e54de5320879863f7b7b9dd..ec261e277e2adc223713c02ee5e12b68a37d2e5f 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
@@ -196,7 +196,6 @@ void reactingOneDim::updatePhiGas()
                 }
             }
         }
-        tHsiGas().clear();
     }
 }
 
diff --git a/src/regionModels/regionModel/Make/options b/src/regionModels/regionModel/Make/options
index a3ae8da833177387e9eecf75b5e2675fc7b481f5..bcd627471b97e93547ed691c52bd70b90be9030e 100644
--- a/src/regionModels/regionModel/Make/options
+++ b/src/regionModels/regionModel/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
diff --git a/src/regionModels/regionModel/regionModel1D/regionModel1D.C b/src/regionModels/regionModel/regionModel1D/regionModel1D.C
index 63e8145c5f6105d0b332daef337582aab9c5e417..ba750c9f769d1949c1f366892b85d2ceda79d97b 100644
--- a/src/regionModels/regionModel/regionModel1D/regionModel1D.C
+++ b/src/regionModels/regionModel/regionModel1D/regionModel1D.C
@@ -183,7 +183,7 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
 )
 {
     tmp<labelField> tcellMoveMap(new labelField(regionMesh().nCells(), 0));
-    labelField& cellMoveMap = tcellMoveMap();
+    labelField& cellMoveMap = tcellMoveMap.ref();
 
     if (!moveMesh_)
     {
diff --git a/src/regionModels/surfaceFilmModels/Make/options b/src/regionModels/surfaceFilmModels/Make/options
index a0eb4828eaa8393b4c2a6c690113957d7bd4359b..0cd68b47e811caae0f788214546a2d504ba7db3a 100644
--- a/src/regionModels/surfaceFilmModels/Make/options
+++ b/src/regionModels/surfaceFilmModels/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
diff --git a/src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/nutkFilmWallFunction/nutkFilmWallFunctionFvPatchScalarField.C b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/nutkFilmWallFunction/nutkFilmWallFunctionFvPatchScalarField.C
index 2b6727b91f3df9304d541007f1bf3687c73ebc16..72b853d0de80b5ee7568d4aa7e6222a950a9f9d8 100644
--- a/src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/nutkFilmWallFunction/nutkFilmWallFunctionFvPatchScalarField.C
+++ b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/nutkFilmWallFunction/nutkFilmWallFunctionFvPatchScalarField.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
@@ -49,7 +49,7 @@ tmp<scalarField> nutkFilmWallFunctionFvPatchScalarField::calcUTau
 ) const
 {
     tmp<scalarField> tuTau(new scalarField(patch().size(), 0.0));
-    scalarField& uTau = tuTau();
+    scalarField& uTau = tuTau.ref();
 
     typedef regionModels::surfaceFilmModels::surfaceFilmModel modelType;
 
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
index e87619f67ff64002e9e140b6fdda1b0be8f56171..529b95f964f0011a1e6e0d94de995c84ca31ede2 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
@@ -299,13 +299,13 @@ tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
       + fvm::div(phi_, U_)
      ==
       - USp_
-//      - fvm::SuSp(rhoSp_, U_)
+   // - fvm::SuSp(rhoSp_, U_)
       - rhoSp_*U_
       + forces_.correct(U_)
       + turbulence_->Su(U_)
     );
 
-    fvVectorMatrix& UEqn = tUEqn();
+    fvVectorMatrix& UEqn = tUEqn.ref();
 
     UEqn.relax();
 
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H
index 3c0bd597f2bbb2912ec4614529ead94d55b0ac2b..e7647ec9e0d34be0ca2241fcce1b1c880dcb9892 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H
@@ -249,7 +249,7 @@ inline tmp<volScalarField> kinematicSingleLayer::gNormClipped() const
         )
     );
 
-    volScalarField& gNormClipped = tgNormClipped();
+    volScalarField& gNormClipped = tgNormClipped.ref();
     gNormClipped.min(0.0);
 
     return tgNormClipped;
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/constantFilmThermo/constantFilmThermo.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/constantFilmThermo/constantFilmThermo.C
index cc53371852734dfdaecf1589de158f4baa8b0a4c..3e7602868b52223b195129caa8016fb1396786c9 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/constantFilmThermo/constantFilmThermo.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/constantFilmThermo/constantFilmThermo.C
@@ -275,8 +275,8 @@ tmp<volScalarField> constantFilmThermo::rho() const
         )
     );
 
-    trho().internalField() = this->rho(0, 0);
-    trho().correctBoundaryConditions();
+    trho.ref().internalField() = this->rho(0, 0);
+    trho.ref().correctBoundaryConditions();
 
     return trho;
 }
@@ -302,8 +302,8 @@ tmp<volScalarField> constantFilmThermo::mu() const
         )
     );
 
-    tmu().internalField() = this->mu(0, 0);
-    tmu().correctBoundaryConditions();
+    tmu.ref().internalField() = this->mu(0, 0);
+    tmu.ref().correctBoundaryConditions();
 
     return tmu;
 }
@@ -329,8 +329,8 @@ tmp<volScalarField> constantFilmThermo::sigma() const
         )
     );
 
-    tsigma().internalField() = this->sigma(0, 0);
-    tsigma().correctBoundaryConditions();
+    tsigma.ref().internalField() = this->sigma(0, 0);
+    tsigma.ref().correctBoundaryConditions();
 
     return tsigma;
 }
@@ -356,8 +356,8 @@ tmp<volScalarField> constantFilmThermo::Cp() const
         )
     );
 
-    tCp().internalField() = this->Cp(0, 0);
-    tCp().correctBoundaryConditions();
+    tCp.ref().internalField() = this->Cp(0, 0);
+    tCp.ref().correctBoundaryConditions();
 
     return tCp;
 }
@@ -383,8 +383,8 @@ tmp<volScalarField> constantFilmThermo::kappa() const
         )
     );
 
-    tkappa().internalField() = this->kappa(0, 0);
-    tkappa().correctBoundaryConditions();
+    tkappa.ref().internalField() = this->kappa(0, 0);
+    tkappa.ref().correctBoundaryConditions();
 
     return tkappa;
 }
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/liquidFilmThermo/liquidFilmThermo.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/liquidFilmThermo/liquidFilmThermo.C
index 59b57136178dbfe854d88a87b87db7a8aa949091..13f6937b2a42ab989547ba3900215baed7df4171 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/liquidFilmThermo/liquidFilmThermo.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmThermoModel/liquidFilmThermo/liquidFilmThermo.C
@@ -255,7 +255,7 @@ tmp<volScalarField> liquidFilmThermo::rho() const
         )
     );
 
-    scalarField& rho = trho().internalField();
+    scalarField& rho = trho.ref().internalField();
 
     if (useReferenceValues_)
     {
@@ -277,7 +277,7 @@ tmp<volScalarField> liquidFilmThermo::rho() const
         }
     }
 
-    trho().correctBoundaryConditions();
+    trho.ref().correctBoundaryConditions();
 
     return trho;
 }
@@ -303,7 +303,7 @@ tmp<volScalarField> liquidFilmThermo::mu() const
         )
     );
 
-    scalarField& mu = tmu().internalField();
+    scalarField& mu = tmu.ref().internalField();
 
     if (useReferenceValues_)
     {
@@ -325,7 +325,7 @@ tmp<volScalarField> liquidFilmThermo::mu() const
         }
     }
 
-    tmu().correctBoundaryConditions();
+    tmu.ref().correctBoundaryConditions();
 
     return tmu;
 }
@@ -351,7 +351,7 @@ tmp<volScalarField> liquidFilmThermo::sigma() const
         )
     );
 
-    scalarField& sigma = tsigma().internalField();
+    scalarField& sigma = tsigma.ref().internalField();
 
     if (useReferenceValues_)
     {
@@ -373,7 +373,7 @@ tmp<volScalarField> liquidFilmThermo::sigma() const
         }
     }
 
-    tsigma().correctBoundaryConditions();
+    tsigma.ref().correctBoundaryConditions();
 
     return tsigma;
 }
@@ -399,7 +399,7 @@ tmp<volScalarField> liquidFilmThermo::Cp() const
         )
     );
 
-    scalarField& Cp = tCp().internalField();
+    scalarField& Cp = tCp.ref().internalField();
 
     if (useReferenceValues_)
     {
@@ -421,7 +421,7 @@ tmp<volScalarField> liquidFilmThermo::Cp() const
         }
     }
 
-    tCp().correctBoundaryConditions();
+    tCp.ref().correctBoundaryConditions();
 
     return tCp;
 }
@@ -447,7 +447,7 @@ tmp<volScalarField> liquidFilmThermo::kappa() const
         )
     );
 
-    scalarField& kappa = tkappa().internalField();
+    scalarField& kappa = tkappa.ref().internalField();
 
     if (useReferenceValues_)
     {
@@ -469,7 +469,7 @@ tmp<volScalarField> liquidFilmThermo::kappa() const
         }
     }
 
-    tkappa().correctBoundaryConditions();
+    tkappa.ref().correctBoundaryConditions();
 
     return tkappa;
 }
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C
index df7ccf0845f49f940def6db4ba2b2669997015aa..4a1c6383e8ab816d95c64ebf41478cefcc01fbf8 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C
@@ -89,8 +89,8 @@ tmp<volVectorField> laminar::Us() const
     );
 
     // apply quadratic profile
-    tUs() = Foam::sqrt(2.0)*owner_.U();
-    tUs().correctBoundaryConditions();
+    tUs.ref() = Foam::sqrt(2.0)*owner_.U();
+    tUs.ref().correctBoundaryConditions();
 
     return tUs;
 }
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C
index 756b1d8cc40aeaf7e94f08fcdf60c7709b4133f9..c8d111579e7ea190f07fe28e0d4e949042752c58 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.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
@@ -153,7 +153,7 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
         )
     );
 
-    vectorField& force = tForce().internalField();
+    vectorField& force = tForce.ref().internalField();
 
     const labelUList& own = owner_.regionMesh().owner();
     const labelUList& nbr = owner_.regionMesh().neighbour();
@@ -229,7 +229,7 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
     tmp<fvVectorMatrix>
         tfvm(new fvVectorMatrix(U, dimForce/dimArea*dimVolume));
 
-    tfvm() += tForce;
+    tfvm.ref() += tForce;
 
     return tfvm;
 }
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/forceList/forceList.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/forceList/forceList.C
index a6b54fe235309c85bfaa8f04674e46cef907bd4b..56a1454c46b35e48fdb6c2fc09a80d4ad6e4cddc 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/forceList/forceList.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/forceList/forceList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,7 +81,7 @@ tmp<fvVectorMatrix> forceList::correct(volVectorField& U)
     (
         new fvVectorMatrix(U, dimForce/dimArea*dimVolume)
     );
-    fvVectorMatrix& result = tResult();
+    fvVectorMatrix& result = tResult.ref();
 
     forAll(*this, i)
     {
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/thermocapillaryForce/thermocapillaryForce.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/thermocapillaryForce/thermocapillaryForce.C
index 0d681da0080f0cbea67a5d7b842cfc04079ac0dd..19614cbcf91b5fe3adc3f07876a8b4ab607dc75f 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/thermocapillaryForce/thermocapillaryForce.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/thermocapillaryForce/thermocapillaryForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,7 +68,7 @@ tmp<fvVectorMatrix> thermocapillaryForce::correct(volVectorField& U)
     tmp<fvVectorMatrix>
         tfvm(new fvVectorMatrix(U, dimForce/dimArea*dimVolume));
 
-    tfvm() += fvc::grad(sigma);
+    tfvm.ref() += fvc::grad(sigma);
 
     return tfvm;
 }
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
index 2709cd08a58320488ad80d403eabab8ffbfc5056..72ec8a44ba9fe0d19a28c3d154321b61021af214 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
@@ -78,7 +78,7 @@ tmp<volScalarField> curvatureSeparation::calcInvR1
     );
 
 
-    scalarField& invR1 = tinvR1().internalField();
+    scalarField& invR1 = tinvR1.ref().internalField();
 
     // apply defined patch radii
     const scalar rMin = 1e-6;
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C
index 4f5c60f1455945f473d372f9d48a7ff1be03b6b2..489b453cb90df3022efb1fb398b940638d5774a4 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C
@@ -124,7 +124,7 @@ tmp<volScalarField> constantRadiation::Shs()
 
     if ((time >= timeStart_) && (time <= timeStart_ + duration_))
     {
-        scalarField& Shs = tShs();
+        scalarField& Shs = tShs.ref();
         const scalarField& Qr = QrConst_.internalField();
         const scalarField& alpha = owner_.alpha().internalField();
 
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C
index b653ee4638b76422400e12315ec3fb9dcc6b12fc..a6350c3ac18131e368eafd92be0a39caa5d3fb68 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C
@@ -107,7 +107,7 @@ tmp<volScalarField> primaryRadiation::Shs()
         )
     );
 
-    scalarField& Shs = tShs();
+    scalarField& Shs = tShs.ref();
     const scalarField& QinP = QinPrimary_.internalField();
     const scalarField& alpha = owner_.alpha().internalField();
 
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
index 3505f9c8a02851f02aae4e63f9d3a9126f8263c9..2d6e8b77b8f484ac4d3884f9a0fb6a909ca27624 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
@@ -124,7 +124,7 @@ tmp<volScalarField> standardRadiation::Shs()
         )
     );
 
-    scalarField& Shs = tShs();
+    scalarField& Shs = tShs.ref();
     const scalarField& QinP = QinPrimary_.internalField();
     const scalarField& delta = owner_.delta().internalField();
     const scalarField& alpha = owner_.alpha().internalField();
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
index f29d5d64ed8b28c04ef3b688f778f402956f9e7d..c4fc973d23cec6b4fcdacc51bed4fe1491143942 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
@@ -751,7 +751,7 @@ tmp<DimensionedField<scalar, volMesh>> thermoSingleLayer::Srho() const
         )
     );
 
-    scalarField& Srho = tSrho();
+    scalarField& Srho = tSrho.ref();
     const scalarField& V = primaryMesh().V();
     const scalar dt = time_.deltaTValue();
 
@@ -805,7 +805,7 @@ tmp<DimensionedField<scalar, volMesh>> thermoSingleLayer::Srho
 
     if (vapId == i)
     {
-        scalarField& Srho = tSrho();
+        scalarField& Srho = tSrho.ref();
         const scalarField& V = primaryMesh().V();
         const scalar dt = time().deltaTValue();
 
@@ -855,7 +855,7 @@ tmp<DimensionedField<scalar, volMesh>> thermoSingleLayer::Sh() const
 /*
     phase change energy fed back into the film...
 
-    scalarField& Sh = tSh();
+    scalarField& Sh = tSh.ref();
     const scalarField& V = primaryMesh().V();
     const scalar dt = time_.deltaTValue();
 
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
index 223254c38e60e4fa95b2df7d6ced18cc138df070..22bd215deb2574b2469cd90e61b730e5264be913 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
@@ -100,8 +100,8 @@ inline tmp<volScalarField> thermoSingleLayer::T
         )
     );
 
-    tT().min(Tmax_);
-    tT().max(Tmin_);
+    tT.ref().min(Tmax_);
+    tT.ref().max(Tmin_);
 
     return tT;
 }
diff --git a/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.C b/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.C
index 360bdfeb1fc246c0c7b8e319730360fcb1721994..1490e558b1ec0af673b143f9926d46f5bdc1fd56 100644
--- a/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.C
+++ b/src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.C
@@ -91,7 +91,7 @@ void thermalBaffle::solveEnergy()
         )
     );
 
-    volScalarField& Q = tQ();
+    volScalarField& Q = tQ.ref();
 
     volScalarField rho("rho", thermo_->rho());
     volScalarField alpha("alpha", thermo_->alpha());
diff --git a/src/sampling/Make/options b/src/sampling/Make/options
index 23e6bc81a7e708d69b5fd14d27073827646e14ef..15ca6b932f8ba3574029ca7ec291cad669d9b7a3 100644
--- a/src/sampling/Make/options
+++ b/src/sampling/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
diff --git a/src/sampling/meshToMesh/meshToMeshTemplates.C b/src/sampling/meshToMesh/meshToMeshTemplates.C
index 24d01d89a3bb05d4cd031ab0cac628aa0779189c..958b932be59bde13d6a0e2c24c7549084a61323a 100644
--- a/src/sampling/meshToMesh/meshToMeshTemplates.C
+++ b/src/sampling/meshToMesh/meshToMeshTemplates.C
@@ -161,7 +161,7 @@ Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapSrcToTgt
         )
     );
 
-    mapSrcToTgt(srcField, cop, tresult());
+    mapSrcToTgt(srcField, cop, tresult.ref());
 
     return tresult;
 }
@@ -280,7 +280,7 @@ Foam::tmp<Foam::Field<Type>> Foam::meshToMesh::mapTgtToSrc
         )
     );
 
-    mapTgtToSrc(tgtField, cop, tresult());
+    mapTgtToSrc(tgtField, cop, tresult.ref());
 
     return tresult;
 }
@@ -466,7 +466,7 @@ Foam::meshToMesh::mapSrcToTgt
         )
     );
 
-    mapSrcToTgt(field, cop, tresult());
+    mapSrcToTgt(field, cop, tresult.ref());
 
     return tresult;
 }
@@ -655,7 +655,7 @@ Foam::meshToMesh::mapTgtToSrc
         )
     );
 
-    mapTgtToSrc(field, cop, tresult());
+    mapTgtToSrc(field, cop, tresult.ref());
 
     return tresult;
 }
diff --git a/src/sampling/probes/patchProbesTemplates.C b/src/sampling/probes/patchProbesTemplates.C
index 325ff16a0ef25a1aea41281b39aa10b035c7c87f..59dfd0f50a5a78cf89a493d5de85cbdd60d06f2c 100644
--- a/src/sampling/probes/patchProbesTemplates.C
+++ b/src/sampling/probes/patchProbesTemplates.C
@@ -202,7 +202,7 @@ Foam::patchProbes::sample
         new Field<Type>(this->size(), unsetVal)
     );
 
-    Field<Type>& values = tValues();
+    Field<Type>& values = tValues.ref();
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
@@ -253,7 +253,7 @@ Foam::patchProbes::sample
         new Field<Type>(this->size(), unsetVal)
     );
 
-    Field<Type>& values = tValues();
+    Field<Type>& values = tValues.ref();
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
diff --git a/src/sampling/probes/probesTemplates.C b/src/sampling/probes/probesTemplates.C
index f035d5fe03a0cb123b81182464bf9f8e4a668731..37d0bd98d56e7a8ee97743896ae78e75d90431a5 100644
--- a/src/sampling/probes/probesTemplates.C
+++ b/src/sampling/probes/probesTemplates.C
@@ -225,7 +225,7 @@ Foam::probes::sample
         new Field<Type>(this->size(), unsetVal)
     );
 
-    Field<Type>& values = tValues();
+    Field<Type>& values = tValues.ref();
 
     if (fixedLocations_)
     {
@@ -295,7 +295,7 @@ Foam::probes::sample
         new Field<Type>(this->size(), unsetVal)
     );
 
-    Field<Type>& values = tValues();
+    Field<Type>& values = tValues.ref();
 
     forAll(*this, probeI)
     {
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
index ac276779242e2f2714dfda783639cfd0b05cecea..1eba2063c40cd243b69be80afc329c503aa8a35a 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
@@ -596,7 +596,7 @@ Foam::isoSurfaceCell::interpolate
 
     // One value per point
     tmp<Field<Type>> tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
 
     forAll(triPoints, i)
     {
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index 58950e006bcf559aa0901f291b6438a22e0231b9..ea92b65ef74ff5f6bde42292c18947a0adfae9e5 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -70,7 +70,7 @@ Foam::isoSurface::adaptPatchFields
             true        // preserveCouples
         )
     );
-    FieldType& sliceFld = tsliceFld();
+    FieldType& sliceFld = tsliceFld.ref();
 
     const fvMesh& mesh = fld.mesh();
 
@@ -737,7 +737,7 @@ Foam::isoSurface::interpolate
     (
         new Field<Type>(points().size(), pTraits<Type>::zero)
     );
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
     labelList nValues(values.size(), 0);
 
     forAll(triPoints, i)
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
index a80bd81a3ba5bd05e13e3b0cc06ce169f5f43e98..72bcbbf701bcdc11770a3755ab1649cc4350d6f3 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
@@ -57,7 +57,7 @@ Foam::sampledIsoSurfaceCell::interpolateField
 
     // One value per point
     tmp<Field<Type>> tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
 
     boolList pointDone(points().size(), false);
 
diff --git a/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C b/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C
index dc7f0cbb6e1eaa81c6dc6b2ef6e235e416ff6ce0..58bd02fcb0fa29f9608038cc1722c69191ecfe78 100644
--- a/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C
+++ b/src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C
@@ -36,7 +36,7 @@ Foam::sampledPatch::sampleField
 {
     // One value per face
     tmp<Field<Type>> tvalues(new Field<Type>(patchFaceLabels_.size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
     forAll(patchFaceLabels_, i)
     {
         label patchI = patchIDs_[patchIndex_[i]];
@@ -57,7 +57,7 @@ Foam::sampledPatch::sampleField
 {
     // One value per face
     tmp<Field<Type>> tvalues(new Field<Type>(patchFaceLabels_.size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
 
     forAll(patchFaceLabels_, i)
     {
@@ -78,7 +78,7 @@ Foam::sampledPatch::interpolateField
 {
     // One value per vertex
     tmp<Field<Type>> tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
 
     const labelList& own = mesh().faceOwner();
 
diff --git a/src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalFieldTemplates.C b/src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalFieldTemplates.C
index 56852b4ec2a1274dfe4961f4f5049fa44df56f2b..be2348729f6a95419d6844d59a4af7b49bad5b03 100644
--- a/src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalFieldTemplates.C
+++ b/src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalFieldTemplates.C
@@ -38,7 +38,7 @@ Foam::sampledPatchInternalField::sampleField
 {
     // One value per face
     tmp<Field<Type>> tvalues(new Field<Type>(patchFaceLabels().size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
 
     forAll(patchStart(), i)
     {
diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C b/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C
index 52e830d75fab7d70845ed1a1e2fa7e02a236a758..60e6b5a489b1f6647c8c6e0d993b1f097afff97f 100644
--- a/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C
+++ b/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C
@@ -47,7 +47,7 @@ Foam::sampledPlane::interpolateField
 {
     // One value per point
     tmp<Field<Type>> tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
 
     boolList pointDone(points().size(), false);
 
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurface.C b/src/sampling/sampledSurface/sampledSurface/sampledSurface.C
index 59a34de598076e0fc0a0cc89fbe39ad47546aca8..38c478edc09ef183409670e0e959b1f54d84256d 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurface.C
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurface.C
@@ -291,7 +291,7 @@ Foam::tmp<Foam::Field<Foam::scalar>>
 Foam::sampledSurface::project(const Field<scalar>& field) const
 {
     tmp<Field<scalar>> tRes(new Field<scalar>(faces().size()));
-    Field<scalar>& res = tRes();
+    Field<scalar>& res = tRes.ref();
 
     forAll(faces(), faceI)
     {
@@ -306,7 +306,7 @@ Foam::tmp<Foam::Field<Foam::scalar>>
 Foam::sampledSurface::project(const Field<vector>& field) const
 {
     tmp<Field<scalar>> tRes(new Field<scalar>(faces().size()));
-    project(tRes(), field);
+    project(tRes.ref(), field);
     return tRes;
 }
 
@@ -315,7 +315,7 @@ Foam::tmp<Foam::Field<Foam::vector>>
 Foam::sampledSurface::project(const Field<sphericalTensor>& field) const
 {
     tmp<Field<vector>> tRes(new Field<vector>(faces().size()));
-    project(tRes(), field);
+    project(tRes.ref(), field);
     return tRes;
 }
 
@@ -324,7 +324,7 @@ Foam::tmp<Foam::Field<Foam::vector>>
 Foam::sampledSurface::project(const Field<symmTensor>& field) const
 {
     tmp<Field<vector>> tRes(new Field<vector>(faces().size()));
-    project(tRes(), field);
+    project(tRes.ref(), field);
     return tRes;
 }
 
@@ -333,7 +333,7 @@ Foam::tmp<Foam::Field<Foam::vector>>
 Foam::sampledSurface::project(const Field<tensor>& field) const
 {
     tmp<Field<vector>> tRes(new Field<vector>(faces().size()));
-    project(tRes(), field);
+    project(tRes.ref(), field);
     return tRes;
 }
 
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C b/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
index d361cdd0294d6a89315d5add628205b446239b14..dcf176b56a5e90663e17912165ff51392e5b75c5 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C
@@ -177,7 +177,7 @@ Foam::sampledSurface::pointAverage
             dimensioned<Type>("zero", dimless, pTraits<Type>::zero)
         )
     );
-    GeometricField<Type, fvPatchField, volMesh>& cellAvg = tcellAvg();
+    GeometricField<Type, fvPatchField, volMesh>& cellAvg = tcellAvg.ref();
 
     labelField nPointCells(mesh.nCells(), 0);
     {
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
index 40190c9ae36a32deb86bf29e2a4bc884f60a3836..c286dde965cf82fd4e3b34fc8ea77cdac2afc07a 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
@@ -36,7 +36,7 @@ Foam::sampledTriSurfaceMesh::sampleField
 {
     // One value per face
     tmp<Field<Type>> tvalues(new Field<Type>(sampleElements_.size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
 
     if (sampleSource_ == cells || sampleSource_ == insideCells)
     {
@@ -92,7 +92,7 @@ Foam::sampledTriSurfaceMesh::interpolateField
 {
     // One value per vertex
     tmp<Field<Type>> tvalues(new Field<Type>(sampleElements_.size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
 
     if (sampleSource_ == cells || sampleSource_ == insideCells)
     {
diff --git a/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFacesTemplates.C b/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFacesTemplates.C
index d78da368378846133a7aa8458b7953478f09dbd1..03c2fffb4adeb32ae980338cd2e2b235c19349ab 100644
--- a/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFacesTemplates.C
+++ b/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFacesTemplates.C
@@ -58,7 +58,7 @@ Foam::sampledThresholdCellFaces::interpolateField
 
     // One value per point
     tmp<Field<Type>> tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
+    Field<Type>& values = tvalues.ref();
 
     boolList pointDone(points().size(), false);
 
diff --git a/src/sixDoFRigidBodyMotion/Make/options b/src/sixDoFRigidBodyMotion/Make/options
index f2367d2e1d47e8b701978c06a03ad5696c8c8668..fd1f401690d3eb816fb956ec3114c8b18d0d4e00 100644
--- a/src/sixDoFRigidBodyMotion/Make/options
+++ b/src/sixDoFRigidBodyMotion/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/postProcessing/functionObjects/forces/lnInclude \
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
index e7a8b92d339880ba8820abd4946cd0bd8abf96f6..c3a94123d1eb85b593ce6cc0fd2a3a086c1a4989 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.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
@@ -356,7 +356,7 @@ Foam::tmp<Foam::pointField> Foam::sixDoFRigidBodyMotion::transform
     );
 
     tmp<pointField> tpoints(new pointField(initialPoints));
-    pointField& points = tpoints();
+    pointField& points = tpoints.ref();
 
     forAll(points, pointi)
     {
diff --git a/src/surfMesh/Make/options b/src/surfMesh/Make/options
index 7e207d0dbeaca258e5a72af8b4eb7bacefc0dee8..079607db058a7cc9acdcca3a7aa9b0f20a0a36a0 100644
--- a/src/surfMesh/Make/options
+++ b/src/surfMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/fileFormats/lnInclude
 
 LIB_LIBS = \
diff --git a/src/thermophysicalModels/SLGThermo/Make/options b/src/thermophysicalModels/SLGThermo/Make/options
index c4f0b2a8a95fa2420cb555ebad3e0f96387baca4..3dac244d79650919f883a3d0123751ab278f633a 100644
--- a/src/thermophysicalModels/SLGThermo/Make/options
+++ b/src/thermophysicalModels/SLGThermo/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 71b7873964d544eddf96d22aa40f4c3372c23c9c..287318da1ff31d9067bc255cafbdd1fa8dce4787 100644
--- a/src/thermophysicalModels/barotropicCompressibilityModel/Make/options
+++ b/src/thermophysicalModels/barotropicCompressibilityModel/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
diff --git a/src/thermophysicalModels/basic/Make/options b/src/thermophysicalModels/basic/Make/options
index b5c859baf1457da06327d60192105d4b8fbaad53..95f178b455134e96eb3b0c998c3bab5429182cdc 100644
--- a/src/thermophysicalModels/basic/Make/options
+++ b/src/thermophysicalModels/basic/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/src/thermophysicalModels/basic/heThermo/heThermo.C b/src/thermophysicalModels/basic/heThermo/heThermo.C
index 0cc1354ffa30e7999a6581b680bed8eff9c61e3e..1db5c4ac7ca00ef6ebcad689b4ab7b00155dd3b3 100644
--- a/src/thermophysicalModels/basic/heThermo/heThermo.C
+++ b/src/thermophysicalModels/basic/heThermo/heThermo.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
@@ -184,7 +184,7 @@ Foam::tmp<Foam::volScalarField> Foam::heThermo<BasicThermo, MixtureType>::he
         )
     );
 
-    volScalarField& he = the();
+    volScalarField& he = the.ref();
     scalarField& heCells = he.internalField();
     const scalarField& pCells = p.internalField();
     const scalarField& TCells = T.internalField();
@@ -221,7 +221,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
 ) const
 {
     tmp<scalarField> the(new scalarField(T.size()));
-    scalarField& he = the();
+    scalarField& he = the.ref();
 
     forAll(T, celli)
     {
@@ -241,7 +241,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::he
 ) const
 {
     tmp<scalarField> the(new scalarField(T.size()));
-    scalarField& he = the();
+    scalarField& he = the.ref();
 
     forAll(T, facei)
     {
@@ -277,7 +277,7 @@ Foam::heThermo<BasicThermo, MixtureType>::hc() const
         )
     );
 
-    volScalarField& hcf = thc();
+    volScalarField& hcf = thc.ref();
     scalarField& hcCells = hcf.internalField();
 
     forAll(hcCells, celli)
@@ -308,7 +308,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cp
 ) const
 {
     tmp<scalarField> tCp(new scalarField(T.size()));
-    scalarField& cp = tCp();
+    scalarField& cp = tCp.ref();
 
     forAll(T, facei)
     {
@@ -344,7 +344,7 @@ Foam::heThermo<BasicThermo, MixtureType>::Cp() const
         )
     );
 
-    volScalarField& cp = tCp();
+    volScalarField& cp = tCp.ref();
 
     forAll(this->T_, celli)
     {
@@ -379,7 +379,7 @@ Foam::heThermo<BasicThermo, MixtureType>::Cv
 ) const
 {
     tmp<scalarField> tCv(new scalarField(T.size()));
-    scalarField& cv = tCv();
+    scalarField& cv = tCv.ref();
 
     forAll(T, facei)
     {
@@ -415,7 +415,7 @@ Foam::heThermo<BasicThermo, MixtureType>::Cv() const
         )
     );
 
-    volScalarField& cv = tCv();
+    volScalarField& cv = tCv.ref();
 
     forAll(this->T_, celli)
     {
@@ -446,7 +446,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::gamma
 ) const
 {
     tmp<scalarField> tgamma(new scalarField(T.size()));
-    scalarField& cpv = tgamma();
+    scalarField& cpv = tgamma.ref();
 
     forAll(T, facei)
     {
@@ -482,7 +482,7 @@ Foam::heThermo<BasicThermo, MixtureType>::gamma() const
         )
     );
 
-    volScalarField& cpv = tgamma();
+    volScalarField& cpv = tgamma.ref();
 
     forAll(this->T_, celli)
     {
@@ -519,7 +519,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cpv
 ) const
 {
     tmp<scalarField> tCpv(new scalarField(T.size()));
-    scalarField& cpv = tCpv();
+    scalarField& cpv = tCpv.ref();
 
     forAll(T, facei)
     {
@@ -555,7 +555,7 @@ Foam::heThermo<BasicThermo, MixtureType>::Cpv() const
         )
     );
 
-    volScalarField& cpv = tCpv();
+    volScalarField& cpv = tCpv.ref();
 
     forAll(this->T_, celli)
     {
@@ -589,7 +589,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::CpByCpv
 ) const
 {
     tmp<scalarField> tCpByCpv(new scalarField(T.size()));
-    scalarField& cpByCpv = tCpByCpv();
+    scalarField& cpByCpv = tCpByCpv.ref();
 
     forAll(T, facei)
     {
@@ -625,7 +625,7 @@ Foam::heThermo<BasicThermo, MixtureType>::CpByCpv() const
         )
     );
 
-    volScalarField& cpByCpv = tCpByCpv();
+    volScalarField& cpByCpv = tCpByCpv.ref();
 
     forAll(this->T_, celli)
     {
@@ -666,7 +666,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::THE
 ) const
 {
     tmp<scalarField> tT(new scalarField(h.size()));
-    scalarField& T = tT();
+    scalarField& T = tT.ref();
 
     forAll(h, celli)
     {
@@ -689,7 +689,7 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::THE
 {
 
     tmp<scalarField> tT(new scalarField(h.size()));
-    scalarField& T = tT();
+    scalarField& T = tT.ref();
     forAll(h, facei)
     {
         T[facei] = this->patchFaceMixture
@@ -708,7 +708,7 @@ Foam::tmp<Foam::volScalarField>
 Foam::heThermo<BasicThermo, MixtureType>::kappa() const
 {
     tmp<Foam::volScalarField> kappa(Cp()*this->alpha_);
-    kappa().rename("kappa");
+    kappa.ref().rename("kappa");
     return kappa;
 }
 
@@ -737,7 +737,7 @@ Foam::heThermo<BasicThermo, MixtureType>::kappaEff
 ) const
 {
     tmp<Foam::volScalarField> kappaEff(Cp()*(this->alpha_ + alphat));
-    kappaEff().rename("kappaEff");
+    kappaEff.ref().rename("kappaEff");
     return kappaEff;
 }
 
@@ -772,7 +772,7 @@ Foam::heThermo<BasicThermo, MixtureType>::alphaEff
 ) const
 {
     tmp<Foam::volScalarField> alphaEff(this->CpByCpv()*(this->alpha_ + alphat));
-    alphaEff().rename("alphaEff");
+    alphaEff.ref().rename("alphaEff");
     return alphaEff;
 }
 
diff --git a/src/thermophysicalModels/chemistryModel/Make/options b/src/thermophysicalModels/chemistryModel/Make/options
index de52d7e6ae8750b7e48f9e7dc740e3c2da538042..7a1e7f43d2b62cfdf78057946e1634fa19f1b922 100644
--- a/src/thermophysicalModels/chemistryModel/Make/options
+++ b/src/thermophysicalModels/chemistryModel/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
index ee7a8974d61d37218dd39009b8e5c8cd72b43260..7c2d9beac2c270b91ec1fd5f5cf803cfa077cfdf 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
@@ -104,7 +104,7 @@ Foam::chemistryModel<CompType, ThermoType>::omega
     label lRef, rRef;
 
     tmp<scalarField> tom(new scalarField(nEqns(), 0.0));
-    scalarField& om = tom();
+    scalarField& om = tom.ref();
 
     forAll(reactions_, i)
     {
@@ -500,7 +500,7 @@ Foam::chemistryModel<CompType, ThermoType>::tc() const
         )
     );
 
-    scalarField& tc = ttc();
+    scalarField& tc = ttc.ref();
     const scalarField& T = this->thermo().T();
     const scalarField& p = this->thermo().p();
 
@@ -540,7 +540,7 @@ Foam::chemistryModel<CompType, ThermoType>::tc() const
     }
 
 
-    ttc().correctBoundaryConditions();
+    ttc.ref().correctBoundaryConditions();
 
     return ttc;
 }
@@ -570,7 +570,7 @@ Foam::chemistryModel<CompType, ThermoType>::Sh() const
 
     if (this->chemistry_)
     {
-        scalarField& Sh = tSh();
+        scalarField& Sh = tSh.ref();
 
         forAll(Y_, i)
         {
@@ -610,7 +610,7 @@ Foam::chemistryModel<CompType, ThermoType>::dQ() const
 
     if (this->chemistry_)
     {
-        volScalarField& dQ = tdQ();
+        volScalarField& dQ = tdQ.ref();
         dQ.dimensionedInternalField() = this->mesh_.V()*Sh()();
     }
 
@@ -668,7 +668,7 @@ Foam::chemistryModel<CompType, ThermoType>::calculateRR
         )
     );
 
-    DimensionedField<scalar, volMesh>& RR = tRR();
+    DimensionedField<scalar, volMesh>& RR = tRR.ref();
 
     const scalarField& T = this->thermo().T();
     const scalarField& p = this->thermo().p();
diff --git a/src/thermophysicalModels/laminarFlameSpeed/Gulders/Gulders.C b/src/thermophysicalModels/laminarFlameSpeed/Gulders/Gulders.C
index a1563921c2bd340d46323c9893f4311f2e606a0e..8ed9f358537b1f8a56afa092781ab26842b5df3a 100644
--- a/src/thermophysicalModels/laminarFlameSpeed/Gulders/Gulders.C
+++ b/src/thermophysicalModels/laminarFlameSpeed/Gulders/Gulders.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
@@ -128,7 +128,7 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::Gulders::Su0pTphi
         )
     );
 
-    volScalarField& Su0 = tSu0();
+    volScalarField& Su0 = tSu0.ref();
 
     forAll(Su0, celli)
     {
@@ -179,7 +179,7 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::Gulders::Su0pTphi
         )
     );
 
-    volScalarField& Su0 = tSu0();
+    volScalarField& Su0 = tSu0.ref();
 
     forAll(Su0, celli)
     {
diff --git a/src/thermophysicalModels/laminarFlameSpeed/GuldersEGR/GuldersEGR.C b/src/thermophysicalModels/laminarFlameSpeed/GuldersEGR/GuldersEGR.C
index 173e3d8841515f2f36c23836fb409f97a8482bde..74e62e71143010b26b21bfd40edc51b7f18c219d 100644
--- a/src/thermophysicalModels/laminarFlameSpeed/GuldersEGR/GuldersEGR.C
+++ b/src/thermophysicalModels/laminarFlameSpeed/GuldersEGR/GuldersEGR.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
@@ -128,7 +128,7 @@ Foam::laminarFlameSpeedModels::GuldersEGR::Su0pTphi
         )
     );
 
-    volScalarField& Su0 = tSu0();
+    volScalarField& Su0 = tSu0.ref();
 
     forAll(Su0, celli)
     {
@@ -181,7 +181,7 @@ Foam::laminarFlameSpeedModels::GuldersEGR::Su0pTphi
         )
     );
 
-    volScalarField& Su0 = tSu0();
+    volScalarField& Su0 = tSu0.ref();
 
     forAll(Su0, celli)
     {
diff --git a/src/thermophysicalModels/laminarFlameSpeed/Make/options b/src/thermophysicalModels/laminarFlameSpeed/Make/options
index f97edb8458c3d408db331fb5aa10930eadd2fc73..925d4a58aab37f576d449679c9073303acf44228 100644
--- a/src/thermophysicalModels/laminarFlameSpeed/Make/options
+++ b/src/thermophysicalModels/laminarFlameSpeed/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
diff --git a/src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.C b/src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.C
index e57be1b68f90511f86b64def90c3715b183cf944..c3eaf9c0e908baf207a2536274d35710d6fc1402 100644
--- a/src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.C
+++ b/src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.C
@@ -343,7 +343,7 @@ Foam::laminarFlameSpeedModels::RaviPetersen::operator()() const
         )
     );
 
-    volScalarField& Su0 = tSu0();
+    volScalarField& Su0 = tSu0.ref();
 
     forAll(Su0, cellI)
     {
diff --git a/src/thermophysicalModels/properties/liquidMixtureProperties/Make/options b/src/thermophysicalModels/properties/liquidMixtureProperties/Make/options
index 138f439935dc242e7002b938b31644854d43b28b..690761730fada84af3821927c633c9bdd8f77e9b 100644
--- a/src/thermophysicalModels/properties/liquidMixtureProperties/Make/options
+++ b/src/thermophysicalModels/properties/liquidMixtureProperties/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
@@ -8,4 +8,3 @@ EXE_INC = \
 LIB_LIBS = \
     -lliquidProperties \
     -lthermophysicalFunctions
-
diff --git a/src/thermophysicalModels/properties/liquidProperties/Make/options b/src/thermophysicalModels/properties/liquidProperties/Make/options
index b964b61294c787fe650b03faee1b09bbb9b48256..e1fbd3ed45354782b62c37ba685a67beea951632 100644
--- a/src/thermophysicalModels/properties/liquidProperties/Make/options
+++ b/src/thermophysicalModels/properties/liquidProperties/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 f2d8f809d5d12dac81a95c25e5a85a21cc9adac1..97cc50921233c7a9d0cee3aab90633b64e1b17e6 100644
--- a/src/thermophysicalModels/properties/solidMixtureProperties/Make/options
+++ b/src/thermophysicalModels/properties/solidMixtureProperties/Make/options
@@ -1,3 +1,3 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 848cce789f2aad33d9ed6eda979ea31c32377425..f8779f538f675e83239fbc61890c5653da29ec84 100644
--- a/src/thermophysicalModels/properties/solidProperties/Make/options
+++ b/src/thermophysicalModels/properties/solidProperties/Make/options
@@ -1,2 +1,2 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
diff --git a/src/thermophysicalModels/radiation/Make/options b/src/thermophysicalModels/radiation/Make/options
index 726b76e76442077c0a0c4f77334e54d9bfe6a6ad..ea189b0b63517b6784f7380df21d424531b5b4b9 100644
--- a/src/thermophysicalModels/radiation/Make/options
+++ b/src/thermophysicalModels/radiation/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C
index 4b07238462100d7a127606b944bbd04e3f15421a..b2fdb85d1632b48b63af64a36ac1dff199fefc2e 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/blackBodyEmission/blackBodyEmission.C
@@ -227,18 +227,15 @@ Foam::radiation::blackBodyEmission::EbDeltaLambdaT
     }
     else
     {
+        scalarField& Ebif = Eb.ref();
+
         forAll(T, i)
         {
             scalar T1 = fLambdaT(band[1]*T[i]);
             scalar T2 = fLambdaT(band[0]*T[i]);
-            dimensionedScalar fLambdaDelta
-            (
-                "fLambdaDelta",
-                dimless,
-                T1 - T2
-            );
-            Eb()[i] = Eb()[i]*fLambdaDelta.value();
+            Ebif[i] *= T1 - T2;
         }
+
         return Eb;
     }
 }
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
index c27520f4e581d8b2bf0c8254b702dc2a9fe111b4..92f0ae3084bf72e977f2f9aef29de83f4f9ddc01 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.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
@@ -253,11 +253,11 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
             );
         }
 
-        IiEq().relax();
+        IiEq.ref().relax();
 
         const solverPerformance ILambdaSol = solve
         (
-            IiEq(),
+            IiEq.ref(),
             mesh_.solver("Ii")
         );
 
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
index 208425212fd696a58855151dbc772dd076e9ae56..5d7e9070ef02a48010648a7b830863e29e714088 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
@@ -213,7 +213,7 @@ Foam::radiation::greyMeanAbsorptionEmission::aCont(const label bandI) const
         )
     );
 
-    scalarField& a = ta().internalField();
+    scalarField& a = ta.ref().internalField();
 
     forAll(a, cellI)
     {
@@ -261,7 +261,7 @@ Foam::radiation::greyMeanAbsorptionEmission::aCont(const label bandI) const
                 );
         }
     }
-    ta().correctBoundaryConditions();
+    ta.ref().correctBoundaryConditions();
     return ta;
 }
 
@@ -300,11 +300,11 @@ Foam::radiation::greyMeanAbsorptionEmission::ECont(const label bandI) const
 
         if (dQ.dimensions() == dimEnergy/dimTime)
         {
-            E().internalField() = EhrrCoeff_*dQ/mesh_.V();
+            E.ref().internalField() = EhrrCoeff_*dQ/mesh_.V();
         }
         else if (dQ.dimensions() == dimEnergy/dimTime/dimVolume)
         {
-            E().internalField() = EhrrCoeff_*dQ;
+            E.ref().internalField() = EhrrCoeff_*dQ;
         }
         else
         {
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
index d3da12eb593959773f843a706fff04ee6b874e76..0e2238eeecf3a8f5328bfb1c153dd73bade7c861 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
@@ -55,10 +55,10 @@ greyMeanSolidAbsorptionEmission::X(const word specie) const
     const volScalarField& p = thermo_.p();
 
     tmp<scalarField> tXj(new scalarField(T.internalField().size(), 0.0));
-    scalarField& Xj = tXj();
+    scalarField& Xj = tXj.ref();
 
     tmp<scalarField> tRhoInv(new scalarField(T.internalField().size(), 0.0));
-    scalarField& rhoInv = tRhoInv();
+    scalarField& rhoInv = tRhoInv.ref();
 
     forAll(mixture_.Y(), specieI)
     {
@@ -162,7 +162,7 @@ calc(const label propertyId) const
         )
     );
 
-    scalarField& a = ta().internalField();
+    scalarField& a = ta.ref().internalField();
 
     forAllConstIter(HashTable<label>, speciesNames_, iter)
     {
@@ -172,7 +172,7 @@ calc(const label propertyId) const
         }
     }
 
-    ta().correctBoundaryConditions();
+    ta.ref().correctBoundaryConditions();
     return ta;
 }
 
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
index 7147a96065aae1213e93327306a88c6dd42a97dc..673043c39bebba79e6343a789259a08af98b6929 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.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
@@ -177,7 +177,7 @@ Foam::radiation::wideBandAbsorptionEmission::aCont(const label bandI) const
         )
     );
 
-    scalarField& a = ta().internalField();
+    scalarField& a = ta.ref().internalField();
 
     forAll(a, i)
     {
@@ -255,7 +255,7 @@ Foam::radiation::wideBandAbsorptionEmission::ECont(const label bandI) const
 
         if (dQ.dimensions() == dimEnergy/dimTime)
         {
-            E().internalField() =
+            E.ref().internalField() =
                 iEhrrCoeffs_[bandI]
                *dQ.internalField()
                *(iBands_[bandI][1] - iBands_[bandI][0])
@@ -264,7 +264,7 @@ Foam::radiation::wideBandAbsorptionEmission::ECont(const label bandI) const
         }
         else if (dQ.dimensions() == dimEnergy/dimTime/dimVolume)
         {
-            E().internalField() =
+            E.ref().internalField() =
                 iEhrrCoeffs_[bandI]
                *dQ.internalField()
                *(iBands_[bandI][1] - iBands_[bandI][0])
diff --git a/src/thermophysicalModels/reactionThermo/Make/options b/src/thermophysicalModels/reactionThermo/Make/options
index f59f44fc8d3a1d08df9534646ba44212a3034ba2..bded2d7118f6d62b9e7b73c14ac287c60c1b4b9b 100644
--- a/src/thermophysicalModels/reactionThermo/Make/options
+++ b/src/thermophysicalModels/reactionThermo/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
diff --git a/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.C b/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.C
index d9062e03f496fcc56644b54e20990b137b9d7898..a3643c0d30b0dab0da6dbeacde183afadbb5d939 100644
--- a/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.C
+++ b/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,7 +68,7 @@ Foam::tmp<Foam::volScalarField> Foam::basicSpecieMixture::W() const
         )
     );
 
-    volScalarField& rW = trW();
+    volScalarField& rW = trW.ref();
 
     forAll(Y, i)
     {
diff --git a/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C b/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C
index 7d51293208e61512d33bad2dc9adfdd81e8dce94..5f963582dca0695fd8f13e1058cf774a240817f3 100644
--- a/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C
+++ b/src/thermophysicalModels/reactionThermo/psiuReactionThermo/heheuPsiThermo.C
@@ -229,7 +229,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::heu
 ) const
 {
     tmp<scalarField> theu(new scalarField(Tu.size()));
-    scalarField& heu = theu();
+    scalarField& heu = theu.ref();
 
     forAll(Tu, celli)
     {
@@ -250,7 +250,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::heu
 ) const
 {
     tmp<scalarField> theu(new scalarField(Tu.size()));
-    scalarField& heu = theu();
+    scalarField& heu = theu.ref();
 
     forAll(Tu, facei)
     {
@@ -283,7 +283,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::Tb() const
         )
     );
 
-    volScalarField& Tb_ = tTb();
+    volScalarField& Tb_ = tTb.ref();
     scalarField& TbCells = Tb_.internalField();
     const scalarField& pCells = this->p_.internalField();
     const scalarField& TCells = this->T_.internalField();
@@ -341,7 +341,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::psiu() const
         )
     );
 
-    volScalarField& psiu = tpsiu();
+    volScalarField& psiu = tpsiu.ref();
     scalarField& psiuCells = psiu.internalField();
     const scalarField& TuCells = this->Tu_.internalField();
     const scalarField& pCells = this->p_.internalField();
@@ -393,7 +393,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::psib() const
         )
     );
 
-    volScalarField& psib = tpsib();
+    volScalarField& psib = tpsib.ref();
     scalarField& psibCells = psib.internalField();
     const volScalarField Tb_(Tb());
     const scalarField& TbCells = Tb_.internalField();
@@ -446,7 +446,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::muu() const
         )
     );
 
-    volScalarField& muu_ = tmuu();
+    volScalarField& muu_ = tmuu.ref();
     scalarField& muuCells = muu_.internalField();
     const scalarField& pCells = this->p_.internalField();
     const scalarField& TuCells = this->Tu_.internalField();
@@ -502,7 +502,7 @@ Foam::heheuPsiThermo<BasicPsiThermo, MixtureType>::mub() const
         )
     );
 
-    volScalarField& mub_ = tmub();
+    volScalarField& mub_ = tmub.ref();
     scalarField& mubCells = mub_.internalField();
     const volScalarField Tb_(Tb());
     const scalarField& pCells = this->p_.internalField();
diff --git a/src/thermophysicalModels/solidChemistryModel/Make/options b/src/thermophysicalModels/solidChemistryModel/Make/options
index e871169b17f52e46e9fdf1ab5878f7c4f4532acc..e922aa724f5f1fa2eb49dd3b39eaebacfa6f38e3 100644
--- a/src/thermophysicalModels/solidChemistryModel/Make/options
+++ b/src/thermophysicalModels/solidChemistryModel/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
index efa78b4e149c85eaceb9f5d9f36da90471e247c4..dee571fe189b28777880b1dc5a3bdf8e8057565c 100644
--- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
@@ -637,7 +637,7 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::gasHs
         )
     );
 
-    volScalarField::InternalField& gasHs = tHs().internalField();
+    volScalarField::InternalField& gasHs = tHs.ref().internalField();
 
     const GasThermo& mixture = gasThermo_[index];
 
diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModelI.H b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModelI.H
index d045a3f894c53d2511fd793e5861d72ba51e77de..b14770004d1877946205549c95c4e29ddcf9ba34 100644
--- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModelI.H
+++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModelI.H
@@ -98,12 +98,15 @@ RRg() const
 
     if (this->chemistry_)
     {
-        DimensionedField<scalar, volMesh>& RRg = tRRg();
+        DimensionedField<scalar, volMesh>& RRg = tRRg.ref();
         for (label i=0; i < nGases_; i++)
         {
             RRg += RRg_[i];
         }
     }
+
     return tRRg;
 }
+
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
index 07cc0c3ae3185008fac145d9f648886b3888357e..8bf5ea932099bfa68411db0522d00f5058258700 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
@@ -136,7 +136,7 @@ Foam::solidChemistryModel<CompType, SolidThermo>::Sh() const
 
     if (this->chemistry_)
     {
-        scalarField& Sh = tSh();
+        scalarField& Sh = tSh.ref();
 
         forAll(Ys_, i)
         {
@@ -176,7 +176,7 @@ Foam::solidChemistryModel<CompType, SolidThermo>::dQ() const
 
     if (this->chemistry_)
     {
-        volScalarField& dQ = tdQ();
+        volScalarField& dQ = tdQ.ref();
         dQ.dimensionedInternalField() = this->mesh_.V()*Sh()();
     }
 
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
index 6223a850869c409e5c754981b872e7486c731ece..43a343b5f87196302ea376ef770a21a7510a52cd 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
@@ -85,7 +85,7 @@ Foam::solidChemistryModel<CompType, SolidThermo>::RRs() const
 
     if (this->chemistry_)
     {
-        DimensionedField<scalar, volMesh>& RRs = tRRs();
+        DimensionedField<scalar, volMesh>& RRs = tRRs.ref();
         for (label i=0; i < nSolids_; i++)
         {
             RRs += RRs_[i];
diff --git a/src/thermophysicalModels/solidSpecie/Make/options b/src/thermophysicalModels/solidSpecie/Make/options
index bdc8b074c03a509ce2f001f40aff6c0c189b0222..85a93d1d170df8cdd0fe01f707948075b80fe868 100644
--- a/src/thermophysicalModels/solidSpecie/Make/options
+++ b/src/thermophysicalModels/solidSpecie/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
 
 LIB_LIBS = \
diff --git a/src/thermophysicalModels/solidThermo/Make/options b/src/thermophysicalModels/solidThermo/Make/options
index d3e1cb87ce5eac43f9599cf03d7d2c861d24bfab..c84b6c55464449cf1fd5e8853db822c7a96f41c2 100644
--- a/src/thermophysicalModels/solidThermo/Make/options
+++ b/src/thermophysicalModels/solidThermo/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
diff --git a/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C b/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
index 75de2bc00314a57260d1d9ee843f493a63553150..7928d9bdff7c38bf4450d1d793aa2d14e6e15063 100644
--- a/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
+++ b/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
@@ -203,7 +203,7 @@ Foam::heSolidThermo<BasicSolidThermo, MixtureType>::Kappa() const
         )
     );
 
-    volVectorField& Kappa = tKappa();
+    volVectorField& Kappa = tKappa.ref();
     vectorField& KappaCells = Kappa.internalField();
     const scalarField& TCells = this->T_.internalField();
     const scalarField& pCells = this->p_.internalField();
@@ -253,7 +253,7 @@ Foam::heSolidThermo<BasicSolidThermo, MixtureType>::Kappa
     const scalarField& Tp = this->T_.boundaryField()[patchi];
     tmp<vectorField> tKappa(new vectorField(pp.size()));
 
-    vectorField& Kappap = tKappa();
+    vectorField& Kappap = tKappa.ref();
 
     forAll(Tp, facei)
     {
diff --git a/src/topoChangerFvMesh/Make/options b/src/topoChangerFvMesh/Make/options
index 44753e64c96546ae740c2d54605959f8ba6890fd..82165ae6e5ec39ca860ac3ef18b61976d679f3ab 100644
--- a/src/topoChangerFvMesh/Make/options
+++ b/src/topoChangerFvMesh/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
diff --git a/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C b/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C
index 02d6a6ae87896e9a8613a3fc1a77cc8fc9963b65..747a1c0edf21695b04533e8d9eadbb36710d07f0 100644
--- a/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C
+++ b/src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C
@@ -59,7 +59,7 @@ Foam::tmp<Foam::scalarField> Foam::movingConeTopoFvMesh::vertexMarkup
         << curLeft << " curRight: " << curRight << endl;
 
     tmp<scalarField> tvertexMarkup(new scalarField(p.size()));
-    scalarField& vertexMarkup = tvertexMarkup();
+    scalarField& vertexMarkup = tvertexMarkup.ref();
 
     forAll(p, pI)
     {
diff --git a/src/transportModels/compressible/Make/options b/src/transportModels/compressible/Make/options
index 71b7873964d544eddf96d22aa40f4c3372c23c9c..287318da1ff31d9067bc255cafbdd1fa8dce4787 100644
--- a/src/transportModels/compressible/Make/options
+++ b/src/transportModels/compressible/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
diff --git a/src/transportModels/immiscibleIncompressibleTwoPhaseMixture/Make/options b/src/transportModels/immiscibleIncompressibleTwoPhaseMixture/Make/options
index 4ae2fbf96644b3ac1443ae0bea87926356d158fd..138dcbf7f55763cb937eb1952e92913782558664 100644
--- a/src/transportModels/immiscibleIncompressibleTwoPhaseMixture/Make/options
+++ b/src/transportModels/immiscibleIncompressibleTwoPhaseMixture/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I.. \
     -I../incompressible/lnInclude \
     -I../interfaceProperties/lnInclude \
diff --git a/src/transportModels/incompressible/Make/options b/src/transportModels/incompressible/Make/options
index 46aa26d4ad8980819f6d6b939101f3f3fef92f64..c55a98a0b64a5facfe106252a638c7a3772b4271 100644
--- a/src/transportModels/incompressible/Make/options
+++ b/src/transportModels/incompressible/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I.. \
     -I../twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/src/transportModels/interfaceProperties/Make/options b/src/transportModels/interfaceProperties/Make/options
index cc63db7dcc45c854ad4a57862b1b9ad635ec5ce7..efbf1b41ca3431db1b4527a8b8f188a73ad239f7 100644
--- a/src/transportModels/interfaceProperties/Make/options
+++ b/src/transportModels/interfaceProperties/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 71b7873964d544eddf96d22aa40f4c3372c23c9c..287318da1ff31d9067bc255cafbdd1fa8dce4787 100644
--- a/src/transportModels/twoPhaseMixture/Make/options
+++ b/src/transportModels/twoPhaseMixture/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
diff --git a/src/transportModels/twoPhaseProperties/Make/options b/src/transportModels/twoPhaseProperties/Make/options
index 32ce36b9594d8e970ba45526b49fce90d0d9fef0..547d45694c84a3d980b5c1b88dc36aaea3d2e224 100644
--- a/src/transportModels/twoPhaseProperties/Make/options
+++ b/src/transportModels/twoPhaseProperties/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -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 9ee5884e5908ccdb6b1ba3acf323b8c524ec996a..948c54f7560fdc0d13aef5830827846f1d0dd459 100644
--- a/src/triSurface/Make/options
+++ b/src/triSurface/Make/options
@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -DCONST_TMP \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude
 
diff --git a/src/triSurface/faceTriangulation/faceTriangulation.C b/src/triSurface/faceTriangulation/faceTriangulation.C
index ddf905f06ba9e29434dd1c520679385be103a9f5..6b44c53065401e24ac4a0db9e056d56bc9175ced 100644
--- a/src/triSurface/faceTriangulation/faceTriangulation.C
+++ b/src/triSurface/faceTriangulation/faceTriangulation.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
@@ -55,7 +55,7 @@ Foam::tmp<Foam::vectorField> Foam::faceTriangulation::calcEdges
 )
 {
     tmp<vectorField> tedges(new vectorField(f.size()));
-    vectorField& edges = tedges();
+    vectorField& edges = tedges.ref();
 
     forAll(f, i)
     {