diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/LaunderSharmaKE/LaunderSharmaKE.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/LaunderSharmaKE/LaunderSharmaKE.C
index 25e33b30ed328da90a28dc687021b8d843067525..c378e30fb13337fb7b90f29e9d9cdb26cc5d6175 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/LaunderSharmaKE/LaunderSharmaKE.C
@@ -7,7 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2007-2019 PCOpt/NTUA
     Copyright (C) 2013-2019 FOSS GP
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,34 +54,12 @@ LaunderSharmaKE::LaunderSharmaKE
 :
     RASModelVariables(mesh, SolverControl)
 {
-    hasTMVar1_ = true;
-    TMVar1Ptr_.reset
-    (
-        new tmp<volScalarField>
-        (
-            mesh_.lookupObjectRef<volScalarField>("k")
-        )
-    );
     TMVar1BaseName_ = "k";
-
-    hasTMVar2_ = true;
-    TMVar2Ptr_.reset
-    (
-        new tmp<volScalarField>
-        (
-            mesh_.lookupObjectRef<volScalarField>("epsilon")
-        )
-    );
     TMVar2BaseName_ = "epsilon";
 
-    hasNut_ = true;
-    nutPtr_.reset
-    (
-        new tmp<volScalarField>
-        (
-            mesh_.lookupObjectRef<volScalarField>("nut")
-        )
-    );
+    TMVar1Ptr_.ref(mesh_.lookupObjectRef<volScalarField>(TMVar1BaseName_));
+    TMVar2Ptr_.ref(mesh_.lookupObjectRef<volScalarField>(TMVar2BaseName_));
+    nutPtr_.ref(mesh_.lookupObjectRef<volScalarField>(nutBaseName_));
 
     allocateInitValues();
     allocateMeanFields();
diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.C
index b460f784baf21ffe26a72658861a131d17bfbcc9..8e132401ee1e79e98da27e3007422de7a167e9e2 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.C
+++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.C
@@ -42,6 +42,7 @@ namespace incompressible
 defineTypeNameAndDebug(RASModelVariables, 0);
 defineRunTimeSelectionTable(RASModelVariables, dictionary);
 
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 void RASModelVariables::allocateInitValues()
@@ -49,36 +50,28 @@ void RASModelVariables::allocateInitValues()
     if (solverControl_.storeInitValues())
     {
         Info<< "Storing initial values of turbulence variables" << endl;
-        if (hasTMVar1_)
+
+        if (hasTMVar1())
         {
             TMVar1InitPtr_.reset
             (
-                new volScalarField
-                (
-                    TMVar1Inst().name()+"Init",TMVar1Inst()
-                )
+                new volScalarField(TMVar1Inst().name()+"Init", TMVar1Inst())
             );
         }
 
-        if (hasTMVar2_)
+        if (hasTMVar2())
         {
             TMVar2InitPtr_.reset
             (
-                new volScalarField
-                (
-                    TMVar2Inst().name()+"Init",TMVar2Inst()
-                )
+                new volScalarField(TMVar2Inst().name()+"Init", TMVar2Inst())
             );
         }
 
-        if (hasNut_)
+        if (hasNut())
         {
             nutInitPtr_.reset
             (
-                new volScalarField
-                (
-                    nutRefInst().name()+"Init",nutRefInst()
-                )
+                new volScalarField(nutRefInst().name()+"Init", nutRefInst())
             );
         }
     }
@@ -90,7 +83,8 @@ void RASModelVariables::allocateMeanFields()
     if (solverControl_.average())
     {
         Info<< "Allocating mean values of turbulence variables" << endl;
-        if (hasTMVar1_)
+
+        if (hasTMVar1())
         {
             TMVar1MeanPtr_.reset
             (
@@ -108,7 +102,8 @@ void RASModelVariables::allocateMeanFields()
                 )
             );
         }
-        if (hasTMVar2_)
+
+        if (hasTMVar2())
         {
             TMVar2MeanPtr_.reset
             (
@@ -127,7 +122,7 @@ void RASModelVariables::allocateMeanFields()
             );
         }
 
-        if (hasNut_)
+        if (hasNut())
         {
             nutMeanPtr_.reset
             (
@@ -149,25 +144,19 @@ void RASModelVariables::allocateMeanFields()
 }
 
 
-RASModelVariables::autoTmp
-RASModelVariables::cloneAutoTmp(const autoTmp& source) const
+Foam::refPtr<Foam::volScalarField>
+RASModelVariables::cloneRefPtr(const refPtr<volScalarField>& obj) const
 {
-    autoTmp returnField(nullptr);
-    if (source && source->valid())
+    if (obj)
     {
-        const volScalarField& sf = source()();
-        DebugInfo
-            << "Cloning " << sf.name() << endl;
+        const volScalarField& sf = obj();
+
         const word timeName = mesh_.time().timeName();
-        returnField.reset
-        (
-            new tmp<volScalarField>
-            (
-                new volScalarField(sf.name() + timeName, sf)
-            )
-        );
+
+        return refPtr<volScalarField>::New(sf.name() + timeName, sf);
     }
-    return returnField;
+
+    return nullptr;
 }
 
 
@@ -198,20 +187,20 @@ RASModelVariables::RASModelVariables
 :
     mesh_(mesh),
     solverControl_(SolverControl),
-    hasTMVar1_(false),
-    hasTMVar2_(false),
-    hasNut_(false),
-    hasDist_(false),
+
+    TMVar1BaseName_(),
+    TMVar2BaseName_(),
+    nutBaseName_("nut"),
+
     TMVar1Ptr_(nullptr),
     TMVar2Ptr_(nullptr),
     nutPtr_(nullptr),
-    dPtr_(nullptr),
-    TMVar1BaseName_(word::null),
-    TMVar2BaseName_(word::null),
-    nutBaseName_("nut"),
+    distPtr_(nullptr),
+
     TMVar1InitPtr_(nullptr),
     TMVar2InitPtr_(nullptr),
     nutInitPtr_(nullptr),
+
     TMVar1MeanPtr_(nullptr),
     TMVar2MeanPtr_(nullptr),
     nutMeanPtr_(nullptr)
@@ -225,25 +214,24 @@ RASModelVariables::RASModelVariables
 :
     mesh_(rmv.mesh_),
     solverControl_(rmv.solverControl_),
-    hasTMVar1_(rmv.hasTMVar1_),
-    hasTMVar2_(rmv.hasTMVar2_),
-    hasNut_(rmv.hasNut_),
-    hasDist_(rmv.hasDist_),
-    TMVar1Ptr_(cloneAutoTmp(rmv.TMVar1Ptr_)),
-    TMVar2Ptr_(cloneAutoTmp(rmv.TMVar2Ptr_)),
-    nutPtr_(cloneAutoTmp(rmv.nutPtr_)),
-    dPtr_(cloneAutoTmp(rmv.dPtr_)),
+
     TMVar1BaseName_(rmv.TMVar1BaseName_),
     TMVar2BaseName_(rmv.TMVar2BaseName_),
     nutBaseName_(rmv.nutBaseName_),
+
+    TMVar1Ptr_(cloneRefPtr(rmv.TMVar1Ptr_)),
+    TMVar2Ptr_(cloneRefPtr(rmv.TMVar2Ptr_)),
+    nutPtr_(cloneRefPtr(rmv.nutPtr_)),
+    distPtr_(cloneRefPtr(rmv.distPtr_)),
+
     TMVar1InitPtr_(nullptr),
     TMVar2InitPtr_(nullptr),
     nutInitPtr_(nullptr),
+
     TMVar1MeanPtr_(nullptr),
     TMVar2MeanPtr_(nullptr),
     nutMeanPtr_(nullptr)
-{
-}
+{}
 
 
 autoPtr<RASModelVariables> RASModelVariables::clone() const
@@ -251,6 +239,7 @@ autoPtr<RASModelVariables> RASModelVariables::clone() const
     return autoPtr<RASModelVariables>::New(*this);
 }
 
+
 // * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
 
 autoPtr<RASModelVariables> RASModelVariables::New
@@ -307,172 +296,6 @@ autoPtr<RASModelVariables> RASModelVariables::New
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool RASModelVariables::hasTMVar1() const
-{
-    return hasTMVar1_;
-}
-
-
-bool RASModelVariables::hasTMVar2() const
-{
-    return hasTMVar2_;
-}
-
-
-bool RASModelVariables::hasNut() const
-{
-    return hasNut_;
-}
-
-
-bool RASModelVariables::hasDist() const
-{
-    return hasDist_;
-}
-
-
-const word& RASModelVariables::TMVar1BaseName() const
-{
-    return TMVar1BaseName_;
-}
-
-
-const word& RASModelVariables::TMVar2BaseName() const
-{
-    return TMVar2BaseName_;
-}
-
-
-const word& RASModelVariables::nutBaseName() const
-{
-    return nutBaseName_;
-}
-
-
-const volScalarField& RASModelVariables::TMVar1() const
-{
-    if (solverControl_.useAveragedFields())
-    {
-        return TMVar1MeanPtr_();
-    }
-    else
-    {
-        return TMVar1Ptr_()();
-    }
-}
-
-
-volScalarField& RASModelVariables::TMVar1()
-{
-    if (solverControl_.useAveragedFields())
-    {
-        return TMVar1MeanPtr_();
-    }
-    else
-    {
-        return TMVar1Ptr_().constCast();
-    }
-}
-
-
-const volScalarField& RASModelVariables::TMVar2() const
-{
-    if (solverControl_.useAveragedFields())
-    {
-        return TMVar2MeanPtr_();
-    }
-    else
-    {
-        return TMVar2Ptr_()();
-    }
-}
-
-volScalarField& RASModelVariables::TMVar2()
-{
-    if (solverControl_.useAveragedFields())
-    {
-        return TMVar2MeanPtr_();
-    }
-    else
-    {
-        return TMVar2Ptr_().constCast();
-    }
-}
-
-const volScalarField& RASModelVariables::nutRef() const
-{
-    if (solverControl_.useAveragedFields() && hasNut_)
-    {
-        return nutMeanPtr_();
-    }
-    else
-    {
-        return nutPtr_()();
-    }
-}
-
-
-volScalarField& RASModelVariables::nutRef()
-{
-    if (solverControl_.useAveragedFields() && hasNut_)
-    {
-        return  nutMeanPtr_();
-    }
-    else
-    {
-        return nutPtr_().constCast();
-    }
-}
-
-
-const volScalarField& RASModelVariables::d() const
-{
-    return dPtr_()();
-}
-
-
-volScalarField& RASModelVariables::d()
-{
-    return dPtr_().constCast();
-}
-
-
-const volScalarField& RASModelVariables::TMVar1Inst() const
-{
-    return TMVar1Ptr_()();
-}
-
-
-volScalarField& RASModelVariables::TMVar1Inst()
-{
-    return TMVar1Ptr_().constCast();
-}
-
-
-const volScalarField& RASModelVariables::TMVar2Inst() const
-{
-    return TMVar2Ptr_()();
-}
-
-
-volScalarField& RASModelVariables::TMVar2Inst()
-{
-    return TMVar2Ptr_().constCast();
-}
-
-
-const volScalarField& RASModelVariables::nutRefInst() const
-{
-    return nutPtr_()();
-}
-
-
-volScalarField& RASModelVariables::nutRefInst()
-{
-    return nutPtr_().constCast();
-}
-
-
 tmp<volScalarField> RASModelVariables::nutJacobianVar1
 (
     const singlePhaseTransportModel& laminarTransport
@@ -482,23 +305,19 @@ tmp<volScalarField> RASModelVariables::nutJacobianVar1
         << "jutJacobianVar1 not implemented for the current turbulence model."
         << "Returning zero field" << endl;
 
-    tmp<volScalarField> nutJacobian
+    return tmp<volScalarField>::New
     (
-        new volScalarField
+        IOobject
         (
-            IOobject
-            (
-                "nutJacobianVar1",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            "nutJacobianVar1",
+            mesh_.time().timeName(),
             mesh_,
-            dimensionedScalar(dimless, Zero)
-        )
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar(dimless, Zero)
     );
-    return nutJacobian;
 }
 
 
@@ -511,38 +330,35 @@ tmp<volScalarField> RASModelVariables::nutJacobianVar2
         << "nutJacobianVar2 not implemented for the current turbulence model."
         << "Returning zero field" << endl;
 
-    tmp<volScalarField> nutJacobian
+    return tmp<volScalarField>::New
     (
-        new volScalarField
+        IOobject
         (
-            IOobject
-            (
-                "nutJacobianVar2",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            "nutJacobianVar2",
+            mesh_.time().timeName(),
             mesh_,
-            dimensionedScalar(dimless, Zero)
-        )
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar(dimless, Zero)
     );
-    return nutJacobian;
 }
 
+
 void RASModelVariables::restoreInitValues()
 {
     if (solverControl_.storeInitValues())
     {
-        if (hasTMVar1_)
+        if (hasTMVar1())
         {
             TMVar1Inst() == TMVar1InitPtr_();
         }
-        if (hasTMVar2_)
+        if (hasTMVar2())
         {
             TMVar2Inst() == TMVar2InitPtr_();
         }
-        if (hasNut_)
+        if (hasNut())
         {
             nutRefInst() == nutInitPtr_();
         }
@@ -557,19 +373,20 @@ void RASModelVariables::resetMeanFields()
         Info<< "Resetting mean turbulent fields to zero" << endl;
 
         // Reset fields to zero
-        if (hasTMVar1_)
+        if (TMVar1Ptr_)
         {
-            TMVar1MeanPtr_() ==
+            TMVar1MeanPtr_.ref() ==
                 dimensionedScalar(TMVar1Inst().dimensions(), Zero);
         }
-        if (hasTMVar2_)
+        if (TMVar2Ptr_)
         {
-            TMVar2MeanPtr_() ==
+            TMVar2MeanPtr_.ref() ==
                 dimensionedScalar(TMVar2Inst().dimensions(), Zero);
         }
-        if (hasNut_)
+        if (nutPtr_)
         {
-            nutMeanPtr_() == dimensionedScalar(nutRefInst().dimensions(), Zero);
+            nutMeanPtr_.ref() ==
+                dimensionedScalar(nutRefInst().dimensions(), Zero);
         }
     }
 }
@@ -580,22 +397,24 @@ void RASModelVariables::computeMeanFields()
     if (solverControl_.doAverageIter())
     {
         const label iAverageIter = solverControl_.averageIter();
-        scalar avIter(iAverageIter);
-        scalar oneOverItP1 = 1./(avIter + 1);
-        scalar mult = avIter*oneOverItP1;
-        if (hasTMVar1_)
+        const scalar avIter(iAverageIter);
+        const scalar oneOverItP1 = 1./(avIter + 1);
+        const scalar mult = avIter*oneOverItP1;
+
+        if (hasTMVar1())
         {
-            TMVar1MeanPtr_() ==
-                TMVar1MeanPtr_()*mult + TMVar1Inst()*oneOverItP1;
+            TMVar1MeanPtr_.ref() ==
+                (TMVar1MeanPtr_()*mult + TMVar1Inst()*oneOverItP1);
         }
-        if (hasTMVar2_)
+        if (hasTMVar2())
         {
-            TMVar2MeanPtr_() ==
-                TMVar2MeanPtr_()*mult + TMVar2Inst()*oneOverItP1;
+            TMVar2MeanPtr_.ref() ==
+                (TMVar2MeanPtr_()*mult + TMVar2Inst()*oneOverItP1);
         }
-        if (hasNut_)
+        if (hasNut())
         {
-            nutMeanPtr_() == nutMeanPtr_()*mult + nutRefInst()*oneOverItP1;
+            nutMeanPtr_.ref() ==
+                (nutMeanPtr_()*mult + nutRefInst()*oneOverItP1);
         }
     }
 }
@@ -607,20 +426,17 @@ tmp<volSymmTensorField> RASModelVariables::devReff
     const volVectorField& U
 ) const
 {
-    return tmp<volSymmTensorField>
+    return tmp<volSymmTensorField>::New
     (
-        new volSymmTensorField
+        IOobject
         (
-            IOobject
-            (
-                "devRhoReff",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-           -(laminarTransport.nu() + nutRef())*dev(twoSymm(fvc::grad(U)))
-        )
+            "devRhoReff",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        -(laminarTransport.nu() + nutRef())*dev(twoSymm(fvc::grad(U)))
     );
 }
 
@@ -632,28 +448,28 @@ void RASModelVariables::correctBoundaryConditions
 {
     if (hasTMVar1())
     {
-        TMVar1Ptr_().constCast().correctBoundaryConditions();
+        TMVar1Inst().correctBoundaryConditions();
         if (solverControl_.average())
         {
-            TMVar1MeanPtr_().correctBoundaryConditions();
+            TMVar1MeanPtr_.ref().correctBoundaryConditions();
         }
     }
 
     if (hasTMVar2())
     {
-        TMVar2Ptr_().constCast().correctBoundaryConditions();
+        TMVar2Inst().correctBoundaryConditions();
         if (solverControl_.average())
         {
-            TMVar2MeanPtr_().correctBoundaryConditions();
+            TMVar2MeanPtr_.ref().correctBoundaryConditions();
         }
     }
 
     if (hasNut())
     {
-        nutPtr_().constCast().correctBoundaryConditions();
+        nutRefInst().correctBoundaryConditions();
         if (solverControl_.average())
         {
-            nutMeanPtr_().correctBoundaryConditions();
+            nutMeanPtr_.ref().correctBoundaryConditions();
         }
     }
 }
@@ -661,22 +477,22 @@ void RASModelVariables::correctBoundaryConditions
 
 void RASModelVariables::transfer(RASModelVariables& rmv)
 {
-    if (rmv.hasTMVar1() && hasTMVar1_)
+    if (rmv.hasTMVar1() && hasTMVar1())
     {
         copyAndRename(TMVar1Inst(), rmv.TMVar1Inst());
     }
 
-    if (rmv.hasTMVar2() && hasTMVar2_)
+    if (rmv.hasTMVar2() && hasTMVar2())
     {
         copyAndRename(TMVar2Inst(), rmv.TMVar2Inst());
     }
 
-    if (rmv.hasNut() && hasNut_)
+    if (rmv.hasNut() && hasNut())
     {
-        copyAndRename(nutRef(), rmv.nutRef());
+        copyAndRename(nutRefInst(), rmv.nutRefInst());
     }
 
-    if (rmv.hasDist() && hasDist_)
+    if (rmv.hasDist() && hasDist())
     {
         copyAndRename(d(), rmv.d());
     }
diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.H
index dae7324f975e1cf91f263c5a4d8bb331f77efa9c..09800b9cb92c9d41e8191b4f88b04620e2a042d6 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.H
+++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.H
@@ -7,7 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2007-2019 PCOpt/NTUA
     Copyright (C) 2013-2019 FOSS GP
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,6 +44,7 @@ SourceFiles
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
 #include "runTimeSelectionTables.H"
+#include "refPtr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,61 +59,44 @@ namespace incompressible
 
 class RASModelVariables
 {
-
-public:
-
-    // Public typedefs
-
-        typedef autoPtr<tmp<volScalarField>> autoTmp;
-
 protected:
 
-    // Protected data
-
+    // Protected Data
 
         const fvMesh& mesh_;
         const solverControl& solverControl_;
 
-        // autoPtrs delete the memory on destruction
-        // Can cause memory mishandling issues in this case
-        // Use regular ptrs instead
-        bool hasTMVar1_;
-        bool hasTMVar2_;
-        bool hasNut_;
-        bool hasDist_;
-        autoTmp TMVar1Ptr_;
-        autoTmp TMVar2Ptr_;
-        autoTmp nutPtr_;
-        autoTmp dPtr_;
         // Base names of the turbulent fields
         word TMVar1BaseName_;
         word TMVar2BaseName_;
         word nutBaseName_;
 
-        // conditionally store initial values
+        refPtr<volScalarField> TMVar1Ptr_;
+        refPtr<volScalarField> TMVar2Ptr_;
+        refPtr<volScalarField> nutPtr_;
+        refPtr<volScalarField> distPtr_;
+
+        // Initial values (demand-driven)
         // For finite differences and optimisation runs
-        autoPtr<volScalarField> TMVar1InitPtr_;
-        autoPtr<volScalarField> TMVar2InitPtr_;
-        autoPtr<volScalarField> nutInitPtr_;
+        refPtr<volScalarField> TMVar1InitPtr_;
+        refPtr<volScalarField> TMVar2InitPtr_;
+        refPtr<volScalarField> nutInitPtr_;
 
-        // conditionally store mean values
-        autoPtr<volScalarField> TMVar1MeanPtr_;
-        autoPtr<volScalarField> TMVar2MeanPtr_;
-        autoPtr<volScalarField> nutMeanPtr_;
+        // Mean values (demand-driven)
+        refPtr<volScalarField> TMVar1MeanPtr_;
+        refPtr<volScalarField> TMVar2MeanPtr_;
+        refPtr<volScalarField> nutMeanPtr_;
 
 
     // Protected Member functions
 
         void allocateInitValues();
         void allocateMeanFields();
-        autoTmp cloneAutoTmp(const autoTmp& source) const;
-        void copyAndRename(volScalarField& f1, volScalarField& f2);
-
-
-private:
 
-    // Private Member Functions
+        refPtr<volScalarField>
+        cloneRefPtr(const refPtr<volScalarField>& obj) const;
 
+        void copyAndRename(volScalarField& f1, volScalarField& f2);
 
         //- No copy assignment
         void operator=(const RASModelVariables&) = delete;
@@ -186,43 +170,44 @@ public:
 
     // Member Functions
 
+        //- Turbulence field names
+        inline const word& TMVar1BaseName() const;
+        inline const word& TMVar2BaseName() const;
+        inline const word& nutBaseName() const;
 
         //- Bools to idenify which turbulent fields are present
-        bool hasTMVar1() const;
-        bool hasTMVar2() const;
-        bool hasNut() const;
-        bool hasDist() const;
-
-        //- Turbulence field names
-        const word& TMVar1BaseName() const;
-        const word& TMVar2BaseName() const;
-        const word& nutBaseName() const;
+        inline bool hasTMVar1() const;
+        inline bool hasTMVar2() const;
+        inline bool hasNut() const;
+        inline bool hasDist() const;
 
         //- Return references to turbulence fields
         //  will return the mean field if it exists, otherwise the
         //  instantaneous one
-        const volScalarField& TMVar1() const;
-              volScalarField& TMVar1();
-        const volScalarField& TMVar2() const;
-              volScalarField& TMVar2();
-        const volScalarField& nutRef() const;
-              volScalarField& nutRef();
-        const volScalarField& d() const;
-              volScalarField& d();
+        inline const volScalarField& TMVar1() const;
+        inline       volScalarField& TMVar1();
+        inline const volScalarField& TMVar2() const;
+        inline       volScalarField& TMVar2();
+        inline const volScalarField& nutRef() const;
+        inline       volScalarField& nutRef();
+        inline const volScalarField& d() const;
+        inline       volScalarField& d();
 
         //- return references to instantaneous turbulence fields
-        const volScalarField& TMVar1Inst() const;
-              volScalarField& TMVar1Inst();
-        const volScalarField& TMVar2Inst() const;
-              volScalarField& TMVar2Inst();
-        const volScalarField& nutRefInst() const;
-              volScalarField& nutRefInst();
+        inline const volScalarField& TMVar1Inst() const;
+        inline       volScalarField& TMVar1Inst();
+        inline const volScalarField& TMVar2Inst() const;
+        inline       volScalarField& TMVar2Inst();
+        inline const volScalarField& nutRefInst() const;
+        inline       volScalarField& nutRefInst();
+
 
         //- Return nut Jacobian wrt the TM vars
         virtual tmp<volScalarField> nutJacobianVar1
         (
             const singlePhaseTransportModel& laminarTransport
         ) const;
+
         virtual tmp<volScalarField> nutJacobianVar2
         (
             const singlePhaseTransportModel& laminarTransport
@@ -264,6 +249,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "RASModelVariablesI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariablesI.H b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariablesI.H
new file mode 100644
index 0000000000000000000000000000000000000000..784f0ef2c887cb8491d59cfda3572ecb48278f93
--- /dev/null
+++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariablesI.H
@@ -0,0 +1,197 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2007-2019 PCOpt/NTUA
+    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const word& RASModelVariables::TMVar1BaseName() const
+{
+    return TMVar1BaseName_;
+}
+
+
+inline const word& RASModelVariables::TMVar2BaseName() const
+{
+    return TMVar2BaseName_;
+}
+
+
+inline const word& RASModelVariables::nutBaseName() const
+{
+    return nutBaseName_;
+}
+
+
+inline bool RASModelVariables::hasTMVar1() const
+{
+    return bool(TMVar1Ptr_);
+}
+
+
+inline bool RASModelVariables::hasTMVar2() const
+{
+    return bool(TMVar2Ptr_);
+}
+
+
+inline bool RASModelVariables::hasNut() const
+{
+    return bool(nutPtr_);
+}
+
+
+inline bool RASModelVariables::hasDist() const
+{
+    return bool(distPtr_);
+}
+
+
+inline const volScalarField& RASModelVariables::TMVar1() const
+{
+    if (solverControl_.useAveragedFields())
+    {
+        return TMVar1MeanPtr_.cref();
+    }
+
+    return TMVar1Ptr_.cref();
+}
+
+
+inline volScalarField& RASModelVariables::TMVar1()
+{
+    if (solverControl_.useAveragedFields())
+    {
+        return TMVar1MeanPtr_.ref();
+    }
+
+    return TMVar1Ptr_.ref();
+}
+
+
+inline const volScalarField& RASModelVariables::TMVar2() const
+{
+    if (solverControl_.useAveragedFields())
+    {
+        return TMVar2MeanPtr_.cref();
+    }
+
+    return TMVar2Ptr_.cref();
+}
+
+inline volScalarField& RASModelVariables::TMVar2()
+{
+    if (solverControl_.useAveragedFields())
+    {
+        return TMVar2MeanPtr_.ref();
+    }
+
+    return TMVar2Ptr_.ref();
+}
+
+inline const volScalarField& RASModelVariables::nutRef() const
+{
+    if (solverControl_.useAveragedFields() && nutMeanPtr_)
+    {
+        return nutMeanPtr_.cref();
+    }
+
+    return nutPtr_.cref();
+}
+
+
+inline volScalarField& RASModelVariables::nutRef()
+{
+    if (solverControl_.useAveragedFields() && nutMeanPtr_)
+    {
+        return nutMeanPtr_.ref();
+    }
+
+    return nutPtr_.ref();
+}
+
+
+inline const volScalarField& RASModelVariables::d() const
+{
+    return distPtr_.cref();
+}
+
+
+inline volScalarField& RASModelVariables::d()
+{
+    return distPtr_.ref();
+}
+
+
+inline const volScalarField& RASModelVariables::TMVar1Inst() const
+{
+    return TMVar1Ptr_.cref();
+}
+
+
+inline volScalarField& RASModelVariables::TMVar1Inst()
+{
+    return TMVar1Ptr_.ref();
+}
+
+
+inline const volScalarField& RASModelVariables::TMVar2Inst() const
+{
+    return TMVar2Ptr_.cref();
+}
+
+
+inline volScalarField& RASModelVariables::TMVar2Inst()
+{
+    return TMVar2Ptr_.ref();
+}
+
+
+inline const volScalarField& RASModelVariables::nutRefInst() const
+{
+    return nutPtr_.cref();
+}
+
+
+inline volScalarField& RASModelVariables::nutRefInst()
+{
+    return nutPtr_.ref();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //
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 40b9f9110e3b9f0aa871f1990dbb0908f65db969..7ffd13e1e11bae1d291d2c3eceac140fe12f3f1b 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/SpalartAllmaras/SpalartAllmaras.C
@@ -7,7 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2007-2019 PCOpt/NTUA
     Copyright (C) 2013-2019 FOSS GP
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,54 +55,35 @@ SpalartAllmaras::SpalartAllmaras
 :
     RASModelVariables(mesh, SolverControl)
 {
-    hasTMVar1_ = true;
-    TMVar1Ptr_.reset
-    (
-        new tmp<volScalarField>
-        (
-            mesh_.lookupObjectRef<volScalarField>("nuTilda")
-        )
-    );
     TMVar1BaseName_ = "nuTilda";
 
+    TMVar1Ptr_.ref(mesh_.lookupObjectRef<volScalarField>(TMVar1BaseName_));
+
     TMVar2Ptr_.reset
     (
-        new tmp<volScalarField>
+        new volScalarField
         (
-            new volScalarField
+            IOobject
             (
-                IOobject
-                (
-                    "dummySpalartAllmarasVar2",
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
+                "dummySpalartAllmarasVar2",
+                mesh.time().timeName(),
                 mesh,
-                dimensionedScalar(dimless, Zero)
-            )
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar(dimless, Zero)
         )
     );
 
-    hasNut_ = true;
-    nutPtr_.reset
-    (
-        new tmp<volScalarField>
-        (
-            mesh_.lookupObjectRef<volScalarField>("nut")
-        )
-    );
+    nutPtr_.ref(mesh_.lookupObjectRef<volScalarField>(nutBaseName_));
 
-    hasDist_ = true;
-    dPtr_.reset
+    // The wall dist name can vary depending on how wallDist was
+    // constructed. Grab the field directly from wallDist
+
+    distPtr_.ref
     (
-        new tmp<volScalarField>
-        (
-            // The wall dist name can vary depending on how wallDist was
-            // constructed. Grab the field directly from wallDist
-            const_cast<volScalarField&>(wallDist::New(mesh_).y())
-        )
+        const_cast<volScalarField&>(wallDist::New(mesh_).y())
     );
 
     allocateInitValues();
@@ -117,32 +98,32 @@ tmp<volScalarField> SpalartAllmaras::nutJacobianVar1
     const singlePhaseTransportModel& laminarTransport
 ) const
 {
-    tmp<volScalarField> tnutJacobian
+    auto tnutJacobian = tmp<volScalarField>::New
     (
-        new volScalarField
+        IOobject
         (
-            IOobject
-            (
-                "nutJacobianVar1",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            "nutJacobianVar1",
+            mesh_.time().timeName(),
             mesh_,
-            dimensionedScalar(dimless, Zero)
-        )
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar(dimless, Zero)
     );
-    volScalarField& nutJacobian = tnutJacobian.ref();
+    auto& nutJacobian = tnutJacobian.ref();
 
     const volScalarField& nu = laminarTransport.nu();
     const volScalarField& nuTilda = TMVar1();
+
     volScalarField chi(nuTilda/nu);
     volScalarField chi3(pow3(chi));
-    scalar Cv13 = pow3(7.1);
+
+    const scalar Cv13 = pow3(7.1);
     volScalarField fv1(chi3/(chi3 + Cv13));
     volScalarField fv1ByChi2Sqr(sqr(chi/(chi3 + Cv13)));
     volScalarField Cdfv1(3.0*Cv13*fv1ByChi2Sqr);
+
     nutJacobian = Cdfv1*chi + fv1;
 
     return tnutJacobian;
diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kEpsilon/kEpsilon.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kEpsilon/kEpsilon.C
index 1aefa84200ada0c8d39b369953f437487cc119c1..a05a88cd4937cd9dde8dcd770ef760d2bae4cd7b 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kEpsilon/kEpsilon.C
+++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kEpsilon/kEpsilon.C
@@ -7,7 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2007-2019 PCOpt/NTUA
     Copyright (C) 2013-2019 FOSS GP
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,34 +54,12 @@ kEpsilon::kEpsilon
 :
     RASModelVariables(mesh, SolverControl)
 {
-    hasTMVar1_ = true;
-    TMVar1Ptr_.reset
-    (
-        new tmp<volScalarField>
-        (
-            mesh_.lookupObjectRef<volScalarField>("k")
-        )
-    );
     TMVar1BaseName_ = "k";
-
-    hasTMVar2_ = true;
-    TMVar2Ptr_.reset
-    (
-        new tmp<volScalarField>
-        (
-            mesh_.lookupObjectRef<volScalarField>("epsilon")
-        )
-    );
     TMVar2BaseName_ = "epsilon";
 
-    hasNut_ = true;
-    nutPtr_.reset
-    (
-        new tmp<volScalarField>
-        (
-            mesh_.lookupObjectRef<volScalarField>("nut")
-        )
-    );
+    TMVar1Ptr_.ref(mesh_.lookupObjectRef<volScalarField>(TMVar1BaseName_));
+    TMVar2Ptr_.ref(mesh_.lookupObjectRef<volScalarField>(TMVar2BaseName_));
+    nutPtr_.ref(mesh_.lookupObjectRef<volScalarField>(nutBaseName_));
 
     allocateInitValues();
     allocateMeanFields();
diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.C
index 8b96bd5d5e69afe2447863f4966e59680b595c60..da8eed8cc13714cbadb0040ea8d41ddb977ca14c 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/kOmegaSST/kOmegaSST.C
@@ -7,7 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2007-2019 PCOpt/NTUA
     Copyright (C) 2013-2019 FOSS GP
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,34 +54,12 @@ kOmegaSST::kOmegaSST
 :
     RASModelVariables(mesh, SolverControl)
 {
-    hasTMVar1_ = true;
-    TMVar1Ptr_.reset
-    (
-        new tmp<volScalarField>
-        (
-            mesh_.lookupObjectRef<volScalarField>("k")
-        )
-    );
     TMVar1BaseName_ = "k";
-
-    hasTMVar2_ = true;
-    TMVar2Ptr_.reset
-    (
-        new tmp<volScalarField>
-        (
-            mesh_.lookupObjectRef<volScalarField>("omega")
-        )
-    );
     TMVar2BaseName_ = "omega";
 
-    hasNut_ = true;
-    nutPtr_.reset
-    (
-        new tmp<volScalarField>
-        (
-            mesh_.lookupObjectRef<volScalarField>("nut")
-        )
-    );
+    TMVar1Ptr_.ref(mesh_.lookupObjectRef<volScalarField>(TMVar1BaseName_));
+    TMVar2Ptr_.ref(mesh_.lookupObjectRef<volScalarField>(TMVar2BaseName_));
+    nutPtr_.ref(mesh_.lookupObjectRef<volScalarField>(nutBaseName_));
 
     allocateInitValues();
     allocateMeanFields();
@@ -96,8 +74,9 @@ void kOmegaSST::correctBoundaryConditions
 )
 {
     // The presence of G is required to update the boundary value of omega
-    const volVectorField& U(turbulence.U());
+    const volVectorField& U = turbulence.U();
     const volScalarField S2(2*magSqr(symm(fvc::grad(U))));
+
     volScalarField G(turbulence.GName(), nutRef() * S2);
     RASModelVariables::correctBoundaryConditions(turbulence);
 }
diff --git a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/laminar/laminar.C b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/laminar/laminar.C
index 091f4fbd7dba23b26605c37bddcd8974cbe0a033..3560e42a166ba312e6bf394b3e6b6cc850ded557 100644
--- a/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/laminar/laminar.C
+++ b/src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/laminar/laminar.C
@@ -7,7 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2007-2019 PCOpt/NTUA
     Copyright (C) 2013-2019 FOSS GP
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,63 +56,55 @@ laminar::laminar
 {
     TMVar1Ptr_.reset
     (
-        new tmp<volScalarField>
+        new volScalarField
         (
-            new volScalarField
+            IOobject
             (
-                IOobject
-                (
-                    "dummylaminarVar1",
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
+                "dummylaminarVar1",
+                mesh.time().timeName(),
                 mesh,
-                dimensionedScalar(dimless, Zero)
-            )
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar(dimless, Zero)
         )
     );
 
     TMVar2Ptr_.reset
     (
-        new tmp<volScalarField>
+        new volScalarField
         (
-            new volScalarField
+            IOobject
             (
-                IOobject
-                (
-                    "dummylaminarVar2",
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
+                "dummylaminarVar2",
+                mesh.time().timeName(),
                 mesh,
-                dimensionedScalar(dimless, Zero)
-            )
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar(dimless, Zero)
         )
     );
 
     nutPtr_.reset
     (
-        new tmp<volScalarField>
+        new volScalarField
         (
-            new volScalarField
+            IOobject
             (
-                IOobject
-                (
-                    "dummylaminarNut",
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
+                "dummylaminarNut",
+                mesh.time().timeName(),
                 mesh,
-                dimensionedScalar(sqr(dimLength)/dimTime, Zero)
-            )
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar(sqr(dimLength)/dimTime, Zero)
         )
     );
+
     allocateInitValues();
 }