diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.C b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
index b98f60e66a2df840e92cce0a12121d73655f8634..6667bfcefaede11601e54747727a493c2f6d7ed9 100644
--- a/src/finiteArea/faMatrices/faMatrix/faMatrix.C
+++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
@@ -49,7 +49,7 @@ void Foam::faMatrix<Type>::addToInternalField
     {
         FatalErrorInFunction
             << "addressing (" << addr.size()
-            << ") and field (" << pf.size() << ") are different sizes" << endl
+            << ") and field (" << pf.size() << ") are different sizes"
             << abort(FatalError);
     }
 
@@ -87,7 +87,7 @@ void Foam::faMatrix<Type>::subtractFromInternalField
     {
         FatalErrorInFunction
             << "addressing (" << addr.size()
-            << ") and field (" << pf.size() << ") are different sizes" << endl
+            << ") and field (" << pf.size() << ") are different sizes"
             << abort(FatalError);
     }
 
diff --git a/src/finiteArea/faMesh/faMeshDemandDrivenData.C b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
index 6215a62bc2ef4bfdffb3f2dbcf9a8c8a71d36347..7d63495708101f71a4acec4469b51c3c0c8a5c07 100644
--- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C
+++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
@@ -1852,14 +1852,14 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
     }
 
 
-    forAll(boundary(), patchI)
+    forAll(boundary(), patchi)
     {
-        const faPatch& fap = boundary()[patchI];
+        const faPatch& fap = boundary()[patchi];
 
         if (Pstream::parRun() && isA<processorFaPatch>(fap))
         {
             const processorFaPatch& procPatch =
-                refCast<const processorFaPatch>(boundary()[patchI]);
+                refCast<const processorFaPatch>(boundary()[patchi]);
 
             const labelList& patchPointLabels = procPatch.pointLabels();
 
diff --git a/src/finiteArea/finiteArea/convectionSchemes/faConvectionScheme/faConvectionScheme.H b/src/finiteArea/finiteArea/convectionSchemes/faConvectionScheme/faConvectionScheme.H
index ae2673c8fc6a9d1e9382d42b94747f9105091436..5170627186aadc33d4c076c8ef645d767e5d9000 100644
--- a/src/finiteArea/finiteArea/convectionSchemes/faConvectionScheme/faConvectionScheme.H
+++ b/src/finiteArea/finiteArea/convectionSchemes/faConvectionScheme/faConvectionScheme.H
@@ -48,6 +48,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 template<class Type>
 class faMatrix;
 
@@ -67,8 +68,9 @@ class convectionScheme
 :
     public refCount
 {
-    // Private data
+    // Private Data
 
+        //- Reference to the mesh
         const faMesh& mesh_;
 
 
@@ -129,11 +131,8 @@ public:
 
     // Member Functions
 
-        //- Return mesh reference
-        const faMesh& mesh() const
-        {
-            return mesh_;
-        }
+        //- Return const reference to mesh
+        const faMesh& mesh() const noexcept { return mesh_; }
 
         virtual tmp<GeometricField<Type, faePatchField, edgeMesh>> flux
         (
diff --git a/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.C b/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.C
index d953d01939ec135df5a1e2f2ba4358c692789bc3..337a4b49178bce3fd765c16af59cf22a870584a6 100644
--- a/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.C
+++ b/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.C
@@ -64,13 +64,10 @@ gaussConvectionScheme<Type>::famDiv
     tmp<edgeScalarField> tweights = tinterpScheme_().weights(vf);
     const edgeScalarField& weights = tweights();
 
-    tmp<faMatrix<Type>> tfam
+    auto tfam = tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            faceFlux.dimensions()*vf.dimensions()
-        )
+        vf,
+        faceFlux.dimensions()*vf.dimensions()
     );
     faMatrix<Type>& fam = tfam.ref();
 
diff --git a/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.H b/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.H
index be22ae9dbc1c689226ee67302eed430d168ee2b1..13302310c9e7583574c29f2825dddb5abadfc365 100644
--- a/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.H
+++ b/src/finiteArea/finiteArea/convectionSchemes/gaussFaConvectionScheme/gaussFaConvectionScheme.H
@@ -60,8 +60,9 @@ class gaussConvectionScheme
 :
     public fa::convectionScheme<Type>
 {
-    // Private data
+    // Private Data
 
+        //- Edge-interpolation scheme
         tmp<edgeInterpolationScheme<Type>> tinterpScheme_;
 
 
diff --git a/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.C b/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.C
index 1777a0f112ea98f6f86f995d08cb389dca63dac2..bf7955f54d644c1c41f8a0e52d9c6c0cb54f2abb 100644
--- a/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.C
+++ b/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.C
@@ -62,15 +62,14 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
     const dimensioned<Type> dt
 )
 {
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_();
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_();
-
-    dimensionedScalar rDeltaT2 =
+    const dimensionedScalar rDeltaT2 =
         4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
 
-    scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
-    scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
+    const scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
+    const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
 
     const IOobject d2dt2IOobject
     (
@@ -79,19 +78,16 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
 
     if (mesh().moving())
     {
-        scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
+        const scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
 
-        scalarField SS0 = mesh().S() + mesh().S0();
-        scalarField S0S00 = mesh().S0() + mesh().S00();
+        const scalarField SS0(mesh().S() + mesh().S0());
+        const scalarField S0S00(mesh().S0() + mesh().S00());
 
-        tmp<GeometricField<Type, faPatchField, areaMesh>> tdt2dt2
+        auto tdt2dt2 = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                d2dt2IOobject,
-                mesh(),
-                dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero)
-            )
+            d2dt2IOobject,
+            mesh(),
+            dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero)
         );
 
         tdt2dt2.ref().primitiveFieldRef() =
@@ -103,15 +99,12 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                d2dt2IOobject,
-                mesh(),
-                dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero),
-                calculatedFaPatchField<Type>::typeName
-            )
+            d2dt2IOobject,
+            mesh(),
+            dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero),
+            calculatedFaPatchField<Type>::typeName
         );
     }
 }
@@ -124,7 +117,7 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT2 =
+    const dimensionedScalar rDeltaT2 =
         4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
 
     const IOobject d2dt2IOobject
@@ -132,58 +125,50 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
         this->fieldIOobject("d2dt2("+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_();
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_();
 
-    scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
-    scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
+    const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
+    const scalar coefft0  = coefft + coefft00;
 
     if (mesh().moving())
     {
-        scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
+        const scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
 
-        scalarField SS0 = mesh().S() + mesh().S0();
-        scalarField S0S00 = mesh().S0() + mesh().S00();
+        const scalarField SS0(mesh().S() + mesh().S0());
+        const scalarField S0S00(mesh().S0() + mesh().S00());
 
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            d2dt2IOobject,
+            mesh(),
+            rDeltaT2.dimensions()*vf.dimensions(),
+            halfRdeltaT2*
+            (
+                coefft*SS0*vf.primitiveField()
+              - (coefft*SS0 + coefft00*S0S00)
+               *vf.oldTime().primitiveField()
+              + (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
+            )/mesh().S(),
+            rDeltaT2.value()*
             (
-                d2dt2IOobject,
-                mesh(),
-                rDeltaT2.dimensions()*vf.dimensions(),
-                halfRdeltaT2*
-                (
-                    coefft*SS0*vf.primitiveField()
-
-                  - (coefft*SS0 + coefft00*S0S00)
-                   *vf.oldTime().primitiveField()
-
-                  + (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
-                )/mesh().S(),
-                rDeltaT2.value()*
-                (
-                    coefft*vf.boundaryField()
-                  - coefft0*vf.oldTime().boundaryField()
-                  + coefft00*vf.oldTime().oldTime().boundaryField()
-                )
+                coefft*vf.boundaryField()
+              - coefft0*vf.oldTime().boundaryField()
+              + coefft00*vf.oldTime().oldTime().boundaryField()
             )
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            d2dt2IOobject,
+            rDeltaT2*
             (
-                d2dt2IOobject,
-                rDeltaT2*
-                (
-                    coefft*vf
-                  - coefft0*vf.oldTime()
-                  + coefft00*vf.oldTime().oldTime()
-                )
+                coefft*vf
+              - coefft0*vf.oldTime()
+              + coefft00*vf.oldTime().oldTime()
             )
         );
     }
@@ -198,7 +183,7 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT2 =
+    const dimensionedScalar rDeltaT2 =
         4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
 
     const IOobject d2dt2IOobject
@@ -206,65 +191,56 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
         this->fieldIOobject("d2dt2("+rho.name()+','+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_();
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_();
 
-    scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
-    scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
+    const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
+    const scalar coefft0  = coefft + coefft00;
 
     if (mesh().moving())
     {
-        scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
+        const scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
 
-        scalarField SS0rhoRho0 =
-            (mesh().S() + mesh().S0())
-           *rho.value();
+        const scalarField SS0rhoRho0((mesh().S() + mesh().S0())*rho.value());
 
-        scalarField S0S00rho0Rho00 =
+        const scalarField S0S00rho0Rho00
+        (
             (mesh().S0() + mesh().S00())
-           *rho.value();
+           *rho.value()
+        );
 
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            d2dt2IOobject,
+            mesh(),
+            rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
+            halfRdeltaT2*
             (
-                d2dt2IOobject,
-                mesh(),
-                rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
-                halfRdeltaT2*
-                (
-                    coefft*SS0rhoRho0*vf.primitiveField()
-
-                  - (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
-                   *vf.oldTime().primitiveField()
-
-                  + (coefft00*S0S00rho0Rho00)
-                   *vf.oldTime().oldTime().primitiveField()
-                )/mesh().S(),
-                rDeltaT2.value()*rho.value()*
-                (
-                    coefft*vf.boundaryField()
-                  - coefft0*vf.oldTime().boundaryField()
-                  + coefft00*vf.oldTime().oldTime().boundaryField()
-                )
+                coefft*SS0rhoRho0*vf.primitiveField()
+              - (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
+               *vf.oldTime().primitiveField()
+              + (coefft00*S0S00rho0Rho00)
+               *vf.oldTime().oldTime().primitiveField()
+            )/mesh().S(),
+            rDeltaT2.value()*rho.value()*
+            (
+                coefft*vf.boundaryField()
+              - coefft0*vf.oldTime().boundaryField()
+              + coefft00*vf.oldTime().oldTime().boundaryField()
             )
         );
     }
     else
     {
-
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            d2dt2IOobject,
+            rDeltaT2 * rho.value() *
             (
-                d2dt2IOobject,
-                rDeltaT2 * rho.value() *
-                (
-                    coefft*vf
-                  - coefft0*vf.oldTime()
-                  + coefft00*vf.oldTime().oldTime()
-                )
+                coefft*vf
+              - coefft0*vf.oldTime()
+              + coefft00*vf.oldTime().oldTime()
             )
         );
     }
@@ -279,7 +255,7 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT2 =
+    const dimensionedScalar rDeltaT2 =
         4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
 
     const IOobject d2dt2IOobject
@@ -287,24 +263,24 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
         this->fieldIOobject("d2dt2("+rho.name()+','+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_();
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_();
 
-    scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
-    scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
+    const scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
+    const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
 
     if (mesh().moving())
     {
-        scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
-        scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();
+        const scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
+        const scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();
 
-        scalarField SS0rhoRho0
+        const scalarField SS0rhoRho0
         (
             (mesh().S() + mesh().S0())
            *(rho.primitiveField() + rho.oldTime().primitiveField())
         );
 
-        scalarField S0S00rho0Rho00
+        const scalarField S0S00rho0Rho00
         (
             (mesh().S0() + mesh().S00())
            *(
@@ -313,69 +289,63 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
             )
         );
 
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            d2dt2IOobject,
+            mesh(),
+            rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
+            quarterRdeltaT2*
             (
-                d2dt2IOobject,
-                mesh(),
-                rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
-                quarterRdeltaT2*
-                (
-                    coefft*SS0rhoRho0*vf.primitiveField()
-
-                  - (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
-                   *vf.oldTime().primitiveField()
-
-                  + (coefft00*S0S00rho0Rho00)
-                   *vf.oldTime().oldTime().primitiveField()
-                )/mesh().S(),
-                halfRdeltaT2*
-                (
-                    coefft
-                   *(rho.boundaryField() + rho.oldTime().boundaryField())
-                   *vf.boundaryField()
-
-                  - (
-                        coefft
-                       *(
-                           rho.boundaryField()
-                         + rho.oldTime().boundaryField()
-                        )
-                      + coefft00
-                       *(
-                           rho.oldTime().boundaryField()
-                         + rho.oldTime().oldTime().boundaryField()
-                        )
-                    )*vf.oldTime().boundaryField()
+                coefft*SS0rhoRho0*vf.primitiveField()
 
+              - (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
+               *vf.oldTime().primitiveField()
+
+              + (coefft00*S0S00rho0Rho00)
+               *vf.oldTime().oldTime().primitiveField()
+            )/mesh().S(),
+            halfRdeltaT2*
+            (
+                coefft
+               *(rho.boundaryField() + rho.oldTime().boundaryField())
+               *vf.boundaryField()
+
+              - (
+                    coefft
+                   *(
+                        rho.boundaryField()
+                      + rho.oldTime().boundaryField()
+                    )
                   + coefft00
                    *(
-                       rho.oldTime().boundaryField()
-                     + rho.oldTime().oldTime().boundaryField()
-                    )*vf.oldTime().oldTime().boundaryField()
-                )
+                        rho.oldTime().boundaryField()
+                      + rho.oldTime().oldTime().boundaryField()
+                    )
+                )*vf.oldTime().boundaryField()
+
+              + coefft00
+               *(
+                    rho.oldTime().boundaryField()
+                  + rho.oldTime().oldTime().boundaryField()
+                )*vf.oldTime().oldTime().boundaryField()
             )
         );
     }
     else
     {
-        dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
+        const dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
 
-        areaScalarField rhoRho0(rho + rho.oldTime());
-        areaScalarField rho0Rho00(rho.oldTime() + rho.oldTime().oldTime());
+        const areaScalarField rhoRho0(rho + rho.oldTime());
+        const areaScalarField rho0Rho00(rho.oldTime() + rho.oldTime().oldTime());
 
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            d2dt2IOobject,
+            halfRdeltaT2*
             (
-                d2dt2IOobject,
-                halfRdeltaT2*
-                (
-                    coefft*rhoRho0*vf
-                  - (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
-                  + coefft00*rho0Rho00*vf.oldTime().oldTime()
-                )
+                coefft*rhoRho0*vf
+              - (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
+              + coefft00*rho0Rho00*vf.oldTime().oldTime()
             )
         );
     }
@@ -389,32 +359,28 @@ EulerFaD2dt2Scheme<Type>::famD2dt2
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    auto tfam = tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            vf.dimensions()*dimArea/dimTime/dimTime
-        )
+        vf,
+        vf.dimensions()*dimArea/dimTime/dimTime
     );
-
     faMatrix<Type>& fam = tfam.ref();
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_();
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_();
 
-    scalar coefft = (deltaT + deltaT0)/(2*deltaT);
-    scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
-    scalar coefft0 = coefft + coefft00;
+    const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
+    const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
+    const scalar coefft0 = coefft + coefft00;
 
-    scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
+    const scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
 
     if (mesh().moving())
     {
-        scalar halfRdeltaT2 = rDeltaT2/2.0;
+        const scalar halfRdeltaT2 = rDeltaT2/2.0;
 
-        scalarField SS0(mesh().S() + mesh().S0());
-        scalarField S0S00(mesh().S0() + mesh().S00());
+        const scalarField SS0(mesh().S() + mesh().S0());
+        const scalarField S0S00(mesh().S0() + mesh().S00());
 
         fam.diag() = (coefft*halfRdeltaT2)*SS0;
 
@@ -449,33 +415,27 @@ EulerFaD2dt2Scheme<Type>::famD2dt2
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    auto tfam = tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            rho.dimensions()*vf.dimensions()*dimArea
-            /dimTime/dimTime
-        )
+        vf,
+        rho.dimensions()*vf.dimensions()*dimArea/dimTime/dimTime
     );
-
     faMatrix<Type>& fam = tfam.ref();
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_();
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_();
 
-    scalar coefft = (deltaT + deltaT0)/(2*deltaT);
-    scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
+    const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
+    const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
 
-    scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
+    const scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
 
     if (mesh().moving())
     {
-        scalar halfRdeltaT2 = 0.5*rDeltaT2;
+        const scalar halfRdeltaT2 = 0.5*rDeltaT2;
 
-        scalarField SS0(mesh().S() + mesh().S0());
-
-        scalarField S0S00(mesh().S0() + mesh().S00());
+        const scalarField SS0(mesh().S() + mesh().S0());
+        const scalarField S0S00(mesh().S0() + mesh().S00());
 
         fam.diag() = rho.value()*(coefft*halfRdeltaT2)*SS0;
 
@@ -510,35 +470,32 @@ EulerFaD2dt2Scheme<Type>::famD2dt2
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    auto tfam = tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            rho.dimensions()*vf.dimensions()*dimArea/dimTime/dimTime
-        )
+        vf,
+        rho.dimensions()*vf.dimensions()*dimArea/dimTime/dimTime
     );
     faMatrix<Type>& fam = tfam.ref();
 
-    scalar deltaT =deltaT_();
-    scalar deltaT0 = deltaT0_();
+    const scalar deltaT =deltaT_();
+    const scalar deltaT0 = deltaT0_();
 
-    scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
-    scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
+    const scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
+    const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
 
-    scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
+    const scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
 
     if (mesh().moving())
     {
-        scalar quarterRdeltaT2 = 0.25*rDeltaT2;
+        const scalar quarterRdeltaT2 = 0.25*rDeltaT2;
 
-        scalarField SS0rhoRho0
+        const scalarField SS0rhoRho0
         (
             (mesh().S() + mesh().S0())
            *(rho.primitiveField() + rho.oldTime().primitiveField())
         );
 
-        scalarField S0S00rho0Rho00
+        const scalarField S0S00rho0Rho00
         (
             (mesh().S0() + mesh().S00())
            *(
@@ -560,14 +517,14 @@ EulerFaD2dt2Scheme<Type>::famD2dt2
     }
     else
     {
-        scalar halfRdeltaT2 = 0.5*rDeltaT2;
+        const scalar halfRdeltaT2 = 0.5*rDeltaT2;
 
-        scalarField rhoRho0
+        const scalarField rhoRho0
         (
             rho.primitiveField() + rho.oldTime().primitiveField()
         );
 
-        scalarField rho0Rho00
+        const scalarField rho0Rho00
         (
             rho.oldTime().primitiveField()
           + rho.oldTime().oldTime().primitiveField()
diff --git a/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.H b/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.H
index f25ca65a793220c7c685e88ce464e6ddc116c35c..ffea3feace039b70036c5c74fda3c476bbc01383 100644
--- a/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.H
+++ b/src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.H
@@ -98,7 +98,7 @@ public:
     // Member Functions
 
         //- Return mesh reference
-        const faMesh& mesh() const
+        const faMesh& mesh() const noexcept
         {
             return fa::faD2dt2Scheme<Type>::mesh();
         }
diff --git a/src/finiteArea/finiteArea/d2dt2Schemes/faD2dt2Scheme/faD2dt2Scheme.C b/src/finiteArea/finiteArea/d2dt2Schemes/faD2dt2Scheme/faD2dt2Scheme.C
index abe9ac38acddbc633353af9ee485666b51814af9..10ca27733a83a833396e4f0e551a9683172affb7 100644
--- a/src/finiteArea/finiteArea/d2dt2Schemes/faD2dt2Scheme/faD2dt2Scheme.C
+++ b/src/finiteArea/finiteArea/d2dt2Schemes/faD2dt2Scheme/faD2dt2Scheme.C
@@ -24,9 +24,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Abstract base class for finite area d2dt2 schemes.
-
 \*---------------------------------------------------------------------------*/
 
 #include "fa.H"
diff --git a/src/finiteArea/finiteArea/d2dt2Schemes/faD2dt2Scheme/faD2dt2Scheme.H b/src/finiteArea/finiteArea/d2dt2Schemes/faD2dt2Scheme/faD2dt2Scheme.H
index ff9ce8ec963be21399c50abfaecf318bceec3ec8..e6a72484db17cc664f82f95b937646928898e33a 100644
--- a/src/finiteArea/finiteArea/d2dt2Schemes/faD2dt2Scheme/faD2dt2Scheme.H
+++ b/src/finiteArea/finiteArea/d2dt2Schemes/faD2dt2Scheme/faD2dt2Scheme.H
@@ -49,6 +49,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 template<class Type>
 class faMatrix;
 
@@ -71,8 +72,9 @@ class faD2dt2Scheme
 {
 protected:
 
-    // Protected data
+    // Protected Data
 
+        //- Reference to the mesh
         const faMesh& mesh_;
 
 
@@ -138,10 +140,7 @@ public:
     // Member Functions
 
         //- Return mesh reference
-        const faMesh& mesh() const
-        {
-            return mesh_;
-        }
+        const faMesh& mesh() const noexcept { return mesh_; }
 
         virtual tmp<GeometricField<Type, faPatchField, areaMesh>> facD2dt2
         (
diff --git a/src/finiteArea/finiteArea/ddtSchemes/EulerFaDdtScheme/EulerFaDdtScheme.C b/src/finiteArea/finiteArea/ddtSchemes/EulerFaDdtScheme/EulerFaDdtScheme.C
index 037bf1fcc230026bea29e9a8adcaff8196946664..26233280df13eb125eb2f025bdaacd273372d5fa 100644
--- a/src/finiteArea/finiteArea/ddtSchemes/EulerFaDdtScheme/EulerFaDdtScheme.C
+++ b/src/finiteArea/finiteArea/ddtSchemes/EulerFaDdtScheme/EulerFaDdtScheme.C
@@ -48,7 +48,7 @@ EulerFaDdtScheme<Type>::facDdt
     const dimensioned<Type> dt
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
@@ -57,14 +57,11 @@ EulerFaDdtScheme<Type>::facDdt
 
     if (mesh().moving())
     {
-        tmp<GeometricField<Type, faPatchField, areaMesh>> tdtdt
+        auto tdtdt = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<Type>(dt.dimensions()/dimTime, Zero)
-            )
+            ddtIOobject,
+            mesh(),
+            dimensioned<Type>(dt.dimensions()/dimTime, Zero)
         );
 
         tdtdt.ref().primitiveFieldRef() =
@@ -91,21 +88,18 @@ EulerFaDdtScheme<Type>::facDdt0
     const dimensioned<Type> dt
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt("+dt.name()+')')
     );
 
-    tmp<GeometricField<Type, faPatchField, areaMesh>> tdtdt0
+    auto tdtdt0 = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
     (
-        new GeometricField<Type, faPatchField, areaMesh>
-        (
-            ddtIOobject,
-            mesh(),
-            -rDeltaT*dt
-        )
+        ddtIOobject,
+        mesh(),
+       -rDeltaT*dt
     );
 
     if (mesh().moving())
@@ -125,7 +119,7 @@ EulerFaDdtScheme<Type>::facDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
@@ -134,34 +128,28 @@ EulerFaDdtScheme<Type>::facDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*vf.dimensions(),
+            rDeltaT.value()*
+            (
+                vf()
+              - vf.oldTime()()*mesh().S0()/mesh().S()
+            ),
+            rDeltaT.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*vf.dimensions(),
-                rDeltaT.value()*
-                (
-                    vf()
-                  - vf.oldTime()()*mesh().S0()/mesh().S()
-                ),
-                rDeltaT.value()*
-                (
-                    vf.boundaryField() - vf.oldTime().boundaryField()
-                )
+                vf.boundaryField() - vf.oldTime().boundaryField()
             )
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                ddtIOobject,
-                rDeltaT*(vf - vf.oldTime())
-            )
+            ddtIOobject,
+            rDeltaT*(vf - vf.oldTime())
         );
     }
 }
@@ -174,7 +162,7 @@ EulerFaDdtScheme<Type>::facDdt0
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
@@ -183,27 +171,21 @@ EulerFaDdtScheme<Type>::facDdt0
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*vf.dimensions(),
-                (-rDeltaT.value())*vf.oldTime().internalField(),
-                (-rDeltaT.value())*vf.oldTime().boundaryField()
-            )
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*vf.dimensions(),
+            (-rDeltaT.value())*vf.oldTime().internalField(),
+            (-rDeltaT.value())*vf.oldTime().boundaryField()
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                ddtIOobject,
-                (-rDeltaT)*vf.oldTime()
-            )
+            ddtIOobject,
+            (-rDeltaT)*vf.oldTime()
         );
     }
 }
