diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C
index d2117d1334..7e5fa4f286 100644
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C
@@ -69,7 +69,8 @@ Foam::XiEqModels::Gulder::~Gulder()
 Foam::tmp<Foam::volScalarField> Foam::XiEqModels::Gulder::XiEq() const
 {
     volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
-    const volScalarField& epsilon = turbulence_.epsilon();
+    const tmp<volScalarField> tepsilon(turbulence_.epsilon());
+    const volScalarField& epsilon = tepsilon();
 
     if (subGridSchelkin_)
     {
diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
index d482f37983..b80b72a4c7 100644
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
@@ -75,8 +75,10 @@ Foam::XiEqModels::SCOPEXiEq::~SCOPEXiEq()
 
 Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
 {
-    const volScalarField& k = turbulence_.k();
-    const volScalarField& epsilon = turbulence_.epsilon();
+    const tmp<volScalarField> tk(turbulence_.k());
+    const volScalarField& k = tk();
+    const tmp<volScalarField> tepsilon(turbulence_.epsilon());
+    const volScalarField& epsilon = tepsilon();
 
     volScalarField up(sqrt((2.0/3.0)*k));
     if (subGridSchelkin_)
diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.C b/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.C
index c586ce2914..1b44e52069 100644
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.C
@@ -66,7 +66,8 @@ Foam::XiGModels::KTS::~KTS()
 Foam::tmp<Foam::volScalarField> Foam::XiGModels::KTS::G() const
 {
     volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
-    const volScalarField& epsilon = turbulence_.epsilon();
+    const tmp<volScalarField> tepsilon(turbulence_.epsilon());
+    const volScalarField& epsilon = tepsilon();
 
     volScalarField tauEta(sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon))));
 
diff --git a/applications/solvers/compressible/rhoCentralFoam/createFieldRefs.H b/applications/solvers/compressible/rhoCentralFoam/createFieldRefs.H
index 5522ccd6aa..f5a13cc982 100644
--- a/applications/solvers/compressible/rhoCentralFoam/createFieldRefs.H
+++ b/applications/solvers/compressible/rhoCentralFoam/createFieldRefs.H
@@ -1,10 +1,9 @@
 volScalarField& p = thermo.p();
 const volScalarField& T = thermo.T();
 const volScalarField& psi = thermo.psi();
-const volScalarField& mu = thermo.mu();
 
 bool inviscid(true);
-if (max(mu.primitiveField()) > 0.0)
+if (max(thermo.mu().cref().primitiveField()) > 0.0)
 {
     inviscid = false;
 }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C
index f3af8611d3..b0c52f0f62 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C
@@ -160,7 +160,8 @@ void VoFPatchTransfer::correct
     const volScalarField& heVoF = thermo.thermo1().he();
     const volScalarField& TVoF = thermo.thermo1().T();
     const volScalarField CpVoF(thermo.thermo1().Cp());
-    const volScalarField& rhoVoF = thermo.thermo1().rho()();
+    const tmp<volScalarField> trhoVoF(thermo.thermo1().rho());
+    const volScalarField& rhoVoF = trhoVoF();
     const volScalarField& alphaVoF = thermo.alpha1();
 
     forAll(patchIDs_, pidi)
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
index 107fe3c0cc..a447b1eca8 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
@@ -28,7 +28,8 @@
     forAllConstIters(mixture.phases(), phase)
     {
         const rhoThermo& thermo = phase().thermo();
-        const volScalarField& rho = thermo.rho()();
+        const tmp<volScalarField> trho(thermo.rho());
+        const volScalarField& rho = trho();
 
         p_rghEqnComps.set
         (
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C
index c8b8dc2280..438125651d 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C
@@ -482,8 +482,8 @@ void Foam::searchableSurfaceControl::cellSizeFunctionVertices
     DynamicList<scalar>& sizes
 ) const
 {
-    const tmp<pointField> tmpPoints = searchableSurface_.points();
-    const pointField& points = tmpPoints();
+    const tmp<pointField> tpoints(searchableSurface_.points());
+    const pointField& points = tpoints();
 
     const scalar nearFeatDistSqrCoeff = 1e-8;
 
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C
index 412ef9bb55..5d216568d3 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C
@@ -110,7 +110,8 @@ Foam::scalar Foam::nonUniformField::interpolate
 {
     const face& faceHitByPt = surfaceTriMesh_.triSurface::operator[](index);
 
-    const pointField& pts = surfaceTriMesh_.points();
+    const tmp<pointField> tpoints(surfaceTriMesh_.points());
+    const pointField& pts = tpoints();
 //    const Map<label>& pMap = surfaceTriMesh_.meshPointMap();
 
     triPointRef tri
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index c1cebc026a..a90fa8ac51 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -982,7 +982,7 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::findOffsetPatchFaces
         offsetBoundaryCells.write();
     }
 
-    return std::move(offsetBoundaryCells);
+    return offsetBoundaryCells;
 }
 
 
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 03f37ea67d..ccc5b21fe9 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -1371,7 +1371,7 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::findRemainingProtrusionSet
         protrudingCells.write();
     }
 
-    return std::move(protrudingCells);
+    return protrudingCells;
 }
 
 