@@ -217,7 +199,7 @@ EulerFaDdtScheme<Type>::facDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
@@ -226,34 +208,28 @@ EulerFaDdtScheme<Type>::facDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            rDeltaT.value()*rho.value()*
+            (
+                vf()
+              - vf.oldTime()()*mesh().S0()/mesh().S()
+            ),
+            rDeltaT.value()*rho.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                rDeltaT.value()*rho.value()*
-                (
-                    vf()
-                  - vf.oldTime()()*mesh().S0()/mesh().S()
-                ),
-                rDeltaT.value()*rho.value()*
-                (
-                    vf.boundaryField() - vf.oldTime().boundaryField()
-                )
+                vf.boundaryField() - vf.oldTime().boundaryField()
             )
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                ddtIOobject,
-                rDeltaT*rho*(vf - vf.oldTime())
-            )
+            ddtIOobject,
+            rDeltaT*rho*(vf - vf.oldTime())
         );
     }
 }
@@ -266,7 +242,7 @@ EulerFaDdtScheme<Type>::facDdt0
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
@@ -275,29 +251,23 @@ EulerFaDdtScheme<Type>::facDdt0
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                (-rDeltaT.value())*rho.value()*
-                    vf.oldTime()()*mesh().S0()/mesh().S(),
-                (-rDeltaT.value())*rho.value()*
-                    vf.oldTime().boundaryField()
-            )
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            (-rDeltaT.value())*rho.value()*
+                vf.oldTime()()*mesh().S0()/mesh().S(),
+            (-rDeltaT.value())*rho.value()*
+                vf.oldTime().boundaryField()
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                ddtIOobject,
-                (-rDeltaT)*rho*vf.oldTime()
-            )
+            ddtIOobject,
+            (-rDeltaT)*rho*vf.oldTime()
         );
     }
 }
@@ -311,7 +281,7 @@ EulerFaDdtScheme<Type>::facDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
@@ -320,37 +290,31 @@ EulerFaDdtScheme<Type>::facDdt
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            rDeltaT.value()*
+            (
+                rho()*vf()
+              - rho.oldTime()()
+               *vf.oldTime()()*mesh().S0()/mesh().S()
+            ),
+            rDeltaT.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                rDeltaT.value()*
-                (
-                    rho()*vf()
-                  - rho.oldTime()()
-                   *vf.oldTime()()*mesh().S0()/mesh().S()
-                ),
-                rDeltaT.value()*
-                (
-                    rho.boundaryField()*vf.boundaryField()
-                  - rho.oldTime().boundaryField()
-                   *vf.oldTime().boundaryField()
-                )
+                rho.boundaryField()*vf.boundaryField()
+              - rho.oldTime().boundaryField()
+               *vf.oldTime().boundaryField()
             )
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                ddtIOobject,
-                rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
-            )
+            ddtIOobject,
+            rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
         );
     }
 }
@@ -364,7 +328,7 @@ EulerFaDdtScheme<Type>::facDdt0
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
@@ -373,35 +337,29 @@ EulerFaDdtScheme<Type>::facDdt0
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            rDeltaT.value()*
+            (
+              - rho.oldTime()()
+               *vf.oldTime()()*mesh().S0()/mesh().S()
+            ),
+            rDeltaT.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                rDeltaT.value()*
-                (
-                  - rho.oldTime()()
-                   *vf.oldTime()()*mesh().S0()/mesh().S()
-                ),
-                rDeltaT.value()*
-                (
-                  - rho.oldTime().boundaryField()
-                   *vf.oldTime().boundaryField()
-                )
+              - rho.oldTime().boundaryField()
+               *vf.oldTime().boundaryField()
             )
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                ddtIOobject,
-                (-rDeltaT)*rho.oldTime()*vf.oldTime()
-            )
+            ddtIOobject,
+            (-rDeltaT)*rho.oldTime()*vf.oldTime()
         );
     }
 }
@@ -414,18 +372,14 @@ EulerFaDdtScheme<Type>::famDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    auto tfam = tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        vf.dimensions()*dimArea/dimTime
     );
-
     faMatrix<Type>& fam = tfam.ref();
 
-    scalar rDeltaT = 1.0/mesh().time().deltaT().value();
+    const scalar rDeltaT = 1.0/mesh().time().deltaT().value();
 
     fam.diag() = rDeltaT*mesh().S();
 
@@ -450,17 +404,14 @@ EulerFaDdtScheme<Type>::famDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    auto tfam = tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            rho.dimensions()*vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        rho.dimensions()*vf.dimensions()*dimArea/dimTime
     );
     faMatrix<Type>& fam = tfam.ref();
 
-    scalar rDeltaT = 1.0/mesh().time().deltaT().value();
+    const scalar rDeltaT = 1.0/mesh().time().deltaT().value();
 
     fam.diag() = rDeltaT*rho.value()*mesh().S();
 
@@ -487,17 +438,14 @@ EulerFaDdtScheme<Type>::famDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    auto tfam = tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            rho.dimensions()*vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        rho.dimensions()*vf.dimensions()*dimArea/dimTime
     );
     faMatrix<Type>& fam = tfam.ref();
 
-    scalar rDeltaT = 1.0/mesh().time().deltaT().value();
+    const scalar rDeltaT = 1.0/mesh().time().deltaT().value();
 
     fam.diag() = rDeltaT*rho.primitiveField()*mesh().S();
 
diff --git a/src/finiteArea/finiteArea/ddtSchemes/backwardFaDdtScheme/backwardFaDdtScheme.C b/src/finiteArea/finiteArea/ddtSchemes/backwardFaDdtScheme/backwardFaDdtScheme.C
index 6dccbad1c35e0f2fcb8e48b3693be41ad24e3109..194c977a7ae4a895d52f3bb2112181ead0686d20 100644
--- a/src/finiteArea/finiteArea/ddtSchemes/backwardFaDdtScheme/backwardFaDdtScheme.C
+++ b/src/finiteArea/finiteArea/ddtSchemes/backwardFaDdtScheme/backwardFaDdtScheme.C
@@ -63,10 +63,8 @@ scalar backwardFaDdtScheme<Type>::deltaT0_(const GeoField& vf) const
     {
         return GREAT;
     }
-    else
-    {
-        return deltaT0_();
-    }
+
+    return deltaT0_();
 }
 
 
@@ -79,30 +77,27 @@ backwardFaDdtScheme<Type>::facDdt
     const dimensioned<Type> dt
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt("+dt.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_();
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_();
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
     if (mesh().moving())
     {
-        tmp<GeometricField<Type, faPatchField, areaMesh>> tdtdt
+        auto tdtdt = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<Type>(dt.dimensions()/dimTime, Zero)
-            )
+            ddtIOobject,
+            mesh(),
+            dimensioned<Type>(dt.dimensions()/dimTime, Zero)
         );
 
         tdtdt.ref().primitiveFieldRef() = rDeltaT.value()*dt.value()*
@@ -130,28 +125,25 @@ backwardFaDdtScheme<Type>::facDdt0
     const dimensioned<Type> dt
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt("+dt.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_();
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_();
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
-    tmp<GeometricField<Type, faPatchField, areaMesh>> tdtdt0
+    auto tdtdt0 = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
     (
-        new GeometricField<Type, faPatchField, areaMesh>
-        (
-            ddtIOobject,
-            mesh(),
-            -rDeltaT*(coefft0 - coefft00)*dt
-        )
+        ddtIOobject,
+        mesh(),
+       -rDeltaT*(coefft0 - coefft00)*dt
     );
 
     if (mesh().moving())
@@ -173,62 +165,56 @@ backwardFaDdtScheme<Type>::facDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt("+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*vf.dimensions(),
+            rDeltaT.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*vf.dimensions(),
-                rDeltaT.value()*
+                coefft*vf() -
                 (
-                    coefft*vf() -
-                    (
-                        coefft0*vf.oldTime()()*mesh().S0()
-                      - coefft00*vf.oldTime().oldTime()()
-                       *mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*
+                    coefft0*vf.oldTime()()*mesh().S0()
+                  - coefft00*vf.oldTime().oldTime()()
+                   *mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*
+            (
+                coefft*vf.boundaryField() -
                 (
-                    coefft*vf.boundaryField() -
-                    (
-                        coefft0*vf.oldTime().boundaryField()
-                      - coefft00*vf.oldTime().oldTime().boundaryField()
-                    )
+                    coefft0*vf.oldTime().boundaryField()
+                  - coefft00*vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            rDeltaT*
             (
-                ddtIOobject,
-                rDeltaT*
-                (
-                    coefft*vf
-                  - coefft0*vf.oldTime()
-                  + coefft00*vf.oldTime().oldTime()
-                )
+                coefft*vf
+              - coefft0*vf.oldTime()
+              + coefft00*vf.oldTime().oldTime()
             )
         );
     }
@@ -242,59 +228,53 @@ backwardFaDdtScheme<Type>::facDdt0
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt0("+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*vf.dimensions(),
+            rDeltaT.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*vf.dimensions(),
-                rDeltaT.value()*
-                (
-                  - (
-                        coefft0*vf.oldTime()()*mesh().S0()
-                      - coefft00*vf.oldTime().oldTime()()
-                       *mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*
-                (
-                  - (
-                        coefft0*vf.oldTime().boundaryField()
-                      - coefft00*vf.oldTime().oldTime().boundaryField()
-                    )
+              - (
+                    coefft0*vf.oldTime()()*mesh().S0()
+                  - coefft00*vf.oldTime().oldTime()()
+                   *mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*
+            (
+              - (
+                    coefft0*vf.oldTime().boundaryField()
+                  - coefft00*vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            rDeltaT*
             (
-                ddtIOobject,
-                rDeltaT*
-                (
-                  - coefft0*vf.oldTime()
-                  + coefft00*vf.oldTime().oldTime()
-                )
+                - coefft0*vf.oldTime()
+                + coefft00*vf.oldTime().oldTime()
             )
         );
     }
@@ -309,62 +289,56 @@ backwardFaDdtScheme<Type>::facDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            rDeltaT.value()*rho.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                rDeltaT.value()*rho.value()*
+                coefft*vf.internalField() -
                 (
-                    coefft*vf.internalField() -
-                    (
-                        coefft0*vf.oldTime()()*mesh().S0()
-                      - coefft00*vf.oldTime().oldTime()()
-                       *mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*rho.value()*
+                    coefft0*vf.oldTime()()*mesh().S0()
+                  - coefft00*vf.oldTime().oldTime()()
+                   *mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*rho.value()*
+            (
+                coefft*vf.boundaryField() -
                 (
-                    coefft*vf.boundaryField() -
-                    (
-                        coefft0*vf.oldTime().boundaryField()
-                      - coefft00*vf.oldTime().oldTime().boundaryField()
-                    )
+                    coefft0*vf.oldTime().boundaryField()
+                  - coefft00*vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            rDeltaT*rho*
             (
-                ddtIOobject,
-                rDeltaT*rho*
-                (
-                    coefft*vf
-                  - coefft0*vf.oldTime()
-                 + coefft00*vf.oldTime().oldTime()
-                )
+                coefft*vf
+              - coefft0*vf.oldTime()
+              + coefft00*vf.oldTime().oldTime()
             )
         );
     }
@@ -379,59 +353,53 @@ backwardFaDdtScheme<Type>::facDdt0
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            rDeltaT.value()*rho.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                rDeltaT.value()*rho.value()*
-                (
-                   -(
-                        coefft0*vf.oldTime()()*mesh().S0()
-                      - coefft00*vf.oldTime().oldTime()()
-                       *mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*rho.value()*
-                (
-                   -(
-                        coefft0*vf.oldTime().boundaryField()
-                      - coefft00*vf.oldTime().oldTime().boundaryField()
-                    )
+              - (
+                    coefft0*vf.oldTime()()*mesh().S0()
+                  - coefft00*vf.oldTime().oldTime()()
+                   *mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*rho.value()*
+            (
+              - (
+                    coefft0*vf.oldTime().boundaryField()
+                  - coefft00*vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            rDeltaT*rho*
             (
-                ddtIOobject,
-                rDeltaT*rho*
-                (
-                  - coefft0*vf.oldTime()
-                 + coefft00*vf.oldTime().oldTime()
-                )
+              - coefft0*vf.oldTime()
+              + coefft00*vf.oldTime().oldTime()
             )
         );
     }
@@ -446,65 +414,59 @@ backwardFaDdtScheme<Type>::facDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            rDeltaT.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                rDeltaT.value()*
+                coefft*rho.internalField()*vf.internalField() -
                 (
-                    coefft*rho.internalField()*vf.internalField() -
-                    (
-                        coefft0*rho.oldTime()()
-                       *vf.oldTime()()*mesh().S0()
-                      - coefft00*rho.oldTime().oldTime()()
-                       *vf.oldTime().oldTime()()*mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*
+                    coefft0*rho.oldTime()()
+                   *vf.oldTime()()*mesh().S0()
+                  - coefft00*rho.oldTime().oldTime()()
+                   *vf.oldTime().oldTime()()*mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*
+            (
+                coefft*vf.boundaryField() -
                 (
-                    coefft*vf.boundaryField() -
-                    (
-                        coefft0*rho.oldTime().boundaryField()
-                       *vf.oldTime().boundaryField()
-                      - coefft00*rho.oldTime().oldTime().boundaryField()
-                       *vf.oldTime().oldTime().boundaryField()
-                    )
+                    coefft0*rho.oldTime().boundaryField()
+                   *vf.oldTime().boundaryField()
+                  - coefft00*rho.oldTime().oldTime().boundaryField()
+                   *vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            rDeltaT*
             (
-                ddtIOobject,
-                rDeltaT*
-                (
-                    coefft*rho*vf
-                  - coefft0*rho.oldTime()*vf.oldTime()
-                  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
-                )
+                coefft*rho*vf
+              - coefft0*rho.oldTime()*vf.oldTime()
+              + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
             )
         );
     }
@@ -519,62 +481,56 @@ backwardFaDdtScheme<Type>::facDdt0
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
     if (mesh().moving())
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            rDeltaT.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                rDeltaT.value()*
-                (
-                  - (
-                        coefft0*rho.oldTime()()
-                       *vf.oldTime()()*mesh().S0()
-                      - coefft00*rho.oldTime().oldTime()()
-                       *vf.oldTime().oldTime()()*mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*
-                (
-                  - (
-                        coefft0*rho.oldTime().boundaryField()
-                       *vf.oldTime().boundaryField()
-                      - coefft00*rho.oldTime().oldTime().boundaryField()
-                       *vf.oldTime().oldTime().boundaryField()
-                    )
+              - (
+                    coefft0*rho.oldTime()()
+                   *vf.oldTime()()*mesh().S0()
+                  - coefft00*rho.oldTime().oldTime()()
+                   *vf.oldTime().oldTime()()*mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*
+            (
+              - (
+                    coefft0*rho.oldTime().boundaryField()
+                   *vf.oldTime().boundaryField()
+                  - coefft00*rho.oldTime().oldTime().boundaryField()
+                   *vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<GeometricField<Type, faPatchField, areaMesh>>
+        return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
-            new GeometricField<Type, faPatchField, areaMesh>
+            ddtIOobject,
+            rDeltaT*
             (
-                ddtIOobject,
-                rDeltaT*
-                (
-                  - coefft0*rho.oldTime()*vf.oldTime()
-                  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
-                )
+              - coefft0*rho.oldTime()*vf.oldTime()
+              + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
             )
         );
     }
@@ -588,25 +544,22 @@ backwardFaDdtScheme<Type>::famDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    auto tfam = tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        vf.dimensions()*dimArea/dimTime
     );
 
     faMatrix<Type>& fam = tfam.ref();
 
-    scalar rDeltaT = 1.0/deltaT_();
+    const scalar rDeltaT = 1.0/deltaT_();
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
     fam.diag() = (coefft*rDeltaT)*mesh().S();
 
@@ -640,24 +593,21 @@ backwardFaDdtScheme<Type>::famDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    auto tfam = tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            rho.dimensions()*vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        rho.dimensions()*vf.dimensions()*dimArea/dimTime
     );
     faMatrix<Type>& fam = tfam.ref();
 
-    scalar rDeltaT = 1.0/deltaT_();
+    const scalar rDeltaT = 1.0/deltaT_();
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
     fam.diag() = (coefft*rDeltaT*rho.value())*mesh().S();
 
@@ -691,24 +641,21 @@ backwardFaDdtScheme<Type>::famDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    auto tfam = tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            rho.dimensions()*vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        rho.dimensions()*vf.dimensions()*dimArea/dimTime
     );
     faMatrix<Type>& fam = tfam.ref();
 
-    scalar rDeltaT = 1.0/deltaT_();
+    const scalar rDeltaT = 1.0/deltaT_();
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
     fam.diag() = (coefft*rDeltaT)*rho.primitiveField()*mesh().S();
 
diff --git a/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.C b/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.C
index 2fb16a27838c5981f85a2408295a881f757bc530..84058600b46fad8632f6a961a85ebd607f867a2b 100644
--- a/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.C
+++ b/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.C
@@ -62,30 +62,27 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
 {
     // No change compared to backward
 
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt("+dt.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_();
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_();
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
     if (mesh().moving())
     {
-        tmp<areaScalarField> tdtdt
+        auto tdtdt = tmp<areaScalarField>::New
         (
-            new areaScalarField
-            (
-                ddtIOobject,
-                mesh(),
-                dimensionedScalar(dt.dimensions()/dimTime, Zero)
-            )
+            ddtIOobject,
+            mesh(),
+            dimensionedScalar(dt.dimensions()/dimTime, Zero)
         );
 
         tdtdt.ref().primitiveFieldRef() = rDeltaT.value()*dt.value()*
@@ -113,28 +110,25 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
 {
     // No change compared to backward
 
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt("+dt.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_();
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_();
 
-    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
-    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
-    scalar coefft0  = coefft + coefft00;
+    const scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    const scalar coefft0  = coefft + coefft00;
 
-    tmp<areaScalarField> tdtdt0
+    auto tdtdt0 = tmp<areaScalarField>::New
     (
-        new areaScalarField
-        (
-            ddtIOobject,
-            mesh(),
-            -rDeltaT*(coefft0 - coefft00)*dt
-        )
+        ddtIOobject,
+        mesh(),
+       -rDeltaT*(coefft0 - coefft00)*dt
     );
 
     if (mesh().moving())
@@ -154,20 +148,20 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
     const areaScalarField& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt("+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
     // Calculate unboundedness indicator
     // Note: all times moved by one because access to internal field
     // copies current field into the old-time level.
-    areaScalarField phict
+    const areaScalarField phict
     (
         mag
         (
@@ -180,62 +174,59 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
                 vf.oldTime()
               - vf.oldTime().oldTime()
             )
-          + dimensionedScalar("small", vf.dimensions(), SMALL)
+          + dimensionedScalar(vf.dimensions(), SMALL)
         )
     );
 
-    areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
+    const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
 
-    areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
-    areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
-    areaScalarField coefft0(coefft + coefft00);
+    const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
+    const areaScalarField coefft00
+    (
+        limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
+    );
+    const areaScalarField coefft0(coefft + coefft00);
 
     if (mesh().moving())
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*vf.dimensions(),
+            rDeltaT.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*vf.dimensions(),
-                rDeltaT.value()*
+                coefft*vf.primitiveField() -
                 (
-                    coefft*vf.primitiveField() -
-                    (
-                        coefft0.primitiveField()
-                        *vf.oldTime().primitiveField()*mesh().S0()
-                      - coefft00.primitiveField()
-                        *vf.oldTime().oldTime().primitiveField()
-                       *mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*
+                    coefft0.primitiveField()
+                    *vf.oldTime().primitiveField()*mesh().S0()
+                  - coefft00.primitiveField()
+                    *vf.oldTime().oldTime().primitiveField()
+                    *mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*
+            (
+                coefft.boundaryField()*vf.boundaryField() -
                 (
-                    coefft.boundaryField()*vf.boundaryField() -
-                    (
-                        coefft0.boundaryField()*
-                        vf.oldTime().boundaryField()
-                      - coefft00.boundaryField()*
-                        vf.oldTime().oldTime().boundaryField()
-                    )
+                    coefft0.boundaryField()*
+                    vf.oldTime().boundaryField()
+                  - coefft00.boundaryField()*
+                    vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            rDeltaT*
             (
-                ddtIOobject,
-                rDeltaT*
-                (
-                    coefft*vf
-                  - coefft0*vf.oldTime()
-                  + coefft00*vf.oldTime().oldTime()
-                )
+                coefft*vf
+              - coefft0*vf.oldTime()
+              + coefft00*vf.oldTime().oldTime()
             )
         );
     }
@@ -247,20 +238,20 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
     const areaScalarField& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt0("+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
     // Calculate unboundedness indicator
     // Note: all times moved by one because access to internal field
     // copies current field into the old-time level.
-    areaScalarField phict
+    const areaScalarField phict
     (
         mag
         (
@@ -273,59 +264,56 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
                 vf.oldTime()
               - vf.oldTime().oldTime()
             )
-          + dimensionedScalar("small", vf.dimensions(), SMALL)
+          + dimensionedScalar(vf.dimensions(), SMALL)
         )
     );
 
-    areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
+    const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
 
-    areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
-    areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
-    areaScalarField coefft0(coefft + coefft00);
+    const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
+    const areaScalarField coefft00
+    (
+        limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
+    );
+    const areaScalarField coefft0(coefft + coefft00);
 
     if (mesh().moving())
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*vf.dimensions(),
+            rDeltaT.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*vf.dimensions(),
-                rDeltaT.value()*
-                (
-                  - (
-                        coefft0.primitiveField()*
-                        vf.oldTime().primitiveField()*mesh().S0()
-                      - coefft00.primitiveField()*
-                        vf.oldTime().oldTime().primitiveField()
-                       *mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*
-                (
-                  - (
-                        coefft0.boundaryField()*
-                        vf.oldTime().boundaryField()
-                      - coefft00.boundaryField()*
-                        vf.oldTime().oldTime().boundaryField()
-                    )
+              - (
+                    coefft0.primitiveField()*
+                    vf.oldTime().primitiveField()*mesh().S0()
+                  - coefft00.primitiveField()*
+                    vf.oldTime().oldTime().primitiveField()
+                    *mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*
+            (
+              - (
+                    coefft0.boundaryField()*
+                    vf.oldTime().boundaryField()
+                  - coefft00.boundaryField()*
+                    vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            rDeltaT*
             (
-                ddtIOobject,
-                rDeltaT*
-                (
-                  - coefft0*vf.oldTime()
-                  + coefft00*vf.oldTime().oldTime()
-                )
+              - coefft0*vf.oldTime()
+              + coefft00*vf.oldTime().oldTime()
             )
         );
     }
@@ -338,20 +326,20 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
     const areaScalarField& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
     // Calculate unboundedness indicator
     // Note: all times moved by one because access to internal field
     // copies current field into the old-time level.
-    areaScalarField phict
+    const areaScalarField phict
     (
         mag
         (
@@ -368,83 +356,81 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
         )
     );
 
-    areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
+    const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
 
-    areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
-    areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
-    areaScalarField coefft0(coefft + coefft00);
+    const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
+    const areaScalarField coefft00
+    (
+        limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
+    );
+    const areaScalarField coefft0(coefft + coefft00);
 
     if (mesh().moving())
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            rDeltaT.value()*rho.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                rDeltaT.value()*rho.value()*
+                coefft*vf.primitiveField() -
                 (
-                    coefft*vf.primitiveField() -
-                    (
-                        coefft0.primitiveField()*
-                        vf.oldTime().primitiveField()*mesh().S0()
-                      - coefft00.primitiveField()*
-                        vf.oldTime().oldTime().primitiveField()
-                       *mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*rho.value()*
+                    coefft0.primitiveField()*
+                    vf.oldTime().primitiveField()*mesh().S0()
+                  - coefft00.primitiveField()*
+                    vf.oldTime().oldTime().primitiveField()
+                    *mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*rho.value()*
+            (
+                coefft.boundaryField()*vf.boundaryField() -
                 (
-                    coefft.boundaryField()*vf.boundaryField() -
-                    (
-                        coefft0.boundaryField()*
-                        vf.oldTime().boundaryField()
-                      - coefft00.boundaryField()*
-                        vf.oldTime().oldTime().boundaryField()
-                    )
+                    coefft0.boundaryField()*
+                    vf.oldTime().boundaryField()
+                  - coefft00.boundaryField()*
+                    vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            rDeltaT*rho*
             (
-                ddtIOobject,
-                rDeltaT*rho*
-                (
-                    coefft*vf
-                  - coefft0*vf.oldTime()
-                 + coefft00*vf.oldTime().oldTime()
-                )
+                coefft*vf
+              - coefft0*vf.oldTime()
+              + coefft00*vf.oldTime().oldTime()
             )
         );
     }
 }
 
+
 tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
 (
     const dimensionedScalar& rho,
     const areaScalarField& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
     // Calculate unboundedness indicator
     // Note: all times moved by one because access to internal field
     // copies current field into the old-time level.
-    areaScalarField phict
+    const areaScalarField phict
     (
         mag
         (
@@ -457,59 +443,56 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
                 vf.oldTime()
               - vf.oldTime().oldTime()
             )
-          + dimensionedScalar("small", vf.dimensions(), SMALL)
+          + dimensionedScalar(vf.dimensions(), SMALL)
         )
     );
 
-    areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
+    const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
 
-    areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
-    areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
-    areaScalarField coefft0(coefft + coefft00);
+    const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
+    const areaScalarField coefft00
+    (
+        limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
+    );
+    const areaScalarField coefft0(coefft + coefft00);
 
     if (mesh().moving())
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            rDeltaT.value()*rho.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                rDeltaT.value()*rho.value()*
-                (
-                   -(
-                        coefft0.primitiveField()*
-                        vf.oldTime().primitiveField()*mesh().S0()
-                      - coefft00.primitiveField()*
-                        vf.oldTime().oldTime().primitiveField()
-                       *mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*rho.value()*
-                (
-                   -(
-                        coefft0.boundaryField()*
-                        vf.oldTime().boundaryField()
-                      - coefft00.boundaryField()*
-                        vf.oldTime().oldTime().boundaryField()
-                    )
+              -(
+                    coefft0.primitiveField()*
+                    vf.oldTime().primitiveField()*mesh().S0()
+                  - coefft00.primitiveField()*
+                    vf.oldTime().oldTime().primitiveField()
+                    *mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*rho.value()*
+            (
+              -(
+                    coefft0.boundaryField()*
+                    vf.oldTime().boundaryField()
+                  - coefft00.boundaryField()*
+                    vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            rDeltaT*rho*
             (
-                ddtIOobject,
-                rDeltaT*rho*
-                (
-                  - coefft0*vf.oldTime()
-                 + coefft00*vf.oldTime().oldTime()
-                )
+              - coefft0*vf.oldTime()
+              + coefft00*vf.oldTime().oldTime()
             )
         );
     }
@@ -522,20 +505,20 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
     const areaScalarField& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
     // Calculate unboundedness indicator
     // Note: all times moved by one because access to internal field
     // copies current field into the old-time level.
-    areaScalarField phict
+    const areaScalarField phict
     (
         mag
         (
@@ -548,65 +531,62 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
                 rho.oldTime()*vf.oldTime()
               - rho.oldTime().oldTime()*vf.oldTime().oldTime()
             )
-          + dimensionedScalar("small", rho.dimensions()*vf.dimensions(), SMALL)
+          + dimensionedScalar(rho.dimensions()*vf.dimensions(), SMALL)
         )
     );
 
-    areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
+    const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
 
-    areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
-    areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
-    areaScalarField coefft0(coefft + coefft00);
+    const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
+    const areaScalarField coefft00
+    (
+        limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
+    );
+    const areaScalarField coefft0(coefft + coefft00);
 
     if (mesh().moving())
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            rDeltaT.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                rDeltaT.value()*
+                coefft*rho.primitiveField()*vf.primitiveField() -
                 (
-                    coefft*rho.primitiveField()*vf.primitiveField() -
-                    (
-                        coefft0.primitiveField()*
-                        rho.oldTime().primitiveField()*
-                        vf.oldTime().primitiveField()*mesh().S0()
-                      - coefft00.primitiveField()*
-                        rho.oldTime().oldTime().primitiveField()
-                       *vf.oldTime().oldTime().primitiveField()*mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*
+                    coefft0.primitiveField()*
+                    rho.oldTime().primitiveField()*
+                    vf.oldTime().primitiveField()*mesh().S0()
+                  - coefft00.primitiveField()*
+                    rho.oldTime().oldTime().primitiveField()
+                    *vf.oldTime().oldTime().primitiveField()*mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*
+            (
+                coefft.boundaryField()*vf.boundaryField() -
                 (
-                    coefft.boundaryField()*vf.boundaryField() -
-                    (
-                        coefft0.boundaryField()*
-                        rho.oldTime().boundaryField()*
-                        vf.oldTime().boundaryField()
-                      - coefft00.boundaryField()*
-                        rho.oldTime().oldTime().boundaryField()*
-                        vf.oldTime().oldTime().boundaryField()
-                    )
+                    coefft0.boundaryField()*
+                    rho.oldTime().boundaryField()*
+                    vf.oldTime().boundaryField()
+                  - coefft00.boundaryField()*
+                    rho.oldTime().oldTime().boundaryField()*
+                    vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            rDeltaT*
             (
-                ddtIOobject,
-                rDeltaT*
-                (
-                    coefft*rho*vf
-                  - coefft0*rho.oldTime()*vf.oldTime()
-                  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
-                )
+                coefft*rho*vf
+              - coefft0*rho.oldTime()*vf.oldTime()
+              + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
             )
         );
     }
@@ -619,20 +599,20 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
     const areaScalarField& vf
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+    const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
 
     const IOobject ddtIOobject
     (
         this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')')
     );
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
     // Calculate unboundedness indicator
     // Note: all times moved by one because access to internal field
     // copies current field into the old-time level.
-    areaScalarField phict
+    const areaScalarField phict
     (
         mag
         (
@@ -645,62 +625,59 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
                 rho.oldTime()*vf.oldTime()
               - rho.oldTime().oldTime()*vf.oldTime().oldTime()
             )
-          + dimensionedScalar("small", rho.dimensions()*vf.dimensions(), SMALL)
+          + dimensionedScalar(rho.dimensions()*vf.dimensions(), SMALL)
         )
     );
 
-    areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
+    const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
 
-    areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
-    areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
-    areaScalarField coefft0(coefft + coefft00);
+    const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
+    const areaScalarField coefft00
+    (
+        limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
+    );
+    const areaScalarField coefft0(coefft + coefft00);
 
     if (mesh().moving())
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            mesh(),
+            rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+            rDeltaT.value()*
             (
-                ddtIOobject,
-                mesh(),
-                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
-                rDeltaT.value()*
-                (
-                  - (
-                        coefft0.primitiveField()*
-                        rho.oldTime().primitiveField()*
-                        vf.oldTime().primitiveField()*mesh().S0()
-                      - coefft00.primitiveField()*
-                        rho.oldTime().oldTime().primitiveField()*
-                        vf.oldTime().oldTime().primitiveField()*mesh().S00()
-                    )/mesh().S()
-                ),
-                rDeltaT.value()*
-                (
-                  - (
-                        coefft0.boundaryField()*
-                        rho.oldTime().boundaryField()*
-                        vf.oldTime().boundaryField()
-                      - coefft00.boundaryField()*
-                        rho.oldTime().oldTime().boundaryField()*
-                        vf.oldTime().oldTime().boundaryField()
-                    )
+              - (
+                    coefft0.primitiveField()*
+                    rho.oldTime().primitiveField()*
+                    vf.oldTime().primitiveField()*mesh().S0()
+                  - coefft00.primitiveField()*
+                    rho.oldTime().oldTime().primitiveField()*
+                    vf.oldTime().oldTime().primitiveField()*mesh().S00()
+                )/mesh().S()
+            ),
+            rDeltaT.value()*
+            (
+              - (
+                    coefft0.boundaryField()*
+                    rho.oldTime().boundaryField()*
+                    vf.oldTime().boundaryField()
+                  - coefft00.boundaryField()*
+                    rho.oldTime().oldTime().boundaryField()*
+                    vf.oldTime().oldTime().boundaryField()
                 )
             )
         );
     }
     else
     {
-        return tmp<areaScalarField>
+        return tmp<areaScalarField>::New
         (
-            new areaScalarField
+            ddtIOobject,
+            rDeltaT*
             (
-                ddtIOobject,
-                rDeltaT*
-                (
-                  - coefft0*rho.oldTime()*vf.oldTime()
-                  + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
-                )
+              - coefft0*rho.oldTime()*vf.oldTime()
+              + coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
             )
         );
     }
@@ -712,26 +689,22 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
     const areaScalarField& vf
 )
 {
-    tmp<faScalarMatrix> tfam
+    auto tfam = tmp<faScalarMatrix>::New
     (
-        new faScalarMatrix
-        (
-            vf,
-            vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        vf.dimensions()*dimArea/dimTime
     );
-
     faScalarMatrix& fam = tfam.ref();
 
-    scalar rDeltaT = 1.0/deltaT_();
+    const scalar rDeltaT = 1.0/deltaT_();
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
     // Calculate unboundedness indicator
     // Note: all times moved by one because access to internal field
     // copies current field into the old-time level.
-    scalarField phict
+    const scalarField phict
     (
         mag
         (
@@ -748,11 +721,14 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
         )
     );
 
-    scalarField limiter(pos(phict) - pos(phict - 1.0));
+    const scalarField limiter(pos(phict) - pos(phict - 1.0));
 
-    scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
-    scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)));
-    scalarField coefft0(coefft + coefft00);
+    const scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
+    const scalarField coefft00
+    (
+        limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))
+    );
+    const scalarField coefft0(coefft + coefft00);
 
     fam.diag() = (coefft*rDeltaT)*mesh().S();
 
@@ -784,25 +760,22 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
     const areaScalarField& vf
 )
 {
-    tmp<faScalarMatrix> tfam
+    auto tfam = tmp<faScalarMatrix>::New
     (
-        new faScalarMatrix
-        (
-            vf,
-            rho.dimensions()*vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        rho.dimensions()*vf.dimensions()*dimArea/dimTime
     );
     faScalarMatrix& fam = tfam.ref();
 
-    scalar rDeltaT = 1.0/deltaT_();
+    const scalar rDeltaT = 1.0/deltaT_();
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
     // Calculate unboundedness indicator
     // Note: all times moved by one because access to internal field
     // copies current field into the old-time level.
-    scalarField phict
+    const scalarField phict
     (
         mag
         (
@@ -819,11 +792,14 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
         )
     );
 
-    scalarField limiter(pos(phict) - pos(phict - 1.0));
+    const scalarField limiter(pos(phict) - pos(phict - 1.0));
 
-    scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
-    scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)));
-    scalarField coefft0(coefft + coefft00);
+    const scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
+    const scalarField coefft00
+    (
+        limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))
+    );
+    const scalarField coefft0(coefft + coefft00);
 
     fam.diag() = (coefft*rDeltaT*rho.value())*mesh().S();
 
@@ -855,25 +831,22 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
     const areaScalarField& vf
 )
 {
-    tmp<faScalarMatrix> tfam
+    auto tfam = tmp<faScalarMatrix>::New
     (
-        new faScalarMatrix
-        (
-            vf,
-            rho.dimensions()*vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        rho.dimensions()*vf.dimensions()*dimArea/dimTime
     );
     faScalarMatrix& fam = tfam.ref();
 
-    scalar rDeltaT = 1.0/deltaT_();
+    const scalar rDeltaT = 1.0/deltaT_();
 
-    scalar deltaT = deltaT_();
-    scalar deltaT0 = deltaT0_(vf);
+    const scalar deltaT = deltaT_();
+    const scalar deltaT0 = deltaT0_(vf);
 
     // Calculate unboundedness indicator
     // Note: all times moved by one because access to internal field
     // copies current field into the old-time level.
-    scalarField phict
+    const scalarField phict
     (
         mag
         (
@@ -894,11 +867,14 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
         )
     );
 
-    scalarField limiter(pos(phict) - pos(phict - 1.0));
+    const scalarField limiter(pos(phict) - pos(phict - 1.0));
 
-    scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
-    scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)));
-    scalarField coefft0(coefft + coefft00);
+    const scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
+    const scalarField coefft00
+    (
+        limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))
+    );
+    const scalarField coefft0(coefft + coefft00);
 
     fam.diag() = (coefft*rDeltaT)*rho.primitiveField()*mesh().S();
 
diff --git a/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.H b/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.H
index 02163a746be1ac4a0befd31e91fc23df5e88df39..d284c03f1327a1815e24a7dabf78d558ce53d113 100644
--- a/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.H
+++ b/src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.H
@@ -76,10 +76,8 @@ class boundedBackwardFaDdtScheme
             {
                 return GREAT;
             }
-            else
-            {
-                return deltaT0_();
-            }
+
+            return deltaT0_();
         }
 
 
@@ -114,7 +112,7 @@ public:
     // Member Functions
 
         //- Return mesh reference
-        const faMesh& mesh() const
+        const faMesh& mesh() const noexcept
         {
             return fa::faDdtScheme<scalar>::mesh();
         }
diff --git a/src/finiteArea/finiteArea/ddtSchemes/faDdtScheme/faDdtScheme.H b/src/finiteArea/finiteArea/ddtSchemes/faDdtScheme/faDdtScheme.H
index 867dee9963133480f699d3a43f28a971bebfc1ff..d7137be1eeb1f95ab38e27a181cc451d99f35459 100644
--- a/src/finiteArea/finiteArea/ddtSchemes/faDdtScheme/faDdtScheme.H
+++ b/src/finiteArea/finiteArea/ddtSchemes/faDdtScheme/faDdtScheme.H
@@ -49,6 +49,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 template<class Type>
 class faMatrix;
 
@@ -71,8 +72,9 @@ class faDdtScheme
 {
 protected:
 
-    // Protected data
+    // Protected Data
 
+        //- Reference to the mesh
         const faMesh& mesh_;
 
 
@@ -138,10 +140,7 @@ public:
     // Member Functions
 
         //- Return mesh reference
-        const faMesh& mesh() const
-        {
-            return mesh_;
-        }
+        const faMesh& mesh() const noexcept { return mesh_; }
 
         virtual tmp<GeometricField<Type, faPatchField, areaMesh>> facDdt
         (
diff --git a/src/finiteArea/finiteArea/ddtSchemes/steadyStateFaDdtScheme/steadyStateFaDdtScheme.C b/src/finiteArea/finiteArea/ddtSchemes/steadyStateFaDdtScheme/steadyStateFaDdtScheme.C
index f4d403278f72af6e192dbf2a2497cc51bfa0e9b6..b007ed3e515b4c7d00d250f63312ec2c1bac7211 100644
--- a/src/finiteArea/finiteArea/ddtSchemes/steadyStateFaDdtScheme/steadyStateFaDdtScheme.C
+++ b/src/finiteArea/finiteArea/ddtSchemes/steadyStateFaDdtScheme/steadyStateFaDdtScheme.C
@@ -180,16 +180,11 @@ steadyStateFaDdtScheme<Type>::famDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    return tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        vf.dimensions()*dimArea/dimTime
     );
-
-    return tfam;
 }
 
 
@@ -201,16 +196,11 @@ steadyStateFaDdtScheme<Type>::famDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    return tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            rho.dimensions()*vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        rho.dimensions()*vf.dimensions()*dimArea/dimTime
     );
-
-    return tfam;
 }
 
 
@@ -222,16 +212,11 @@ steadyStateFaDdtScheme<Type>::famDdt
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-    tmp<faMatrix<Type>> tfam
+    return tmp<faMatrix<Type>>::New
     (
-        new faMatrix<Type>
-        (
-            vf,
-            rho.dimensions()*vf.dimensions()*dimArea/dimTime
-        )
+        vf,
+        rho.dimensions()*vf.dimensions()*dimArea/dimTime
     );
-
-    return tfam;
 }
 
 
diff --git a/src/finiteArea/finiteArea/divSchemes/faDivScheme/faDivScheme.H b/src/finiteArea/finiteArea/divSchemes/faDivScheme/faDivScheme.H
index 56f4abc5342b453db51a2bbf03fa2b6e88447add..905a21da3d84bc221419d42f657cec5c63d241d2 100644
--- a/src/finiteArea/finiteArea/divSchemes/faDivScheme/faDivScheme.H
+++ b/src/finiteArea/finiteArea/divSchemes/faDivScheme/faDivScheme.H
@@ -50,6 +50,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 template<class Type>
 class faMatrix;
 
@@ -69,16 +70,18 @@ class divScheme
 :
     public refCount
 {
-
 protected:
 
     // Protected data
 
+        //- Reference to the mesh
         const faMesh& mesh_;
+
+        //- Edge-interpolation scheme
         tmp<edgeInterpolationScheme<Type>> tinterpScheme_;
 
 
-    // Private Member Functions
+    // Protected Member Functions
 
         //- No copy construct
         divScheme(const divScheme&) = delete;
@@ -135,10 +138,7 @@ public:
     // Member Functions
 
         //- Return mesh reference
-        const faMesh& mesh() const
-        {
-            return mesh_;
-        }
+        const faMesh& mesh() const noexcept { return mesh_; }
 
         virtual tmp
         <
diff --git a/src/finiteArea/finiteArea/fac/facAverage.C b/src/finiteArea/finiteArea/fac/facAverage.C
index 57c4e78bbb2a4decaa0b0cf71c6cf88789135b1c..1079ae2c1c3e30b2c9dabb9ae16d2291d1d1aa46 100644
--- a/src/finiteArea/finiteArea/fac/facAverage.C
+++ b/src/finiteArea/finiteArea/fac/facAverage.C
@@ -50,23 +50,19 @@ average
 {
     const faMesh& mesh = ssf.mesh();
 
-    tmp<GeometricField<Type, faPatchField, areaMesh>> taverage
+    auto taverage = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
     (
-        new GeometricField<Type, faPatchField, areaMesh>
+        IOobject
         (
-            IOobject
-            (
-                "average("+ssf.name()+')',
-                ssf.instance(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            ssf.dimensions()
-        )
+            "average("+ssf.name()+')',
+            ssf.instance(),
+            mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        ssf.dimensions()
     );
-
     GeometricField<Type, faPatchField, areaMesh>& av = taverage.ref();
 
     av.ref() =
diff --git a/src/finiteArea/finiteArea/fac/facAverage.H b/src/finiteArea/finiteArea/fac/facAverage.H
index 5ff0a1993af14114a13ccc8d68d649e674c53d12..1aac0cdb2728235835baa8b8c247ae3dfe5e4f8f 100644
--- a/src/finiteArea/finiteArea/fac/facAverage.H
+++ b/src/finiteArea/finiteArea/fac/facAverage.H
@@ -34,7 +34,6 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-
 #ifndef facAverage_H
 #define facAverage_H
 
diff --git a/src/finiteArea/finiteArea/fac/facDdt.H b/src/finiteArea/finiteArea/fac/facDdt.H
index 64c6d35d7605cf9bea99fbf7bc2b2e54f9e0e5af..2b9c872477b36a1673f6d4fcaaf655b1dcc0062f 100644
--- a/src/finiteArea/finiteArea/fac/facDdt.H
+++ b/src/finiteArea/finiteArea/fac/facDdt.H
@@ -37,7 +37,6 @@ Author
 
 \*---------------------------------------------------------------------------*/
 
-
 #ifndef facDdt_H
 #define facDdt_H
 
diff --git a/src/finiteArea/finiteArea/fac/facDiv.H b/src/finiteArea/finiteArea/fac/facDiv.H
index d98e2d4ca7b0e28ac0caf7d5f16fdcbcdbb9a237..b3f4ecadba7581ae4eaa51754c4849e9b819c635 100644
--- a/src/finiteArea/finiteArea/fac/facDiv.H
+++ b/src/finiteArea/finiteArea/fac/facDiv.H
@@ -34,7 +34,6 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-
 #ifndef facDiv_H
 #define facDiv_H
 
diff --git a/src/finiteArea/finiteArea/fac/facEdgeIntegrate.C b/src/finiteArea/finiteArea/fac/facEdgeIntegrate.C
index 3cae46f5eecc369a6985269b4cb00e440bbdfa0a..66347759be1f2ebacdb1001078874994348a5169 100644
--- a/src/finiteArea/finiteArea/fac/facEdgeIntegrate.C
+++ b/src/finiteArea/finiteArea/fac/facEdgeIntegrate.C
@@ -50,20 +50,17 @@ edgeIntegrate
 {
     const faMesh& mesh = ssf.mesh();
 
-    tmp<GeometricField<Type, faPatchField, areaMesh>> tvf
+    auto tvf = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
     (
-        new GeometricField<Type, faPatchField, areaMesh>
+        IOobject
         (
-            IOobject
-            (
-                "edgeIntegrate("+ssf.name()+')',
-                ssf.instance(),
-                ssf.db()
-            ),
-            mesh,
-            dimensioned<Type>(ssf.dimensions()/dimArea, Zero),
-            zeroGradientFaPatchField<Type>::typeName
-        )
+            "edgeIntegrate("+ssf.name()+')',
+            ssf.instance(),
+            ssf.db()
+        ),
+        mesh,
+        dimensioned<Type>(ssf.dimensions()/dimArea, Zero),
+        zeroGradientFaPatchField<Type>::typeName
     );
     GeometricField<Type, faPatchField, areaMesh>& vf = tvf.ref();
 
@@ -122,20 +119,17 @@ edgeSum
 {
     const faMesh& mesh = ssf.mesh();
 
-    tmp<GeometricField<Type, faPatchField, areaMesh>> tvf
+    auto tvf = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
     (
-        new GeometricField<Type, faPatchField, areaMesh>
+        IOobject
         (
-            IOobject
-            (
-                "edgeSum("+ssf.name()+')',
-                ssf.instance(),
-                ssf.db()
-            ),
-            mesh,
-            dimensioned<Type>(ssf.dimensions(), Zero),
-            zeroGradientFaPatchField<Type>::typeName
-        )
+            "edgeSum("+ssf.name()+')',
+            ssf.instance(),
+            ssf.db()
+        ),
+        mesh,
+        dimensioned<Type>(ssf.dimensions(), Zero),
+        zeroGradientFaPatchField<Type>::typeName
     );
     GeometricField<Type, faPatchField, areaMesh>& vf = tvf.ref();
 
diff --git a/src/finiteArea/finiteArea/fac/facEdgeIntegrate.H b/src/finiteArea/finiteArea/fac/facEdgeIntegrate.H
index d7acc4231c2fcaa1a4c234416605770ce7965f94..b9a1d887c7a1ad498e7d24f591df599b6d9b2c76 100644
--- a/src/finiteArea/finiteArea/fac/facEdgeIntegrate.H
+++ b/src/finiteArea/finiteArea/fac/facEdgeIntegrate.H
@@ -35,7 +35,6 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-
 #ifndef facEdgeIntegrate_H
 #define facEdgeIntegrate_H
 
diff --git a/src/finiteArea/finiteArea/fac/facGrad.H b/src/finiteArea/finiteArea/fac/facGrad.H
index 29c9d4ef0c8b4d3d864be1625702014bac3d1f1b..89e60708343309da5288feacb7460ecad3b7f8e4 100644
--- a/src/finiteArea/finiteArea/fac/facGrad.H
+++ b/src/finiteArea/finiteArea/fac/facGrad.H
@@ -34,7 +34,6 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-
 #ifndef facGrad_H
 #define facGrad_H
 
diff --git a/src/finiteArea/finiteArea/fac/facLaplacian.H b/src/finiteArea/finiteArea/fac/facLaplacian.H
index 0e29d7c97dbdb274a6c5845f222f21e6afa57c98..53697dcdce6a9461904fab99a9b53c7568824a2b 100644
--- a/src/finiteArea/finiteArea/fac/facLaplacian.H
+++ b/src/finiteArea/finiteArea/fac/facLaplacian.H
@@ -34,7 +34,6 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-
 #ifndef facLaplacian_H
 #define facLaplacian_H
 
diff --git a/src/finiteArea/finiteArea/fac/facLnGrad.H b/src/finiteArea/finiteArea/fac/facLnGrad.H
index 8a136ed844aab1d6a7eeba5bfb4e7bca9c0bc6a7..cfc3cbb11182a64dafbc14e0d9bfe31693fcf993 100644
--- a/src/finiteArea/finiteArea/fac/facLnGrad.H
+++ b/src/finiteArea/finiteArea/fac/facLnGrad.H
@@ -34,7 +34,6 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-
 #ifndef facLnGrad_H
 #define facLnGrad_H
 
diff --git a/src/finiteArea/finiteArea/gradSchemes/faGradScheme/faGradScheme.H b/src/finiteArea/finiteArea/gradSchemes/faGradScheme/faGradScheme.H
index 11696af03bdf98ea7a8da6998d98b13f4e9e9443..e016fc31ed73178266c059b53ccd7107e00daa98 100644
--- a/src/finiteArea/finiteArea/gradSchemes/faGradScheme/faGradScheme.H
+++ b/src/finiteArea/finiteArea/gradSchemes/faGradScheme/faGradScheme.H
@@ -49,6 +49,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class faMesh;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -65,8 +66,9 @@ class gradScheme
 :
     public refCount
 {
-    // Private data
+    // Private Data
 
+        //- Reference to the mesh
         const faMesh& mesh_;
 
 
@@ -119,10 +121,7 @@ public:
     // Member Functions
 
         //- Return mesh reference
-        const faMesh& mesh() const
-        {
-            return mesh_;
-        }
+        const faMesh& mesh() const noexcept { return mesh_; }
 
         //- Calculate and return the grad of the given field
         virtual tmp
diff --git a/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.H b/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.H
index 84303d3fa5bf5be505e909cbb43b4087783f5462..7f4601b57cbd41d88da2757344dfae8a1720234b 100644
--- a/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.H
+++ b/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.H
@@ -62,8 +62,9 @@ class gaussGrad
 :
     public fa::gradScheme<Type>
 {
-    // Private data
+    // Private Data
 
+        //- Edge-interpolation scheme
         tmp<edgeInterpolationScheme<Type>> tinterpScheme_;
 
 
diff --git a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C
index 5893faa9102a8ff0c79089a9b2260f9eded23b09..daff764a36128cb2ab29f4357c726ba95a0825c9 100644
--- a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C
+++ b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C
@@ -63,22 +63,19 @@ leastSquaresFaGrad<Type>::grad
 
     const faMesh& mesh = vsf.mesh();
 
-    tmp<GeometricField<GradType, faPatchField, areaMesh>> tlsGrad
+    auto tlsGrad = tmp<GeometricField<GradType, faPatchField, areaMesh>>::New
     (
-        new GeometricField<GradType, faPatchField, areaMesh>
+        IOobject
         (
-            IOobject
-            (
-                "grad(" + vsf.name() + ')',
-                vsf.instance(),
-                vsf.db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensioned<GradType>(vsf.dimensions()/dimLength, Zero),
-            zeroGradientFaPatchField<GradType>::typeName
-        )
+            "grad(" + vsf.name() + ')',
+            vsf.instance(),
+            vsf.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensioned<GradType>(vsf.dimensions()/dimLength, Zero),
+        zeroGradientFaPatchField<GradType>::typeName
     );
     GeometricField<GradType, faPatchField, areaMesh>& lsGrad = tlsGrad.ref();
 
@@ -93,10 +90,10 @@ leastSquaresFaGrad<Type>::grad
 
     forAll(own, edgei)
     {
-        label ownEdgeI = own[edgei];
-        label neiEdgeI = nei[edgei];
+        const label ownEdgeI = own[edgei];
+        const label neiEdgeI = nei[edgei];
 
-        Type deltaVsf = vsf[neiEdgeI] - vsf[ownEdgeI];
+        const Type deltaVsf(vsf[neiEdgeI] - vsf[ownEdgeI]);
 
         lsGrad[ownEdgeI] += ownLs[edgei]*deltaVsf;
         lsGrad[neiEdgeI] -= neiLs[edgei]*deltaVsf;
@@ -112,7 +109,7 @@ leastSquaresFaGrad<Type>::grad
 
         if (vsf.boundaryField()[patchi].coupled())
         {
-            Field<Type> neiVsf
+            const Field<Type> neiVsf
             (
                 vsf.boundaryField()[patchi].patchNeighbourField()
             );
diff --git a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C
index 94ea4d5975a02bfe108737ec8877c940193109ce..614ae26a3b4da0e7b99fc0099ffaab9a0ad754a4 100644
--- a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C
+++ b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C
@@ -117,10 +117,10 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
 
     forAll(owner, facei)
     {
-        label own = owner[facei];
-        label nei = neighbour[facei];
+        const label own = owner[facei];
+        const label nei = neighbour[facei];
 
-        vector d = C[nei] - C[own];
+        vector d(C[nei] - C[own]);
 
         // Do not allow any mag(val) < SMALL
         if (mag(d) < SMALL)
@@ -128,7 +128,7 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
             d = vector::uniform(SMALL);
         }
 
-        symmTensor wdd = (1.0/magSqr(d))*sqr(d);
+        const symmTensor wdd((1.0/magSqr(d))*sqr(d));
 
         dd[own] += wdd;
         dd[nei] += wdd;
@@ -164,10 +164,10 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
     // Revisit all faces and calculate the lsP and lsN vectors
     forAll(owner, facei)
     {
-        label own = owner[facei];
-        label nei = neighbour[facei];
+        const label own = owner[facei];
+        const label nei = neighbour[facei];
 
-        vector d = C[nei] - C[own];
+        vector d(C[nei] - C[own]);
 
         // Do not allow any mag(val) < SMALL
         if (mag(d) < SMALL)
@@ -175,7 +175,7 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
             d = vector::uniform(SMALL);
         }
 
-        scalar magSfByMagSqrd = 1.0/magSqr(d);
+        const scalar magSfByMagSqrd = 1.0/magSqr(d);
 
         lsP[facei] = magSfByMagSqrd*(invDd[own] & d);
         lsN[facei] = -magSfByMagSqrd*(invDd[nei] & d);
diff --git a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrad.H b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrad.H
index 8705ecafa544233b7d89f3dca18c97f34977231c..fd9f21d5754f26d3be135c93f0316961526d4dc4 100644
--- a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrad.H
+++ b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrad.H
@@ -65,6 +65,7 @@ class edgeLimitedGrad
 {
     // Private Data
 
+        //- Basic gradient scheme
         tmp<fa::gradScheme<Type>> basicGradScheme_;
 
         //- Limiter coefficient
diff --git a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C
index aa6cf862f96310f160411676a0fc47e489ff830f..3f88d529c1a8a9a9137f18717b883921115d4c8e 100644
--- a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C
+++ b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C
@@ -97,19 +97,19 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::grad
     // create limiter
     scalarField limiter(vsf.internalField().size(), 1.0);
 
-    scalar rk = (1.0/k_ - 1.0);
+    const scalar rk = (1.0/k_ - 1.0);
 
     forAll(owner, edgei)
     {
-        label own = owner[edgei];
-        label nei = neighbour[edgei];
+        const label own = owner[edgei];
+        const label nei = neighbour[edgei];
 
-        scalar vsfOwn = vsf[own];
-        scalar vsfNei = vsf[nei];
+        const scalar vsfOwn = vsf[own];
+        const scalar vsfNei = vsf[nei];
 
         scalar maxEdge = max(vsfOwn, vsfNei);
         scalar minEdge = min(vsfOwn, vsfNei);
-        scalar maxMinEdge = rk*(maxEdge - minEdge);
+        const scalar maxMinEdge = rk*(maxEdge - minEdge);
         maxEdge += maxMinEdge;
         minEdge -= maxMinEdge;
 
@@ -145,14 +145,14 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::grad
 
             forAll(pOwner, pEdgei)
             {
-                label own = pOwner[pEdgei];
+                const label own = pOwner[pEdgei];
 
-                scalar vsfOwn = vsf[own];
-                scalar vsfNei = psfNei[pEdgei];
+                const scalar vsfOwn = vsf[own];
+                const scalar vsfNei = psfNei[pEdgei];
 
                 scalar maxEdge = max(vsfOwn, vsfNei);
                 scalar minEdge = min(vsfOwn, vsfNei);
-                scalar maxMinEdge = rk*(maxEdge - minEdge);
+                const scalar maxMinEdge = rk*(maxEdge - minEdge);
                 maxEdge += maxMinEdge;
                 minEdge -= maxMinEdge;
 
@@ -168,14 +168,14 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::grad
         {
             forAll(pOwner, pEdgei)
             {
-                label own = pOwner[pEdgei];
+                const label own = pOwner[pEdgei];
 
-                scalar vsfOwn = vsf[own];
-                scalar vsfNei = psf[pEdgei];
+                const scalar vsfOwn = vsf[own];
+                const scalar vsfNei = psf[pEdgei];
 
                 scalar maxEdge = max(vsfOwn, vsfNei);
                 scalar minEdge = min(vsfOwn, vsfNei);
-                scalar maxMinEdge = rk*(maxEdge - minEdge);
+                const scalar maxMinEdge = rk*(maxEdge - minEdge);
                 maxEdge += maxMinEdge;
                 minEdge -= maxMinEdge;
 
@@ -231,25 +231,25 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::grad
     // create limiter
     scalarField limiter(vvf.internalField().size(), 1.0);
 
-    scalar rk = (1.0/k_ - 1.0);
+    const scalar rk = (1.0/k_ - 1.0);
 
     forAll(owner, edgei)
     {
-        label own = owner[edgei];
-        label nei = neighbour[edgei];
+        const label own = owner[edgei];
+        const label nei = neighbour[edgei];
 
-        vector vvfOwn = vvf[own];
-        vector vvfNei = vvf[nei];
+        const vector vvfOwn(vvf[own]);
+        const vector vvfNei(vvf[nei]);
 
         // owner side
-        vector gradf = (Cf[edgei] - C[own]) & g[own];
+        vector gradf((Cf[edgei] - C[own]) & g[own]);
 
         scalar vsfOwn = gradf & vvfOwn;
         scalar vsfNei = gradf & vvfNei;
 
         scalar maxEdge = max(vsfOwn, vsfNei);
         scalar minEdge = min(vsfOwn, vsfNei);
-        scalar maxMinEdge = rk*(maxEdge - minEdge);
+        const scalar maxMinEdge = rk*(maxEdge - minEdge);
         maxEdge += maxMinEdge;
         minEdge -= maxMinEdge;
 
@@ -294,19 +294,19 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::grad
 
             forAll(pOwner, pEdgei)
             {
-                label own = pOwner[pEdgei];
+                const label own = pOwner[pEdgei];
 
-                vector vvfOwn = vvf[own];
-                vector vvfNei = psfNei[pEdgei];
+                const vector vvfOwn(vvf[own]);
+                const vector vvfNei(psfNei[pEdgei]);
 
-                vector gradf = (pCf[pEdgei] - C[own]) & g[own];
+                const vector gradf((pCf[pEdgei] - C[own]) & g[own]);
 
-                scalar vsfOwn = gradf & vvfOwn;
-                scalar vsfNei = gradf & vvfNei;
+                const scalar vsfOwn = gradf & vvfOwn;
+                const scalar vsfNei = gradf & vvfNei;
 
                 scalar maxEdge = max(vsfOwn, vsfNei);
                 scalar minEdge = min(vsfOwn, vsfNei);
-                scalar maxMinEdge = rk*(maxEdge - minEdge);
+                const scalar maxMinEdge = rk*(maxEdge - minEdge);
                 maxEdge += maxMinEdge;
                 minEdge -= maxMinEdge;
 
@@ -322,19 +322,19 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::grad
         {
             forAll(pOwner, pEdgei)
             {
-                label own = pOwner[pEdgei];
+                const label own = pOwner[pEdgei];
 
-                vector vvfOwn = vvf[own];
-                vector vvfNei = psf[pEdgei];
+                const vector vvfOwn(vvf[own]);
+                const vector vvfNei(psf[pEdgei]);
 
-                vector gradf = (pCf[pEdgei] - C[own]) & g[own];
+                const vector gradf((pCf[pEdgei] - C[own]) & g[own]);
 
-                scalar vsfOwn = gradf & vvfOwn;
-                scalar vsfNei = gradf & vvfNei;
+                const scalar vsfOwn = gradf & vvfOwn;
+                const scalar vsfNei = gradf & vvfNei;
 
                 scalar maxEdge = max(vsfOwn, vsfNei);
                 scalar minEdge = min(vsfOwn, vsfNei);
-                scalar maxMinEdge = rk*(maxEdge - minEdge);
+                const scalar maxMinEdge = rk*(maxEdge - minEdge);
                 maxEdge += maxMinEdge;
                 minEdge -= maxMinEdge;
 
diff --git a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrad.H b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrad.H
index a8a0bcbe9674f84bc0682e14dcc4dbfaf75e83f2..60a87cc48d989551ae720236a17622ca629e2a69 100644
--- a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrad.H
+++ b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrad.H
@@ -68,6 +68,7 @@ class faceLimitedGrad
 {
     // Private Data
 
+        //- Basic gradient scheme
         tmp<fa::gradScheme<Type>> basicGradScheme_;
 
         //- Limiter coefficient
diff --git a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C
index 38bd27711145660def2e934451f153fb03866582..c24cbcf0bd264033c7861a08ee301f94394e00ee 100644
--- a/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C
+++ b/src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C
@@ -121,11 +121,11 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
 
     forAll(owner, facei)
     {
-        label own = owner[facei];
-        label nei = neighbour[facei];
+        const label own = owner[facei];
+        const label nei = neighbour[facei];
 
-        scalar vsfOwn = vsf[own];
-        scalar vsfNei = vsf[nei];
+        const scalar vsfOwn = vsf[own];
+        const scalar vsfNei = vsf[nei];
 
         maxVsf[own] = max(maxVsf[own], vsfNei);
         minVsf[own] = min(minVsf[own], vsfNei);
@@ -149,8 +149,8 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
 
             forAll(pOwner, pFacei)
             {
-                label own = pOwner[pFacei];
-                scalar vsfNei = psfNei[pFacei];
+                const label own = pOwner[pFacei];
+                const scalar vsfNei = psfNei[pFacei];
 
                 maxVsf[own] = max(maxVsf[own], vsfNei);
                 minVsf[own] = min(minVsf[own], vsfNei);
@@ -160,8 +160,8 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
         {
             forAll(pOwner, pFacei)
             {
-                label own = pOwner[pFacei];
-                scalar vsfNei = psf[pFacei];
+                const label own = pOwner[pFacei];
+                const scalar vsfNei = psf[pFacei];
 
                 maxVsf[own] = max(maxVsf[own], vsfNei);
                 minVsf[own] = min(minVsf[own], vsfNei);
@@ -174,7 +174,7 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
 
     if (k_ < 1.0)
     {
-        scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
+        const scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
         maxVsf += maxMinVsf;
         minVsf -= maxMinVsf;
     }
@@ -185,8 +185,8 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
 
     forAll(owner, facei)
     {
-        label own = owner[facei];
-        label nei = neighbour[facei];
+        const label own = owner[facei];
+        const label nei = neighbour[facei];
 
         // owner side
         limitEdge
@@ -214,7 +214,7 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
 
         forAll(pOwner, pFacei)
         {
-            label own = pOwner[pFacei];
+            const label own = pOwner[pFacei];
 
             limitEdge
             (
@@ -270,8 +270,8 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
 
     forAll(owner, facei)
     {
-        label own = owner[facei];
-        label nei = neighbour[facei];
+        const label own = owner[facei];
+        const label nei = neighbour[facei];
 
         const vector& vsfOwn = vsf[own];
         const vector& vsfNei = vsf[nei];
@@ -293,11 +293,11 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
 
         if (psf.coupled())
         {
-            vectorField psfNei(psf.patchNeighbourField());
+            const vectorField psfNei(psf.patchNeighbourField());
 
             forAll(pOwner, pFacei)
             {
-                label own = pOwner[pFacei];
+                const label own = pOwner[pFacei];
                 const vector& vsfNei = psfNei[pFacei];
 
                 maxVsf[own] = max(maxVsf[own], vsfNei);
@@ -308,7 +308,7 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
         {
             forAll(pOwner, pFacei)
             {
-                label own = pOwner[pFacei];
+                const label own = pOwner[pFacei];
                 const vector& vsfNei = psf[pFacei];
 
                 maxVsf[own] = max(maxVsf[own], vsfNei);
@@ -322,7 +322,7 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
 
     if (k_ < 1.0)
     {
-        vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
+        const vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
         maxVsf += maxMinVsf;
         minVsf -= maxMinVsf;
 
@@ -336,8 +336,8 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
 
     forAll(owner, facei)
     {
-        label own = owner[facei];
-        label nei = neighbour[facei];
+        const label own = owner[facei];
+        const label nei = neighbour[facei];
 
         // owner side
         limitEdge
@@ -365,7 +365,7 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
 
         forAll(pOwner, pFacei)
         {
-            label own = pOwner[pFacei];
+            const label own = pOwner[pFacei];
 
             limitEdge
             (
diff --git a/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.C b/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.C
index 6a3c84bec0f92045b818047670feeea7c51a3f57..18c2b1a336653f8724db6e1e9f288f055573df36 100644
--- a/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.C
+++ b/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.C
@@ -59,19 +59,16 @@ correctedLnGrad<Type>::correction
 {
     const faMesh& mesh = this->mesh();
 
-    tmp<GeometricField<Type, faePatchField, edgeMesh>> tssf
+    auto tssf = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
     (
-        new GeometricField<Type, faePatchField, edgeMesh>
+        IOobject
         (
-            IOobject
-            (
-                "lnGradCorr("+vf.name()+')',
-                vf.instance(),
-                vf.db()
-            ),
-            mesh,
-            vf.dimensions()*mesh.deltaCoeffs().dimensions()
-        )
+            "lnGradCorr("+vf.name()+')',
+            vf.instance(),
+            vf.db()
+        ),
+        mesh,
+        vf.dimensions()*mesh.deltaCoeffs().dimensions()
     );
     GeometricField<Type, faePatchField, edgeMesh>& ssf = tssf.ref();
 
diff --git a/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.H b/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.H
index b7128c00a6cb91b6f9a79f6ff09941cf14099385..fbc082fa5fe7449dec18e18a7bc6dadc5fd2731f 100644
--- a/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.H
+++ b/src/finiteArea/finiteArea/lnGradSchemes/correctedLnGrad/correctedLnGrad.H
@@ -102,7 +102,7 @@ public:
         }
 
         //- Return true if this scheme uses an explicit correction
-        virtual bool corrected() const
+        virtual bool corrected() const noexcept
         {
             return !this->mesh().orthogonal();
         }
diff --git a/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.C b/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.C
index b8cef6777752dd7a6ee573cd9466558b33fdb90d..243d6fc38a03708a2c84ac0169e340213104e526 100644
--- a/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.C
+++ b/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.C
@@ -60,23 +60,20 @@ fourthLnGrad<Type>::correction
 {
     const faMesh& mesh = this->mesh();
 
-    tmp<GeometricField<Type, faePatchField, edgeMesh>> tcorr
+    auto tcorr = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
     (
-        new GeometricField<Type, faePatchField, edgeMesh>
+        IOobject
         (
-            IOobject
-            (
-                "lnGradCorr("+vf.name()+')',
-                vf.instance(),
-                vf.db()
-            ),
-            mesh,
-            vf.dimensions()*this->mesh().deltaCoeffs().dimensions()
-        )
+            "lnGradCorr("+vf.name()+')',
+            vf.instance(),
+            vf.db()
+        ),
+        mesh,
+        vf.dimensions()*this->mesh().deltaCoeffs().dimensions()
     );
     GeometricField<Type, faePatchField, edgeMesh>& corr = tcorr.ref();
 
-    edgeVectorField m(mesh.Le()/mesh.magLe());
+    const edgeVectorField m(mesh.Le()/mesh.magLe());
 
     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; ++cmpt)
     {
diff --git a/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.H b/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.H
index 97a250e9c6e672a27805b3b514a56281199d0abb..9d41df0b0060b3e354f70b75992265e95b43bd40 100644
--- a/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.H
+++ b/src/finiteArea/finiteArea/lnGradSchemes/fourthLnGrad/fourthLnGrad.H
@@ -102,10 +102,7 @@ public:
         }
 
         //- Return true if this scheme uses an explicit correction
-        virtual bool corrected() const
-        {
-            return true;
-        }
+        virtual bool corrected() const noexcept { return true; }
 
         //- Return the explicit correction to the fourthLnGrad
         //  for the given field
diff --git a/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.C b/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.C
index 8c3f93a3f6f064b38428f2d0baa1051084854f86..9b4b37afeeef47fa8063b5f51220ae88b36cd6d2 100644
--- a/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.C
+++ b/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.C
@@ -70,17 +70,17 @@ limitedLnGrad<Type>::correction
            *mag(lnGradScheme<Type>::lnGrad(vf, deltaCoeffs(vf), "orthSnGrad"))
            /(
                 (1 - limitCoeff_)*mag(corr)
-              + dimensionedScalar("small", corr.dimensions(), SMALL)
+              + dimensionedScalar(corr.dimensions(), SMALL)
             ),
-            dimensionedScalar("one", dimless, 1.0)
+            dimensionedScalar(dimless, 1.0)
         )
     );
 
     if (fa::debug)
     {
-        Info<< "limitedLnGrad :: limiter min: " << min(limiter.internalField())
-            << " max: "<< max(limiter.internalField())
-            << " avg: " << average(limiter.internalField()) << endl;
+        Info<< "limitedLnGrad :: limiter min: " << gMin(limiter.internalField())
+            << " max: "<< gMax(limiter.internalField())
+            << " avg: " << gAverage(limiter.internalField()) << endl;
     }
 
     return limiter*corr;
diff --git a/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.H b/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.H
index 1e0a6d7720e34a81f4c45d959cbef19d4ca8c58e..dc89586444d992f0b903ac293b38f86567862293 100644
--- a/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.H
+++ b/src/finiteArea/finiteArea/lnGradSchemes/limitedLnGrad/limitedLnGrad.H
@@ -68,7 +68,7 @@ class limitedLnGrad
 :
     public lnGradScheme<Type>
 {
-    // Private data
+    // Private Data
 
         //- Limiter.  0 = no limiting, 1 = full limiting
         scalar limitCoeff_;
@@ -128,7 +128,7 @@ public:
         }
 
         //- Return true if this scheme uses an explicit correction
-        virtual bool corrected() const
+        virtual bool corrected() const noexcept
         {
             return !this->mesh().orthogonal();
         }
diff --git a/src/finiteArea/finiteArea/lnGradSchemes/lnGradScheme/lnGradScheme.C b/src/finiteArea/finiteArea/lnGradSchemes/lnGradScheme/lnGradScheme.C
index 2fbc4fbc81bdfcae9ca5ae23f78a43d622fcbadd..412fdf00d84806c0214705c43fcb8b023f2a46f0 100644
--- a/src/finiteArea/finiteArea/lnGradSchemes/lnGradScheme/lnGradScheme.C
+++ b/src/finiteArea/finiteArea/lnGradSchemes/lnGradScheme/lnGradScheme.C
@@ -106,22 +106,18 @@ lnGradScheme<Type>::lnGrad
 {
     const faMesh& mesh = vf.mesh();
 
-    // construct GeometricField<Type, faePatchField, edgeMesh>
-    tmp<GeometricField<Type, faePatchField, edgeMesh>> tssf
+    auto tssf = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
     (
-        new GeometricField<Type, faePatchField, edgeMesh>
+        IOobject
         (
-            IOobject
-            (
-                lnGradName + "("+vf.name()+')',
-                vf.instance(),
-                vf.db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            vf.dimensions()*tdeltaCoeffs().dimensions()
-        )
+            lnGradName + "("+vf.name()+')',
+            vf.instance(),
+            vf.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        vf.dimensions()*tdeltaCoeffs().dimensions()
     );
     GeometricField<Type, faePatchField, edgeMesh>& ssf = tssf.ref();
 
@@ -132,17 +128,17 @@ lnGradScheme<Type>::lnGrad
     const labelUList& owner = mesh.owner();
     const labelUList& neighbour = mesh.neighbour();
 
-    forAll(owner, faceI)
+    forAll(owner, facei)
     {
-        ssf[faceI] =
-            deltaCoeffs[faceI]*(vf[neighbour[faceI]] - vf[owner[faceI]]);
+        ssf[facei] =
+            deltaCoeffs[facei]*(vf[neighbour[facei]] - vf[owner[facei]]);
     }
 
     auto& ssfb = ssf.boundaryFieldRef();
 
-    forAll(vf.boundaryField(), patchI)
+    forAll(vf.boundaryField(), patchi)
     {
-        ssfb[patchI] = vf.boundaryField()[patchI].snGrad();
+        ssfb[patchi] = vf.boundaryField()[patchi].snGrad();
     }
 
     return tssf;
diff --git a/src/finiteArea/finiteArea/lnGradSchemes/lnGradScheme/lnGradScheme.H b/src/finiteArea/finiteArea/lnGradSchemes/lnGradScheme/lnGradScheme.H
index 63adf7ea78a15263e6614fc3e9e7877fe9d4d7ff..3dfc2d45e1e0cc73018fe1bf91eb5954556f3096 100644
--- a/src/finiteArea/finiteArea/lnGradSchemes/lnGradScheme/lnGradScheme.H
+++ b/src/finiteArea/finiteArea/lnGradSchemes/lnGradScheme/lnGradScheme.H
@@ -48,6 +48,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class faMesh;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -64,7 +65,7 @@ class lnGradScheme
 :
     public refCount
 {
-    // Private data
+    // Private Data
 
         //- Hold reference to mesh
         const faMesh& mesh_;
@@ -116,10 +117,7 @@ public:
     // Member Functions
 
         //- Return mesh reference
-        const faMesh& mesh() const
-        {
-            return mesh_;
-        }
+        const faMesh& mesh() const noexcept { return mesh_; }
 
         //- Return the lnGrad of the given cell field
         //  with the given weighting factors
@@ -138,10 +136,7 @@ public:
         ) const = 0;
 
         //- Return true if this scheme uses an explicit correction
-        virtual bool corrected() const
-        {
-            return false;
-        }
+        virtual bool corrected() const noexcept { return false; }
 
         //- Return the explicit correction to the lnGrad
         //  for the given field
diff --git a/src/finiteArea/finiteArea/lnGradSchemes/uncorrectedLnGrad/uncorrectedLnGrad.H b/src/finiteArea/finiteArea/lnGradSchemes/uncorrectedLnGrad/uncorrectedLnGrad.H
index 8ff48692cecadb96f9ce608b8fd407be67e90096..709920e67aad49c443a5508cebd9464a09a7b53f 100644
--- a/src/finiteArea/finiteArea/lnGradSchemes/uncorrectedLnGrad/uncorrectedLnGrad.H
+++ b/src/finiteArea/finiteArea/lnGradSchemes/uncorrectedLnGrad/uncorrectedLnGrad.H
@@ -101,10 +101,7 @@ public:
         }
 
         //- Return true if this scheme uses an explicit correction
-        virtual bool corrected() const
-        {
-            return false;
-        }
+        virtual bool corrected() const noexcept { return false; }
 
         //- Return the explicit correction to the uncorrectedLnGrad
         //- for the given field
diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolate.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolate.C
index 54dcd271b1460ce156813086109b87c999756221..2fcbe116a1a743c73c4935647750dca34d1f57e5 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolate.C
+++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolate.C
@@ -100,14 +100,14 @@ Foam::fac::interpolate
     Istream& schemeData
 )
 {
-#   ifdef DEBUGInterpolations
+    #ifdef DEBUGInterpolations
     if (edgeInterpolation::debug)
     {
         InfoInFunction
             << "interpolating GeometricField<Type, faPatchField, areaMesh> "
             << endl;
     }
-#   endif
+    #endif
 
     return scheme<Type>(faceFlux, schemeData)().interpolate(vf);
 }
@@ -122,7 +122,7 @@ Foam::fac::interpolate
     const word& name
 )
 {
-#   ifdef DEBUGInterpolations
+    #ifdef DEBUGInterpolations
     if (edgeInterpolation::debug)
     {
         InfoInFunction
@@ -130,7 +130,7 @@ Foam::fac::interpolate
             << "using " << name
             << endl;
     }
-#   endif
+    #endif
 
     return scheme<Type>(faceFlux, name)().interpolate(vf);
 }
@@ -199,14 +199,14 @@ Foam::fac::interpolate
     Istream& schemeData
 )
 {
-#   ifdef DEBUGInterpolations
+    #ifdef DEBUGInterpolations
     if (edgeInterpolation::debug)
     {
         InfoInFunction
             << "interpolating GeometricField<Type, faPatchField, areaMesh> "
             << endl;
     }
-#   endif
+    #endif
 
     return scheme<Type>(vf.mesh(), schemeData)().interpolate(vf);
 }
@@ -220,7 +220,7 @@ Foam::fac::interpolate
     const word& name
 )
 {
-#   ifdef DEBUGInterpolations
+    #ifdef DEBUGInterpolations
     if (edgeInterpolation::debug)
     {
         InfoInFunction
@@ -228,7 +228,7 @@ Foam::fac::interpolate
             << "using " << name
             << endl;
     }
-#   endif
+    #endif
 
     return scheme<Type>(vf.mesh(), name)().interpolate(vf);
 }
@@ -257,7 +257,7 @@ Foam::fac::interpolate
     const GeometricField<Type, faPatchField, areaMesh>& vf
 )
 {
-#   ifdef DEBUGInterpolations
+    #ifdef DEBUGInterpolations
     if (edgeInterpolation::debug)
     {
         InfoInFunction
@@ -265,7 +265,7 @@ Foam::fac::interpolate
             << "using run-time selected scheme"
             << endl;
     }
-#   endif
+    #endif
 
     return interpolate(vf, "interpolate(" + vf.name() + ')');
 }
diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
index 39778175566321d7434270ca5ecb90d19ddacba3..3c397c4b075e7b16c1ec65ca759bb562a1eb22cf 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
+++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
@@ -31,6 +31,7 @@ License
 #include "edgeFields.H"
 #include "demandDrivenData.H"
 #include "faPatchFields.H"
+#include "unitConversion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -170,20 +171,20 @@ bool Foam::edgeInterpolation::movePoints() const
 }
 
 
-const Foam::vector& Foam::edgeInterpolation::skewCorr(const label edgeI) const
+const Foam::vector& Foam::edgeInterpolation::skewCorr(const label edgei) const
 {
     #ifdef FA_SKEW_CORRECTION
 
     return
         (
             skewCorrectionVectorsPtr_
-          ? (*skewCorrectionVectorsPtr_)[edgeI]
+          ? (*skewCorrectionVectorsPtr_)[edgei]
           : pTraits<vector>::zero
         );
 
     #else
 
-    return (*skewCorrectionVectorsPtr_)[edgeI];
+    return (*skewCorrectionVectorsPtr_)[edgei];
 
     #endif
 }
@@ -227,44 +228,44 @@ void Foam::edgeInterpolation::makeLPN() const
     // Calculate skewness correction vectors if necessary
     (void) skew();
 
-    forAll(owner, edgeI)
+    forAll(owner, edgei)
     {
-        const vector& skewCorrEdge = skewCorr(edgeI);
+        const vector& skewCorrEdge = skewCorr(edgei);
 
         scalar lPE =
             mag
             (
-                edgeCentres[edgeI]
+                edgeCentres[edgei]
               - skewCorrEdge
-              - faceCentres[owner[edgeI]]
+              - faceCentres[owner[edgei]]
             );
 
         scalar lEN =
             mag
             (
-                faceCentres[neighbour[edgeI]]
-              - edgeCentres[edgeI]
+                faceCentres[neighbour[edgei]]
+              - edgeCentres[edgei]
               + skewCorrEdge
             );
 
-        lPNIn[edgeI] = (lPE + lEN);
+        lPNIn[edgei] = (lPE + lEN);
 
         // Do not allow any mag(val) < SMALL
-        if (mag(lPNIn[edgeI]) < SMALL)
+        if (mag(lPNIn[edgei]) < SMALL)
         {
-            lPNIn[edgeI] = SMALL;
+            lPNIn[edgei] = SMALL;
         }
     }
 
 
-    forAll(lPN.boundaryField(), patchI)
+    forAll(lPN.boundaryField(), patchi)
     {
-        mesh().boundary()[patchI].makeDeltaCoeffs
+        mesh().boundary()[patchi].makeDeltaCoeffs
         (
-            lPN.boundaryFieldRef()[patchI]
+            lPN.boundaryFieldRef()[patchi]
         );
 
-        lPN.boundaryFieldRef()[patchI] = 1.0/lPN.boundaryField()[patchI];
+        lPN.boundaryFieldRef()[patchi] = 1.0/lPN.boundaryField()[patchi];
     }
 
 
@@ -313,23 +314,23 @@ void Foam::edgeInterpolation::makeWeights() const
     // Calculate skewness correction vectors if necessary
     (void) skew();
 
-    forAll(owner, edgeI)
+    forAll(owner, edgei)
     {
-        const vector& skewCorrEdge = skewCorr(edgeI);
+        const vector& skewCorrEdge = skewCorr(edgei);
 
         scalar lPE =
             mag
             (
-                edgeCentres[edgeI]
+                edgeCentres[edgei]
               - skewCorrEdge
-              - faceCentres[owner[edgeI]]
+              - faceCentres[owner[edgei]]
             );
 
         scalar lEN =
             mag
             (
-                faceCentres[neighbour[edgeI]]
-              - edgeCentres[edgeI]
+                faceCentres[neighbour[edgei]]
+              - edgeCentres[edgei]
               + skewCorrEdge
             );
 
@@ -337,15 +338,15 @@ void Foam::edgeInterpolation::makeWeights() const
         const scalar lPN = lPE + lEN;
         if (mag(lPN) > SMALL)
         {
-            weightsIn[edgeI] = lEN/lPN;
+            weightsIn[edgei] = lEN/lPN;
         }
     }
 
-    forAll(mesh().boundary(), patchI)
+    forAll(mesh().boundary(), patchi)
     {
-        mesh().boundary()[patchI].makeWeights
+        mesh().boundary()[patchi].makeWeights
         (
-            weights.boundaryFieldRef()[patchI]
+            weights.boundaryFieldRef()[patchi]
         );
     }
 
@@ -400,36 +401,36 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
     // Calculate skewness correction vectors if necessary
     (void) skew();
 
-    forAll(owner, edgeI)
+    forAll(owner, edgei)
     {
         // Edge normal - area normal
         vector edgeNormal =
-            normalised(lengths[edgeI] ^ edges[edgeI].vec(points));
+            normalised(lengths[edgei] ^ edges[edgei].vec(points));
 
         // Unit delta vector
         vector unitDelta =
-            faceCentres[neighbour[edgeI]]
-          - faceCentres[owner[edgeI]];
+            faceCentres[neighbour[edgei]]
+          - faceCentres[owner[edgei]];
 
         unitDelta.removeCollinear(edgeNormal);
         unitDelta.normalise();
 
 
-        const vector& skewCorrEdge = skewCorr(edgeI);
+        const vector& skewCorrEdge = skewCorr(edgei);
 
         scalar lPE =
             mag
             (
-                edgeCentres[edgeI]
+                edgeCentres[edgei]
               - skewCorrEdge
-              - faceCentres[owner[edgeI]]
+              - faceCentres[owner[edgei]]
             );
 
         scalar lEN =
             mag
             (
-                faceCentres[neighbour[edgeI]]
-              - edgeCentres[edgeI]
+                faceCentres[neighbour[edgei]]
+              - edgeCentres[edgei]
               + skewCorrEdge
             );
 
@@ -437,22 +438,22 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
 
 
         // Edge normal - area tangent
-        edgeNormal = normalised(lengths[edgeI]);
+        edgeNormal = normalised(lengths[edgei]);
 
         // Do not allow any mag(val) < SMALL
         const scalar alpha = lPN*(unitDelta & edgeNormal);
         if (mag(alpha) > SMALL)
         {
-            dc[edgeI] = scalar(1)/max(alpha, 0.05*lPN);
+            dc[edgei] = scalar(1)/max(alpha, 0.05*lPN);
         }
     }
 
 
-    forAll(deltaCoeffs.boundaryField(), patchI)
+    forAll(deltaCoeffs.boundaryField(), patchi)
     {
-        mesh().boundary()[patchI].makeDeltaCoeffs
+        mesh().boundary()[patchi].makeDeltaCoeffs
         (
-            deltaCoeffs.boundaryFieldRef()[patchI]
+            deltaCoeffs.boundaryFieldRef()[patchi]
         );
     }
 }
@@ -500,65 +501,67 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
 
     vectorField& CorrVecsIn = CorrVecs.primitiveFieldRef();
 
-    forAll(owner, edgeI)
+    forAll(owner, edgei)
     {
         // Edge normal - area normal
         vector edgeNormal =
-            normalised(lengths[edgeI] ^ edges[edgeI].vec(points));
+            normalised(lengths[edgei] ^ edges[edgei].vec(points));
 
         // Unit delta vector
         vector unitDelta =
-            faceCentres[neighbour[edgeI]]
-          - faceCentres[owner[edgeI]];
+            faceCentres[neighbour[edgei]]
+          - faceCentres[owner[edgei]];
 
         unitDelta.removeCollinear(edgeNormal);
         unitDelta.normalise();
 
         // Edge normal - area tangent
-        edgeNormal = normalised(lengths[edgeI]);
+        edgeNormal = normalised(lengths[edgei]);
 
         // Do not allow any mag(val) < SMALL
         const scalar alpha = unitDelta & edgeNormal;
         if (mag(alpha) > SMALL)
         {
-            deltaCoeffs[edgeI] = scalar(1)/alpha;
+            deltaCoeffs[edgei] = scalar(1)/alpha;
         }
 
         // Edge correction vector
-        CorrVecsIn[edgeI] =
+        CorrVecsIn[edgei] =
             edgeNormal
-          - deltaCoeffs[edgeI]*unitDelta;
+          - deltaCoeffs[edgei]*unitDelta;
     }
 
 
     edgeVectorField::Boundary& CorrVecsbf = CorrVecs.boundaryFieldRef();
 
-    forAll(CorrVecs.boundaryField(), patchI)
+    forAll(CorrVecs.boundaryField(), patchi)
     {
-        mesh().boundary()[patchI].makeCorrectionVectors(CorrVecsbf[patchI]);
+        mesh().boundary()[patchi].makeCorrectionVectors(CorrVecsbf[patchi]);
     }
 
-    scalar NonOrthogCoeff = 0.0;
+
+    constexpr scalar maxNonOrthoRatio = 0.1;
+    scalar nonOrthoCoeff = 0;
 
     if (owner.size() > 0)
     {
         scalarField sinAlpha(deltaCoeffs*mag(CorrVecs.internalField()));
 
-        forAll(sinAlpha, edgeI)
+        forAll(sinAlpha, edgei)
         {
-            sinAlpha[edgeI] = max(-1, min(sinAlpha[edgeI], 1));
+            sinAlpha[edgei] = max(-1, min(sinAlpha[edgei], 1));
         }
 
-        NonOrthogCoeff = max(Foam::asin(sinAlpha)*180.0/M_PI);
+        nonOrthoCoeff = max(Foam::asin(sinAlpha)*radToDeg());
     }
 
-    reduce(NonOrthogCoeff, maxOp<scalar>());
+    reduce(nonOrthoCoeff, maxOp<scalar>());
 
     DebugInFunction
-        << "non-orthogonality coefficient = " << NonOrthogCoeff << " deg."
+        << "non-orthogonality coefficient = " << nonOrthoCoeff << " deg."
         << endl;
 
-    if (NonOrthogCoeff < 0.1)
+    if (nonOrthoCoeff < maxNonOrthoRatio)
     {
         nonOrthCorrectionVectorsPtr_.reset(nullptr);
     }
@@ -610,15 +613,15 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
     const edgeList& edges = mesh().edges();
 
 
-    forAll(neighbour, edgeI)
+    forAll(neighbour, edgei)
     {
-        const vector& P = C[owner[edgeI]];
-        const vector& N = C[neighbour[edgeI]];
-        const vector& S = points[edges[edgeI].start()];
+        const vector& P = C[owner[edgei]];
+        const vector& N = C[neighbour[edgei]];
+        const vector& S = points[edges[edgei].start()];
 
         // (T:Eq. 5.4)
         const vector d(N - P);
-        const vector e(edges[edgeI].vec(points));
+        const vector e(edges[edgei].vec(points));
         const vector de(d^e);
         const scalar alpha = magSqr(de);
 
@@ -632,36 +635,36 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
         // (T:Eq. 5.3)
         const vector E(S + beta*e);
 
-        skewCorrVecs[edgeI] = Ce[edgeI] - E;
+        skewCorrVecs[edgei] = Ce[edgei] - E;
     }
 
 
     edgeVectorField::Boundary& bSkewCorrVecs =
         skewCorrVecs.boundaryFieldRef();
 
-    forAll(skewCorrVecs.boundaryField(), patchI)
+    forAll(skewCorrVecs.boundaryField(), patchi)
     {
-        faePatchVectorField& patchSkewCorrVecs = bSkewCorrVecs[patchI];
+        faePatchVectorField& patchSkewCorrVecs = bSkewCorrVecs[patchi];
 
         if (patchSkewCorrVecs.coupled())
         {
             const labelUList& edgeFaces =
-                mesh().boundary()[patchI].edgeFaces();
+                mesh().boundary()[patchi].edgeFaces();
 
             const edgeList::subList patchEdges =
-                mesh().boundary()[patchI].patchSlice(edges);
+                mesh().boundary()[patchi].patchSlice(edges);
 
-            vectorField ngbC(C.boundaryField()[patchI].patchNeighbourField());
+            vectorField ngbC(C.boundaryField()[patchi].patchNeighbourField());
 
-            forAll(patchSkewCorrVecs, edgeI)
+            forAll(patchSkewCorrVecs, edgei)
             {
-                const vector& P = C[edgeFaces[edgeI]];
-                const vector& N = ngbC[edgeI];
-                const vector& S = points[patchEdges[edgeI].start()];
+                const vector& P = C[edgeFaces[edgei]];
+                const vector& N = ngbC[edgei];
+                const vector& S = points[patchEdges[edgei].start()];
 
                 // (T:Eq. 5.4)
                 const vector d(N - P);
-                const vector e(patchEdges[edgeI].vec(points));
+                const vector e(patchEdges[edgei].vec(points));
                 const vector de(d^e);
                 const scalar alpha = magSqr(de);
 
@@ -674,8 +677,8 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
 
                 const vector E(S + beta*e);
 
-                patchSkewCorrVecs[edgeI] =
-                    Ce.boundaryField()[patchI][edgeI] - E;
+                patchSkewCorrVecs[edgei] =
+                    Ce.boundaryField()[patchi][edgei] - E;
             }
         }
     }
@@ -685,22 +688,22 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
     constexpr scalar maxSkewRatio = 0.1;
     scalar skewCoeff = 0;
 
-    forAll(own, edgeI)
+    forAll(own, edgei)
     {
-        const scalar magSkew = mag(skewCorrVecs[edgeI]);
+        const scalar magSkew = mag(skewCorrVecs[edgei]);
 
         const scalar lPN =
             mag
             (
-                Ce[edgeI]
-              - skewCorrVecs[edgeI]
-              - C[owner[edgeI]]
+                Ce[edgei]
+              - skewCorrVecs[edgei]
+              - C[owner[edgei]]
             )
           + mag
             (
-                C[neighbour[edgeI]]
-              - Ce[edgeI]
-              + skewCorrVecs[edgeI]
+                C[neighbour[edgei]]
+              - Ce[edgei]
+              + skewCorrVecs[edgei]
             );
 
         const scalar ratio = magSkew/lPN;
diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H
index d8fdfde511a7a7583c46c780317b65e6ad59081e..f6e2aed4248481c36dc6fc0045f8e42b33505352 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H
+++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H
@@ -142,10 +142,7 @@ public:
     // Member Functions
 
         //- Return mesh reference
-        const faMesh& mesh() const noexcept
-        {
-            return faMesh_;
-        }
+        const faMesh& mesh() const noexcept { return faMesh_; }
 
         //- Return reference to PN geodesic distance
         const edgeScalarField& lPN() const;
diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C
index b348c7ca0f2e8639c020fe1c89427f6d55dc3838..3f9419508b5ee19aeda27cf79bf375327047e55e 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C
+++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C
@@ -163,19 +163,16 @@ Foam::edgeInterpolationScheme<Type>::interpolate
     const labelUList& P = mesh.owner();
     const labelUList& N = mesh.neighbour();
 
-    tmp<GeometricField<Type, faePatchField, edgeMesh>> tsf
+    auto tsf = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
     (
-        new GeometricField<Type, faePatchField, edgeMesh>
+        IOobject
         (
-            IOobject
-            (
-                "interpolate("+vf.name()+')',
-                vf.instance(),
-                vf.db()
-            ),
-            mesh,
-            vf.dimensions()
-        )
+            "interpolate("+vf.name()+')',
+            vf.instance(),
+            vf.db()
+        ),
+        mesh,
+        vf.dimensions()
     );
     GeometricField<Type, faePatchField, edgeMesh>& sf = tsf.ref();
 
@@ -208,8 +205,8 @@ Foam::edgeInterpolationScheme<Type>::interpolate
 
         if (vf.boundaryField()[pi].coupled())
         {
-            label size = vf.boundaryField()[pi].patch().size();
-            label start = vf.boundaryField()[pi].patch().start();
+            const label size = vf.boundaryField()[pi].patch().size();
+            const label start = vf.boundaryField()[pi].patch().start();
 
             Field<Type> pOwnVf = vf.boundaryField()[pi].patchInternalField();
             Field<Type> pNgbVf = vf.boundaryField()[pi].patchNeighbourField();
@@ -233,10 +230,6 @@ Foam::edgeInterpolationScheme<Type>::interpolate
                       + pY[i]*transform(TN, pNgbVf[i])
                     );
             }
-
-//             sf.boundaryFieldRef()[pi] =
-//                 pLambda*vf.boundaryField()[pi].patchInternalField()
-//               + pY*vf.boundaryField()[pi].patchNeighbourField();
         }
         else
         {
@@ -279,19 +272,16 @@ Foam::edgeInterpolationScheme<Type>::interpolate
     const labelUList& P = mesh.owner();
     const labelUList& N = mesh.neighbour();
 
-    tmp<GeometricField<Type, faePatchField, edgeMesh>> tsf
+    auto tsf = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
     (
-        new GeometricField<Type, faePatchField, edgeMesh>
+        IOobject
         (
-            IOobject
-            (
-                "interpolate("+vf.name()+')',
-                vf.instance(),
-                vf.db()
-            ),
-            mesh,
-            vf.dimensions()
-        )
+            "interpolate("+vf.name()+')',
+            vf.instance(),
+            vf.db()
+        ),
+        mesh,
+        vf.dimensions()
     );
     GeometricField<Type, faePatchField, edgeMesh>& sf = tsf.ref();
 
@@ -325,8 +315,8 @@ Foam::edgeInterpolationScheme<Type>::interpolate
 
         if (vf.boundaryField()[pi].coupled())
         {
-            label size = vfb[pi].patch().size();
-            label start = vfb[pi].patch().start();
+            const label size = vfb[pi].patch().size();
+            const label start = vfb[pi].patch().start();
 
             Field<Type> pOwnVf(vfb[pi].patchInternalField());
             Field<Type> pNgbVf(vfb[pi].patchNeighbourField());
@@ -350,10 +340,6 @@ Foam::edgeInterpolationScheme<Type>::interpolate
                       + (1 - pLambda[i])*transform(TN, pNgbVf[i])
                     );
             }
-
-//             tsf().boundaryFieldRef()[pi] =
-//                 pLambda*vf.boundaryField()[pi].patchInternalField()
-//              + (1 - pLambda)*vf.boundaryField()[pi].patchNeighbourField();
         }
         else
         {
@@ -395,19 +381,16 @@ Foam::edgeInterpolationScheme<Type>::euclidianInterpolate
     const labelUList& P = mesh.owner();
     const labelUList& N = mesh.neighbour();
 
-    tmp<GeometricField<Type, faePatchField, edgeMesh>> tsf
+    auto tsf = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
     (
-        new GeometricField<Type, faePatchField, edgeMesh>
+        IOobject
         (
-            IOobject
-            (
-                "interpolate("+vf.name()+')',
-                vf.instance(),
-                vf.db()
-            ),
-            mesh,
-            vf.dimensions()
-        )
+            "interpolate("+vf.name()+')',
+            vf.instance(),
+            vf.db()
+        ),
+        mesh,
+        vf.dimensions()
     );
     GeometricField<Type, faePatchField, edgeMesh>& sf = tsf.ref();
 
diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.H b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.H
index 64bc2a9cef82885f52ff49a0f0f99fd650418a47..4936d2a207adcf10cb11918227d288120351bdee 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.H
+++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.H
@@ -48,6 +48,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class faMesh;
 
 /*---------------------------------------------------------------------------*\
@@ -59,9 +60,9 @@ class edgeInterpolationScheme
 :
     public refCount
 {
-    // Private data
+    // Private Data
 
-        //- Hold reference to mesh
+        //- Reference to the mesh
         const faMesh& mesh_;
 
 
@@ -136,10 +137,7 @@ public:
     // Member Functions
 
         //- Return mesh reference
-        const faMesh& mesh() const
-        {
-            return mesh_;
-        }
+        const faMesh& mesh() const noexcept { return mesh_; }
 
 
         //- Return the face-interpolate of the given cell field
diff --git a/src/finiteArea/interpolation/edgeInterpolation/schemes/Gamma/Gamma.H b/src/finiteArea/interpolation/edgeInterpolation/schemes/Gamma/Gamma.H
index c4bd35163959a804a34d2f8d28ff74b2efb64608..5e15274dd9c3ebb81f7bed3cf6360ec61277386c 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/schemes/Gamma/Gamma.H
+++ b/src/finiteArea/interpolation/edgeInterpolation/schemes/Gamma/Gamma.H
@@ -54,7 +54,11 @@ namespace Foam
 
 class GammaWeight
 {
-    scalar k_;
+    // Private Data
+
+        //- Model coefficient [0,1]
+        scalar k_;
+
 
 public:
 
@@ -72,7 +76,7 @@ public:
 
         // Rescale k_ to be >= 0 and <= 0.5 (TVD conformant)
         // and avoid the /0 when k_ = 0
-        k_ = max(k_/2.0, SMALL);
+        k_ = max(0.5*k_, SMALL);
     }
 
 
@@ -87,10 +91,7 @@ public:
         const vector& d
     ) const
     {
-        scalar magd = mag(d);
-        vector dHat = d/mag(d);
-
-        scalar gradf = (phiN - phiP)/magd;
+        const vector dHat(normalised(d));
 
         scalar gradcf;
         scalar udWeight;
@@ -109,8 +110,10 @@ public:
         // Stabilise for division
         gradcf = stabilise(gradcf, SMALL);
 
-        scalar phict = 1 - 0.5*gradf/gradcf;
-        scalar limiter = min(max(phict/k_, 0), 1);
+        const scalar gradf = (phiN - phiP)/mag(d);
+
+        const scalar phict = 1 - 0.5*gradf/gradcf;
+        const scalar limiter = min(max(phict/k_, 0), 1);
 
         return limiter*cdWeight + (1 - limiter)*udWeight;
     }
diff --git a/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C b/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C
index cc9e9a871cd6f492c3cf576f582346f4ce8a426e..455904beee0daae1b2032795a8a08f7b3d9c38df 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C
+++ b/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C
@@ -62,9 +62,9 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
 {
     const faMesh& mesh = this->mesh();
 
-    tmp<edgeScalarField> tWeightingFactors
+    auto tWeightingFactors = tmp<edgeScalarField>::New
     (
-        new edgeScalarField(mesh.edgeInterpolation::weights())
+        mesh.edgeInterpolation::weights()
     );
     edgeScalarField& weightingFactors = tWeightingFactors.ref();
 
@@ -73,7 +73,8 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
     tmp<areaScalarField> tvf = limiter(phi);
     const areaScalarField& vf = tvf();
 
-    const areaVectorField gradc(fac::grad(vf));
+    tmp<areaVectorField> tgradc = fac::grad(vf);
+    const areaVectorField& gradc = tgradc.cref();
 
 //     edgeVectorField d
 //     (
@@ -130,70 +131,70 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
 
     auto& bWeights = weightingFactors.boundaryFieldRef();
 
-    forAll(bWeights, patchI)
+    forAll(bWeights, patchi)
     {
-        if (bWeights[patchI].coupled())
+        if (bWeights[patchi].coupled())
         {
-            scalarField& pWeights = bWeights[patchI];
+            scalarField& pWeights = bWeights[patchi];
 
-            const scalarField& pEdgeFlux = edgeFlux_.boundaryField()[patchI];
+            const scalarField& pEdgeFlux = edgeFlux_.boundaryField()[patchi];
 
-            scalarField pVfP(vf.boundaryField()[patchI].patchInternalField());
+            scalarField pVfP(vf.boundaryField()[patchi].patchInternalField());
 
-            scalarField pVfN(vf.boundaryField()[patchI].patchNeighbourField());
+            scalarField pVfN(vf.boundaryField()[patchi].patchNeighbourField());
 
             vectorField pGradcP
             (
-                gradc.boundaryField()[patchI].patchInternalField()
+                gradc.boundaryField()[patchi].patchInternalField()
             );
 
             vectorField pGradcN
             (
-                gradc.boundaryField()[patchI].patchNeighbourField()
+                gradc.boundaryField()[patchi].patchNeighbourField()
             );
 
             vectorField CP
             (
-                mesh.areaCentres().boundaryField()[patchI].patchInternalField()
+                mesh.areaCentres().boundaryField()[patchi].patchInternalField()
             );
 
             vectorField CN
             (
-                mesh.areaCentres().boundaryField()[patchI]
+                mesh.areaCentres().boundaryField()[patchi]
                 .patchNeighbourField()
             );
 
             vectorField nP
             (
-                mesh.faceAreaNormals().boundaryField()[patchI]
+                mesh.faceAreaNormals().boundaryField()[patchi]
                .patchInternalField()
             );
 
             vectorField nN
             (
-                mesh.faceAreaNormals().boundaryField()[patchI]
+                mesh.faceAreaNormals().boundaryField()[patchi]
                .patchNeighbourField()
             );
 
             scalarField pLPN
             (
-                mesh.edgeInterpolation::lPN().boundaryField()[patchI]
+                mesh.edgeInterpolation::lPN().boundaryField()[patchi]
             );
 
-            forAll(pWeights, edgeI)
+            forAll(pWeights, edgei)
             {
-                vector d(CN[edgeI] - CP[edgeI]);
+                vector d(CN[edgei] - CP[edgei]);
 
-                if (pEdgeFlux[edgeI] > 0)
+                if (pEdgeFlux[edgei] > 0)
                 {
-                    d.removeCollinear(nP[edgeI]);
+                    d.removeCollinear(nP[edgei]);
                 }
                 else
                 {
-                    d.removeCollinear(nN[edgeI]);
+                    d.removeCollinear(nN[edgei]);
                 }
                 d.normalise();
-                d *= pLPN[edgeI];
+                d *= pLPN[edgei];
 
                 // Do not allow any mag(val) < SMALL
                 if (mag(d) < SMALL)
@@ -201,15 +202,15 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
                     d = vector::uniform(SMALL);
                 }
 
-                pWeights[edgeI] =
+                pWeights[edgei] =
                     this->weight
                     (
-                        pWeights[edgeI],
-                        pEdgeFlux[edgeI],
-                        pVfP[edgeI],
-                        pVfN[edgeI],
-                        pGradcP[edgeI],
-                        pGradcN[edgeI],
+                        pWeights[edgei],
+                        pEdgeFlux[edgei],
+                        pVfP[edgei],
+                        pVfN[edgei],
+                        pGradcP[edgei],
+                        pGradcN[edgei],
                         d
                     );
             }
diff --git a/src/finiteArea/interpolation/edgeInterpolation/schemes/blended/blendedEdgeInterpolation.H b/src/finiteArea/interpolation/edgeInterpolation/schemes/blended/blendedEdgeInterpolation.H
index 41f3636176835e27cb5bbdccc4dca3696775a5d0..bf23c0f0da65be59b1190ad344996ad566305c8c 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/schemes/blended/blendedEdgeInterpolation.H
+++ b/src/finiteArea/interpolation/edgeInterpolation/schemes/blended/blendedEdgeInterpolation.H
@@ -56,8 +56,9 @@ class blendedEdgeInterpolation
     public linearEdgeInterpolation<Type>,
     public upwindEdgeInterpolation<Type>
 {
-    // Private data
+    // Private Data
 
+        //- Blending factor
         const scalar blendingFactor_;
 
 
diff --git a/src/finiteArea/interpolation/edgeInterpolation/schemes/skewCorrected/skewCorrectedEdgeInterpolation.H b/src/finiteArea/interpolation/edgeInterpolation/schemes/skewCorrected/skewCorrectedEdgeInterpolation.H
index da63c64444d2a01ba880adde40e7eb974dfec6de..962f9800a84dd25d8e2c73cea150db96a41bda19 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/schemes/skewCorrected/skewCorrectedEdgeInterpolation.H
+++ b/src/finiteArea/interpolation/edgeInterpolation/schemes/skewCorrected/skewCorrectedEdgeInterpolation.H
@@ -60,6 +60,7 @@ class skewCorrectedEdgeInterpolation
 {
     // Private Data
 
+        //- Edge-interpolation scheme
         tmp<edgeInterpolationScheme<Type>> tScheme_;
 
 
@@ -132,9 +133,8 @@ public:
 
             const edgeVectorField& scv = mesh.skewCorrectionVectors();
 
-            tmp<GeometricField<Type, faePatchField, edgeMesh>> tsfCorr
-            (
-                new GeometricField<Type, faePatchField, edgeMesh>
+            auto tsfCorr =
+                tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
                 (
                     IOobject
                     (
@@ -144,9 +144,7 @@ public:
                     ),
                     mesh,
                     dimensioned<Type>(vf.dimensions(), Zero)
-                )
-            );
-
+                );
             GeometricField<Type, faePatchField, edgeMesh>& corr = tsfCorr.ref();
 
             for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; ++cmpt)
diff --git a/src/finiteArea/interpolation/edgeInterpolation/schemes/upwind/upwindEdgeInterpolation.H b/src/finiteArea/interpolation/edgeInterpolation/schemes/upwind/upwindEdgeInterpolation.H
index df84af7a6d8d15576da24412d5d6aa02b71244b5..508fe844d3d3abad0d9cadb7796db4cabd02fd97 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/schemes/upwind/upwindEdgeInterpolation.H
+++ b/src/finiteArea/interpolation/edgeInterpolation/schemes/upwind/upwindEdgeInterpolation.H
@@ -55,8 +55,9 @@ class upwindEdgeInterpolation
 :
     virtual public edgeInterpolationScheme<Type>
 {
-    // Private data
+    // Private Data
 
+        //- Face flux
         const edgeScalarField& faceFlux_;