diff --git a/applications/utilities/surface/surfaceInflate/surfaceInflate.C b/applications/utilities/surface/surfaceInflate/surfaceInflate.C
index 364b00e0e7..69adf9f171 100644
--- a/applications/utilities/surface/surfaceInflate/surfaceInflate.C
+++ b/applications/utilities/surface/surfaceInflate/surfaceInflate.C
@@ -221,9 +221,10 @@ void detectSelfIntersections
     const edgeList& edges = s.edges();
     const indexedOctree<treeDataTriSurface>& tree = s.tree();
     const labelList& meshPoints = s.meshPoints();
-    const pointField& points = s.points();
+    const tmp<pointField> tpoints(s.points());
+    const pointField& points = tpoints();
 
-    isEdgeIntersecting.setSize(edges.size());
+    isEdgeIntersecting.resize_nocopy(edges.size());
     isEdgeIntersecting = false;
 
     forAll(edges, edgeI)
@@ -311,7 +312,8 @@ label detectIntersectionPoints
         detectSelfIntersections(s, isEdgeIntersecting);
 
         const edgeList& edges = s.edges();
-        const pointField& points = s.points();
+        const tmp<pointField> tpoints(s.points());
+        const pointField& points = tpoints();
 
         forAll(edges, edgeI)
         {
@@ -836,9 +838,10 @@ int main(int argc, char *argv[])
         // Do some smoothing (Lloyds algorithm)
         lloydsSmoothing(nSmooth, s, isFeaturePoint, edgeStat, isAffectedPoint);
 
-
         // Update pointDisplacement
-        const pointField& pts = s.points();
+        const tmp<pointField> tpoints(s.points());
+        const pointField& pts = tpoints();
+
         forAll(meshPoints, i)
         {
             label meshPointI = meshPoints[i];
diff --git a/src/OpenFOAM/fields/cloud/cloud.H b/src/OpenFOAM/fields/cloud/cloud.H
index ed1109cdea..e0256fd6dc 100644
--- a/src/OpenFOAM/fields/cloud/cloud.H
+++ b/src/OpenFOAM/fields/cloud/cloud.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -172,6 +172,21 @@ public:
             {
                 return obr.lookupObject<IOField<Type>>(fieldName);
             }
+
+            //- Lookup an IOField within object registry
+            //  Fatal if not found or wrong type
+            //
+            //  Note: const char signature to avoid spurious
+            //  -Wdangling-reference with gcc-13
+            template<class Type>
+            static const IOField<Type>& lookupIOField
+            (
+                const char* fieldName,
+                const objectRegistry& obr
+            )
+            {
+                return obr.lookupObject<IOField<Type>>(word(fieldName));
+            }
 };
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
index 9341a5ad12..5813f534cf 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
@@ -289,10 +289,12 @@ void Foam::processorPolyPatch::calcGeometry(PstreamBuffers& pBufs)
 
                 label vertI = 0;
 
-                forAll(faceCentres(), facej)
+                const vectorField::subField ownFaceCentres = faceCentres();
+
+                forAll(ownFaceCentres, facej)
                 {
                     const point& c0 = neighbFaceCentres_[facej];
-                    const point& c1 = faceCentres()[facej];
+                    const point& c1 = ownFaceCentres[facej];
 
                     writeOBJ(ccStr, c0, c1, vertI);
                 }
diff --git a/src/atmosphericModels/fvOptions/atmAmbientTurbSource/atmAmbientTurbSourceTemplates.C b/src/atmosphericModels/fvOptions/atmAmbientTurbSource/atmAmbientTurbSourceTemplates.C
index e1dcf80748..d4ef78c3cd 100644
--- a/src/atmosphericModels/fvOptions/atmAmbientTurbSource/atmAmbientTurbSourceTemplates.C
+++ b/src/atmosphericModels/fvOptions/atmAmbientTurbSource/atmAmbientTurbSourceTemplates.C
@@ -45,7 +45,8 @@ void Foam::fv::atmAmbientTurbSource::atmAmbientTurbSourceEpsilon
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& epsilon = turbPtr->epsilon();
+    const tmp<volScalarField> tepsilon(turbPtr->epsilon());
+    const volScalarField& epsilon = tepsilon();
 
     // (Heuristically derived from RS:Eq. 4, rhs-term:5)
     eqn +=
@@ -67,7 +68,9 @@ void Foam::fv::atmAmbientTurbSource::atmAmbientTurbSourceOmega
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& omega = turbPtr->omega();
+    const tmp<volScalarField> tomega(turbPtr->omega());
+    const volScalarField& omega = tomega();
+
     const volScalarField::Internal& beta =
         mesh_.lookupObjectRef<volScalarField::Internal>
         (
@@ -93,7 +96,8 @@ void Foam::fv::atmAmbientTurbSource::atmAmbientTurbSourceK
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& k = turbPtr->k();
+    const tmp<volScalarField> tk(turbPtr->k());
+    const volScalarField& k = tk();
 
     if (isEpsilon_)
     {
diff --git a/src/atmosphericModels/fvOptions/atmBuoyancyTurbSource/atmBuoyancyTurbSourceTemplates.C b/src/atmosphericModels/fvOptions/atmBuoyancyTurbSource/atmBuoyancyTurbSourceTemplates.C
index 468266ee3e..c6a1e7f9b6 100644
--- a/src/atmosphericModels/fvOptions/atmBuoyancyTurbSource/atmBuoyancyTurbSourceTemplates.C
+++ b/src/atmosphericModels/fvOptions/atmBuoyancyTurbSource/atmBuoyancyTurbSourceTemplates.C
@@ -47,8 +47,10 @@ void Foam::fv::atmBuoyancyTurbSource::atmBuoyancyTurbSourceEpsilon
         );
 
     // Fetch required fields from the epsilon-based model
-    const volScalarField& k = turbPtr->k();
-    const volScalarField& epsilon = turbPtr->epsilon();
+    const tmp<volScalarField> tk(turbPtr->k());
+    const volScalarField& k = tk();
+    const tmp<volScalarField> tepsilon(turbPtr->epsilon());
+    const volScalarField& epsilon = tepsilon();
     const volScalarField::Internal& GbyNu =
         mesh_.lookupObjectRef<volScalarField::Internal>
         (
@@ -77,8 +79,10 @@ void Foam::fv::atmBuoyancyTurbSource::atmBuoyancyTurbSourceOmega
         );
 
     // Fetch required fields from the omega-based model
-    const volScalarField& k = turbPtr->k();
-    const volScalarField& omega = turbPtr->omega();
+    const tmp<volScalarField> tk(turbPtr->k());
+    const volScalarField& k = tk();
+    const tmp<volScalarField> tomega(turbPtr->omega());
+    const volScalarField& omega = tomega();
     const volScalarField::Internal& GbyNu =
         mesh_.lookupObjectRef<volScalarField::Internal>
         (
@@ -121,7 +125,8 @@ void Foam::fv::atmBuoyancyTurbSource::atmBuoyancyTurbSourceK
             turbulenceModel::propertiesName
         );
 
-    const volScalarField& k = turbPtr->k();
+    const tmp<volScalarField> tk(turbPtr->k());
+    const volScalarField& k = tk();
 
     eqn += fvm::Sp(alpha()*rho()*B_/k(), k);
 }
diff --git a/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSourceTemplates.C b/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSourceTemplates.C
index 0627313a52..f7125eab8b 100644
--- a/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSourceTemplates.C
+++ b/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSourceTemplates.C
@@ -45,7 +45,8 @@ void Foam::fv::atmPlantCanopyTurbSource::atmPlantCanopyTurbSourceEpsilon
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& epsilon = turbPtr->epsilon();
+    const tmp<volScalarField> tepsilon(turbPtr->epsilon());
+    const volScalarField& epsilon = tepsilon();
     const volVectorField::Internal& U = turbPtr->U()();
 
     eqn -= fvm::Sp(alpha()*rho()*(C1_ - C2_)*calcPlantCanopyTerm(U), epsilon);
@@ -66,7 +67,8 @@ void Foam::fv::atmPlantCanopyTurbSource::atmPlantCanopyTurbSourceOmega
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& omega = turbPtr->omega();
+    const tmp<volScalarField> tomega(turbPtr->omega());
+    const volScalarField& omega = tomega();
     const volVectorField::Internal& U = turbPtr->U()();
     const volScalarField::Internal& gamma =
         mesh_.lookupObjectRef<volScalarField::Internal>
diff --git a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.H b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.H
index cdde0dc187..c5b2a4417d 100644
--- a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.H
+++ b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.H
@@ -38,8 +38,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef singleCellFvMesh_H
-#define singleCellFvMesh_H
+#ifndef Foam_singleCellFvMesh_H
+#define Foam_singleCellFvMesh_H
 
 #include "fvPatchFieldMapper.H"
 #include "fvMesh.H"
@@ -58,7 +58,7 @@ class singleCellFvMesh
 :
     public fvMesh
 {
-    // Private data
+    // Private Data
 
         const labelListIOList patchFaceAgglomeration_;
 
@@ -179,36 +179,39 @@ public:
 
     // Member Functions
 
-        bool agglomerate() const
+        bool agglomerate() const noexcept
         {
-            return patchFaceAgglomeration_.size() > 0;
+            return !patchFaceAgglomeration_.empty();
         }
 
         //- From patchFace on this back to original mesh or agglomeration
-        const labelListList& patchFaceMap() const
+        const labelListList& patchFaceMap() const noexcept
         {
             return patchFaceMap_;
         }
 
         //- From point on this back to original mesh
-        const labelList& pointMap() const
+        const labelList& pointMap() const noexcept
         {
             return pointMap_;
         }
 
         //- From face on original mesh to face on this
-        const labelList& reverseFaceMap() const
+        const labelList& reverseFaceMap() const noexcept
         {
             return reverseFaceMap_;
         }
 
         //- From point on original mesh to point on this (or -1 for removed
         //- points)
-        const labelList& reversePointMap() const
+        const labelList& reversePointMap() const noexcept
         {
             return reversePointMap_;
         }
 
+        //- Interpolate for overset (unused)
+        using fvMesh::interpolate;
+
         //- Map volField. Internal field set to average, patch fields straight
         //- copies.
         template<class Type>
@@ -217,7 +220,6 @@ public:
         (
             const GeometricField<Type, fvPatchField, volMesh>&
         ) const;
-
 };
 
 
diff --git a/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C b/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C
index 24dfae4379..ce710d6c25 100644
--- a/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C
+++ b/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C
@@ -132,9 +132,9 @@ void Foam::patchTransformedInterpolation::interpolate
     labelList pointDisplacementNSum(nPoints, Zero);
     vectorField pointDisplacementSum(nPoints, Zero);
 
-    forAll(patches_, patchI)
+    for (const label patchi : patches_)
     {
-        const polyPatch& patch(mesh().boundaryMesh()[patches_[patchI]]);
+        const polyPatch& patch = mesh().boundaryMesh()[patchi];
 
         forAll(patch, pFaceI)
         {
@@ -145,7 +145,7 @@ void Foam::patchTransformedInterpolation::interpolate
             const labelList cPoints(c.labels(mesh().faces()));
 
             // Consider movement around the face centre
-            const point& xOrigin(patch.faceCentres()[pFaceI]);
+            const point xOrigin(patch.faceCentres()[pFaceI]);
 
             // Mean translation
             const vector uMean(f.average(points, pointDisplacement));
diff --git a/src/fvOptions/sources/derived/buoyancyTurbSource/buoyancyTurbSourceTemplates.C b/src/fvOptions/sources/derived/buoyancyTurbSource/buoyancyTurbSourceTemplates.C
index 10cdc62b8c..db5ef7b800 100644
--- a/src/fvOptions/sources/derived/buoyancyTurbSource/buoyancyTurbSourceTemplates.C
+++ b/src/fvOptions/sources/derived/buoyancyTurbSource/buoyancyTurbSourceTemplates.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd
+    Copyright (C) 2020,2023 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,7 +46,8 @@ void Foam::fv::buoyancyTurbSource::buoyancyTurbSourceK
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& nut = turbPtr->nut();
+    const tmp<volScalarField> tnut(turbPtr->nut());
+    const volScalarField& nut = tnut();
 
     const dictionary& turbDict = turbPtr->coeffDict();
     const scalar Prt
diff --git a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index 12fd8ab036..aab3ce1bd8 100644
--- a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -231,7 +231,7 @@ void Foam::FreeStream<CloudType>::inflow()
 
             label celli = mesh.faceOwner()[globalFaceIndex];
 
-            const vector& fC = patch.faceCentres()[pFI];
+            const vector fC = patch.faceCentres()[pFI];
 
             scalar fA = mag(patch.faceAreas()[pFI]);
 
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
index 8dad26141a..8066ce6951 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
@@ -3805,26 +3805,27 @@ const Foam::dictionary& Foam::meshRefinement::subDict
     enum keyType::option matchOpt
 )
 {
-    const auto finder(dict.csearch(keyword, matchOpt));
+    const dictionary* dictptr = dict.findDict(keyword, matchOpt);
 
-    if (!finder.good())
+    if (!dictptr)
     {
-        auto& err = FatalIOErrorInFunction(dict);
-
-        err << "Entry '" << keyword << "' not found in dictionary "
+        FatalIOErrorInFunction(dict)
+            << "Entry '" << keyword
+            << "' not found (or not a dictionary) in dictionary "
             << dict.relativeName() << nl;
 
         if (noExit)
         {
+            // Dummy return
             return dictionary::null;
         }
         else
         {
-            err << exit(FatalIOError);
+            FatalIOError << exit(FatalIOError);
         }
     }
 
-    return finder.dict();
+    return *dictptr;
 }
 
 
@@ -3840,19 +3841,18 @@ Foam::ITstream& Foam::meshRefinement::lookup
 
     if (!eptr)
     {
-        auto& err = FatalIOErrorInFunction(dict);
-
-        err << "Entry '" << keyword << "' not found in dictionary "
+        FatalIOErrorInFunction(dict)
+            << "Entry '" << keyword << "' not found in dictionary "
             << dict.relativeName() << nl;
 
         if (noExit)
         {
-            // Fake entry
-            return dict.first()->stream();
+            // Dummy return
+            return ITstream::empty_stream();
         }
         else
         {
-            err << exit(FatalIOError);
+            FatalIOError << exit(FatalIOError);
         }
     }
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/adjointSensitivity/adjointEikonalSolver/adjointEikonalSolver.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/adjointSensitivity/adjointEikonalSolver/adjointEikonalSolver.C
index 64735b5080..37e6b9d942 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/adjointSensitivity/adjointEikonalSolver/adjointEikonalSolver.C
+++ b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/adjointSensitivity/adjointEikonalSolver/adjointEikonalSolver.C
@@ -83,7 +83,8 @@ void adjointEikonalSolver::read()
 tmp<surfaceScalarField> adjointEikonalSolver::computeYPhi()
 {
     // Primal distance field
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
 
     volVectorField ny
     (
@@ -198,7 +199,8 @@ void adjointEikonalSolver::solve()
     read();
 
     // Primal distance field
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
 
     // Convecting flux
     tmp<surfaceScalarField> tyPhi = computeYPhi();
@@ -275,7 +277,9 @@ boundaryVectorField& adjointEikonalSolver::distanceSensitivities()
 
     boundaryVectorField& distanceSens = distanceSensPtr_();
 
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
+
     for (const label patchi : sensitivityPatchIDs_)
     {
         vectorField nf(mesh_.boundary()[patchi].nf());
@@ -294,7 +298,8 @@ tmp<volTensorField> adjointEikonalSolver::getFISensitivityTerm() const
 {
     Info<< "Calculating distance sensitivities " << endl;
 
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
     const volVectorField gradD(fvc::grad(d));
 
     auto gradDDa
@@ -354,7 +359,8 @@ tmp<scalarField> adjointEikonalSolver::topologySensitivities
     const word& designVarsName
 ) const
 {
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
 
     auto tres(tmp<scalarField>::New(d.primitiveField().size(), Zero));
     scalarField dSens(d.primitiveField()*da_.primitiveField());
@@ -377,7 +383,9 @@ const volScalarField& adjointEikonalSolver::da()
 
 tmp<volVectorField> adjointEikonalSolver::gradEikonal()
 {
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
+
     volVectorField gradD(fvc::grad(d));
     return tmp<volVectorField>::New("gradEikonal", 2*gradD & fvc::grad(gradD));
 }
diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/SpalartAllmaras/SpalartAllmaras.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/SpalartAllmaras/SpalartAllmaras.C
index 9e17b38e59..c583938e48 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/SpalartAllmaras/SpalartAllmaras.C
@@ -95,7 +95,8 @@ tmp<volScalarField> SpalartAllmaras::nutJacobianVar1
     );
     auto& nutJacobian = tnutJacobian.ref();
 
-    const volScalarField& nu = laminarTransport.nu();
+    const tmp<volScalarField> tnu(laminarTransport.nu());
+    const volScalarField& nu = tnu();
     const volScalarField& nuTilda = TMVar1();
 
     volScalarField chi(nuTilda/nu);
diff --git a/src/phaseSystemModels/multiphaseEuler/multiphaseSystem/phaseModel/phaseModel.C b/src/phaseSystemModels/multiphaseEuler/multiphaseSystem/phaseModel/phaseModel.C
index 7f1270c587..0da2401b99 100644
--- a/src/phaseSystemModels/multiphaseEuler/multiphaseSystem/phaseModel/phaseModel.C
+++ b/src/phaseSystemModels/multiphaseEuler/multiphaseSystem/phaseModel/phaseModel.C
@@ -241,7 +241,8 @@ void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
 {
     surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
     const volScalarField::Boundary& alphaBf = boundaryField();
-    const surfaceScalarField::Boundary& phiBf = phi().boundaryField();
+    const tmp<surfaceScalarField> tphi(phi());
+    const auto& phiBf = tphi().boundaryField();
 
     forAll(alphaPhiBf, patchi)
     {
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index 7623bd2d70..fb7c537721 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -355,7 +355,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer()
                 fvVectorMatrix& eqn = *eqns[phase.name()];
 
                 const volVectorField& U = eqn.psi();
-                const surfaceScalarField& phi = phase.phi();
+                const tmp<surfaceScalarField> tphi(phase.phi());
+                const surfaceScalarField& phi = tphi();
 
                 eqn -=
                     Vm
@@ -408,7 +409,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransferf()
 
         if (!phase.stationary())
         {
-            const volVectorField& U = phase.U();
+            const tmp<volVectorField> tU(phase.U());
+            const volVectorField& U = tU();
 
             UgradUs.set
             (
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/interfacialModels/dragModels/segregated/segregated.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/interfacialModels/dragModels/segregated/segregated.C
index e3b8208670..467f1817ae 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/interfacialModels/dragModels/segregated/segregated.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/interfacialModels/dragModels/segregated/segregated.C
@@ -85,8 +85,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
     const volScalarField& alpha1(pair_.phase1());
     const volScalarField& alpha2(pair_.phase2());
 
-    const volScalarField& rho1(pair_.phase1().rho());
-    const volScalarField& rho2(pair_.phase2().rho());
+    const tmp<volScalarField> trho1(pair_.phase1().rho());
+    const tmp<volScalarField> trho2(pair_.phase2().rho());
+
+    const volScalarField& rho1 = trho1();
+    const volScalarField& rho2 = trho2();
 
     tmp<volScalarField> tnu1(pair_.phase1().nu());
     tmp<volScalarField> tnu2(pair_.phase2().nu());
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
index 54b2f2fbcd..56e2b855db 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
@@ -93,7 +93,8 @@ Foam::tmp<Foam::fvScalarMatrix>
 Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 {
     const volScalarField& alpha = *this;
-    const volScalarField& rho = this->rho();
+    const tmp<volScalarField> trho(this->rho());
+    const volScalarField& rho(trho());
 
     const tmp<volVectorField> tU(this->U());
     const volVectorField& U(tU());
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MovingPhaseModel/MovingPhaseModel.C
index 4e75adec97..c8d030117d 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MovingPhaseModel/MovingPhaseModel.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MovingPhaseModel/MovingPhaseModel.C
@@ -283,7 +283,8 @@ Foam::tmp<Foam::fvVectorMatrix>
 Foam::MovingPhaseModel<BasePhaseModel>::UEqn()
 {
     const volScalarField& alpha = *this;
-    const volScalarField& rho = this->thermo().rho();
+    const tmp<volScalarField> trho = this->thermo().rho();
+    const volScalarField& rho = trho();
 
     return
     (
@@ -303,14 +304,13 @@ Foam::MovingPhaseModel<BasePhaseModel>::UfEqn()
     // As the "normal" U-eqn but without the ddt terms
 
     const volScalarField& alpha = *this;
-    const volScalarField& rho = this->thermo().rho();
 
     return
     (
         fvm::div(alphaRhoPhi_, U_)
       - fvm::Sp(fvc::div(alphaRhoPhi_), U_)
       + fvm::SuSp(- this->continuityErrorSources(), U_)
-      + this->fluid().MRF().DDt(alpha*rho, U_)
+      + this->fluid().MRF().DDt(alpha*this->thermo().rho(), U_)
       + turbulence_->divDevRhoReff(U_)
     );
 }
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
index 9482f0244c..40c5464098 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
@@ -152,7 +152,8 @@ Foam::MultiComponentPhaseModel<BasePhaseModel>::YiEqn(volScalarField& Yi)
 {
     const volScalarField& alpha = *this;
     const surfaceScalarField alphaRhoPhi(this->alphaRhoPhi());
-    const volScalarField& rho = this->thermo().rho();
+    const tmp<volScalarField> trho(this->thermo().rho());
+    const volScalarField& rho = trho();
 
     return
     (
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/phaseModel/phaseModel.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/phaseModel/phaseModel.C
index 9cb69181e4..d697c19a38 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/phaseModel/phaseModel.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/phaseModel/phaseModel.C
@@ -207,7 +207,8 @@ void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
 {
     surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
     const volScalarField::Boundary& alphaBf = boundaryField();
-    const surfaceScalarField::Boundary& phiBf = phi()().boundaryField();
+    const tmp<surfaceScalarField> tphi(phi());
+    const auto& phiBf = tphi().boundaryField();
 
     forAll(alphaPhiBf, patchi)
     {
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C
index a38aaa2824..445162fcb6 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C
@@ -154,7 +154,8 @@ Foam::diameterModels::nucleationModels::wallBoiling::addToNucleationRate
 {
     const sizeGroup& fi = popBal_.sizeGroups()[i];
     const phaseModel& phase = fi.phase();
-    const volScalarField& rho = phase.rho();
+    const tmp<volScalarField> trho(phase.rho());
+    const volScalarField& rho = trho();
     const tmp<volScalarField> talphat(turbulence_.alphat());
     const volScalarField::Boundary& alphatBf = talphat().boundaryField();
 
diff --git a/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
index db4ccbfac4..e3e046cf44 100644
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
@@ -1273,7 +1273,8 @@ void Foam::diameterModels::populationBalanceModel::solve()
                 const phaseModel& phase = fi.phase();
                 const volScalarField& alpha = phase;
                 const dimensionedScalar& residualAlpha = phase.residualAlpha();
-                const volScalarField& rho = phase.thermo().rho();
+                const tmp<volScalarField> trho(phase.thermo().rho());
+                const volScalarField& rho = trho();
 
                 fvScalarMatrix sizeGroupEqn
                 (
diff --git a/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C b/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C
index 4bde1a14c0..cf5b7670bd 100644
--- a/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C
+++ b/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C
@@ -156,7 +156,8 @@ JohnsonJacksonSchaeffer::nu
     }
 
     const fvPatchList& patches = phase.mesh().boundary();
-    const volVectorField& U = phase.U();
+    const tmp<volVectorField> tU(phase.U());
+    const volVectorField& U = tU();
 
     volScalarField::Boundary& nufBf = nuf.boundaryFieldRef();
 
diff --git a/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
index b0aac1aa02..d0c5de3052 100644
--- a/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+++ b/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
@@ -154,7 +154,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu
     }
 
     const fvPatchList& patches = phase.mesh().boundary();
-    const volVectorField& U = phase.U();
+    const tmp<volVectorField> tU(phase.U());
+    const volVectorField& U = tU();
 
     volScalarField::Boundary& nufBf = nuf.boundaryFieldRef();
 
diff --git a/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
index 9ac98b1c21..84a80a1a1c 100644
--- a/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
+++ b/src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
@@ -266,7 +266,8 @@ Foam::RASModels::kineticTheoryModel::R() const
 Foam::tmp<Foam::volScalarField>
 Foam::RASModels::kineticTheoryModel::pPrime() const
 {
-    const volScalarField& rho = phase_.rho();
+    const tmp<volScalarField> trho(phase_.rho());
+    const volScalarField& rho = trho();
 
     tmp<volScalarField> tpPrime
     (
@@ -365,11 +366,13 @@ void Foam::RASModels::kineticTheoryModel::correct()
 {
     // Local references
     volScalarField alpha(max(alpha_, scalar(0)));
-    const volScalarField& rho = phase_.rho();
+    const tmp<volScalarField> trho(phase_.rho());
+    const volScalarField& rho = trho();
     const surfaceScalarField& alphaRhoPhi = alphaRhoPhi_;
     const volVectorField& U = U_;
-    const volVectorField& Uc_ =
+    const tmp<volVectorField> tUc =
         refCast<const twoPhaseSystem>(phase_.fluid()).otherPhase(phase_).U();
+    const volVectorField& Uc = tUc();
 
     const scalar sqrtPi = sqrt(constant::mathematical::pi);
     dimensionedScalar ThetaSmall("ThetaSmall", Theta_.dimensions(), 1e-6);
@@ -421,7 +424,7 @@ void Foam::RASModels::kineticTheoryModel::correct()
         volScalarField J2
         (
             "J2",
-            0.25*sqr(beta)*da*magSqr(U - Uc_)
+            0.25*sqr(beta)*da*magSqr(U - Uc)
            /(
                max(alpha, residualAlpha_)*rho
               *sqrtPi*(ThetaSqrt + ThetaSmallSqrt)
diff --git a/src/phaseSystemModels/reactingEuler/twoPhaseSystem/twoPhaseSystem.C b/src/phaseSystemModels/reactingEuler/twoPhaseSystem/twoPhaseSystem.C
index 51f8dd9898..d554ea70fe 100644
--- a/src/phaseSystemModels/reactingEuler/twoPhaseSystem/twoPhaseSystem.C
+++ b/src/phaseSystemModels/reactingEuler/twoPhaseSystem/twoPhaseSystem.C
@@ -136,8 +136,10 @@ void Foam::twoPhaseSystem::solve()
     word alphaScheme("div(phi," + alpha1.name() + ')');
     word alpharScheme("div(phir," + alpha1.name() + ')');
 
-    const surfaceScalarField& phi1 = phase1_.phi();
-    const surfaceScalarField& phi2 = phase2_.phi();
+    const tmp<surfaceScalarField> tphi1(phase1_.phi());
+    const surfaceScalarField& phi1 = tphi1();
+    const tmp<surfaceScalarField> tphi2(phase2_.phi());
+    const surfaceScalarField& phi2 = tphi2();
 
     // Construct the dilatation rate source term
     tmp<volScalarField::Internal> tdgdt;
diff --git a/src/phaseSystemModels/twoPhaseEuler/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/src/phaseSystemModels/twoPhaseEuler/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
index c621651753..ccb2ff83fc 100644
--- a/src/phaseSystemModels/twoPhaseEuler/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+++ b/src/phaseSystemModels/twoPhaseEuler/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
@@ -154,7 +154,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu
     }
 
     const fvPatchList& patches = phase.mesh().boundary();
-    const volVectorField& U = phase.U();
+    const tmp<volVectorField> tU(phase.U());
+    const volVectorField& U = tU();
 
     volScalarField::Boundary& nufBf = nuf.boundaryFieldRef();
 
diff --git a/src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C b/src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C
index 36d077fae2..f8129e55ed 100644
--- a/src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C
+++ b/src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C
@@ -80,11 +80,12 @@ Foam::diameterModels::IATEsources::turbulentBreakUp::R() const
 
     volScalarField R = tR();
 
-    scalar Cti = Cti_.value();
-    scalar WeCr = WeCr_.value();
+    const scalar Cti = Cti_.value();
+    const scalar WeCr = WeCr_.value();
     volScalarField Ut(this->Ut());
     volScalarField We(this->We());
-    const volScalarField& d(iate_.d()());
+    const tmp<volScalarField> td(iate_.d());
+    const volScalarField& d = td();
 
     forAll(R, celli)
     {
diff --git a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.C b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.C
index 47642607c3..f7a080808d 100644
--- a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.C
+++ b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.C
@@ -192,7 +192,7 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::updateCoeffs()
     }
 
     const scalarField UAvePyr(-phiPyr/patch().magSf());
-    const vectorField& nf = patch().nf();
+    tmp<vectorField> nf(patch().nf());
 
 
     // Evaluate velocity
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.C b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.C
index b1fb235f02..8c02bb90bc 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.C
@@ -61,8 +61,8 @@ multiBandZoneAbsorptionEmission
     emiCoeffs_(maxBands_),
     nBands_(0),
     zoneAbsorptivity_(),
-    zoneEmisivity_(),
-    zoneCells_()
+    zoneEmissivity_(),
+    zoneIds_()
 {
     coeffsDict_.readEntry("absorptivity", absCoeffs_);
     coeffsDict_.readEntry("emissivity", emiCoeffs_);
@@ -71,11 +71,11 @@ multiBandZoneAbsorptionEmission
     const dictionary& zoneDict = coeffsDict_.subDict("zones");
 
     zoneDict.readEntry("absorptivity", zoneAbsorptivity_);
-    zoneDict.readEntry("emissivity", zoneEmisivity_);
+    zoneDict.readEntry("emissivity", zoneEmissivity_);
 
-    zoneCells_.setSize(zoneAbsorptivity_.size(), -1);
+    zoneIds_.resize(zoneAbsorptivity_.size(), -1);
 
-    label i = 0;
+    label numZones = 0;
     forAllConstIters(zoneAbsorptivity_, iter)
     {
         label zoneID = mesh.cellZones().findZoneID(iter.key());
@@ -86,9 +86,10 @@ multiBandZoneAbsorptionEmission
                 << "Valid cellZones are " << mesh.cellZones().names()
                 << exit(FatalError);
         }
-        zoneCells_[i++] = zoneID;
+        zoneIds_[numZones] = zoneID;
+        ++numZones;
     }
-
+    // zoneIds_.resize(numZones);
 }
 
 
@@ -124,24 +125,19 @@ Foam::radiation::multiBandZoneAbsorptionEmission::aCont
         )
     );
 
-    volScalarField& a = ta.ref();
+    scalarField& a = ta.ref().primitiveFieldRef();
 
-    forAll(zoneCells_, zonei)
+    for (const label zonei : zoneIds_)
     {
-        const cellZone& cZone = mesh().cellZones()[zoneCells_[zonei]];
+        const cellZone& zn = mesh().cellZones()[zonei];
+        const auto iter = zoneAbsorptivity_.cfind(zn.name());
 
-        tmp<volScalarField> tzoneAbs(a*0.0);
-        volScalarField& zoneAbs = tzoneAbs.ref();
-
-        const scalarList& abs = zoneAbsorptivity_.find(cZone.name())();
-
-        forAll(cZone, i)
+        if (iter.good())  // Check is redundant (cannot fail)
         {
-            label cellId = cZone[i];
-            zoneAbs[cellId] = abs[bandI] - absCoeffs_[bandI];
-        }
+            const scalarList& absorb = iter.val();
 
-        a += zoneAbs;
+            UIndirectList<scalar>(a, zn) = absorb[bandI];
+        }
     }
 
     return ta;
@@ -171,25 +167,20 @@ Foam::radiation::multiBandZoneAbsorptionEmission::eCont
         )
     );
 
-    volScalarField& e = te.ref();
+    scalarField& e = te.ref().primitiveFieldRef();
 
-    forAll(zoneCells_, zonei)
+    for (const label zonei : zoneIds_)
     {
-        const cellZone& cZone = mesh().cellZones()[zoneCells_[zonei]];
+        const cellZone& zn = mesh().cellZones()[zonei];
+        const auto iter = zoneEmissivity_.cfind(zn.name());
 
-        tmp<volScalarField> tzoneEm(e*0.0);
-        volScalarField& zoneEm = tzoneEm.ref();
-
-        const scalarList& emi = zoneEmisivity_.find(cZone.name())();
-
-        forAll(cZone, i)
+        if (iter.good())  // Check is redundant (cannot fail)
         {
-            label cellId = cZone[i];
-            zoneEm[cellId] = emi[bandI] - emiCoeffs_[bandI];
-        }
-        e += zoneEm;
-    }
+            const scalarList& emit = iter.val();
 
+            UIndirectList<scalar>(e, zn) = emit[bandI];
+        }
+    }
 
     return te;
 }
diff --git a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.H b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.H
index 35eb5e7a1b..be13939944 100644
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.H
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.H
@@ -86,11 +86,11 @@ private:
         //- Cell zones absorptivity
         HashTable<scalarList> zoneAbsorptivity_;
 
-        //- Cell zones emisivity
-        HashTable<scalarList> zoneEmisivity_;
+        //- Cell zones emissivity
+        HashTable<scalarList> zoneEmissivity_;
 
-        //- Cells for each zone
-        labelList zoneCells_;
+        //- The cellZones selected
+        labelList zoneIds_;
 
 
 
diff --git a/src/waveModels/waveModel/waveModel.C b/src/waveModels/waveModel/waveModel.C
index 21c9cf3592..94ac689823 100644
--- a/src/waveModels/waveModel/waveModel.C
+++ b/src/waveModels/waveModel/waveModel.C
@@ -90,8 +90,7 @@ void Foam::waveModel::initialiseGeometry()
     }
 
     // Local face centres
-    const vectorField& Cf = patch_.faceCentres();
-    const vectorField CfLocal(Rgl_ & Cf);
+    const vectorField CfLocal(Rgl_ & patch_.faceCentres());
     z_ = CfLocal.component(2);
 
     // Local face extents in z-direction