From 1bc2ffad9907dea217ad03c61224d102e2123ba8 Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Tue, 3 Nov 2020 10:25:02 +0000
Subject: [PATCH] BUG: atm wall functions: fix double "value" entry issue
 (#1900)

  STYLE: atm wall functions: use auto and bool types wherever possible
  TUT: atmosphericModels: changes for style consistency
---
 .../atmBoundaryLayer/atmBoundaryLayer.C       |   2 +-
 .../atmBoundaryLayer/atmBoundaryLayer.H       |   9 +-
 ...ndaryLayerInletEpsilonFvPatchScalarField.H |  10 +-
 ...atmBoundaryLayerInletKFvPatchScalarField.H |   8 +-
 ...oundaryLayerInletOmegaFvPatchScalarField.H |   8 +-
 ...daryLayerInletVelocityFvPatchVectorField.H |  11 +-
 ...entHeatFluxTemperatureFvPatchScalarField.C |   4 +-
 ...entHeatFluxTemperatureFvPatchScalarField.H |  12 +-
 ...atmAlphatkWallFunctionFvPatchScalarField.C |  23 ++--
 ...atmAlphatkWallFunctionFvPatchScalarField.H |  10 +-
 ...atmEpsilonWallFunctionFvPatchScalarField.C |   2 +-
 ...atmEpsilonWallFunctionFvPatchScalarField.H |   8 +-
 .../atmNutUWallFunctionFvPatchScalarField.C   |  30 ++---
 .../atmNutUWallFunctionFvPatchScalarField.H   |  16 +--
 .../atmNutWallFunctionFvPatchScalarField.C    |  26 ++--
 .../atmNutWallFunctionFvPatchScalarField.H    |  12 +-
 .../atmNutkWallFunctionFvPatchScalarField.C   |  30 ++---
 .../atmNutkWallFunctionFvPatchScalarField.H   |  14 +--
 .../atmOmegaWallFunctionFvPatchScalarField.C  |   5 +-
 .../atmOmegaWallFunctionFvPatchScalarField.H  |   8 +-
 .../{Allrun.par => Allrun-parallel}           |   2 +
 .../HargreavesWright_2007/README              | 117 ------------------
 .../HargreavesWright_2007/README.md           | 112 +++++++++++++++++
 .../HargreavesWright_2007/system/controlDict  |   1 -
 .../system/decomposeParDict                   |   5 +-
 .../HargreavesWright_2007/system/fvSchemes    |  10 +-
 .../HargreavesWright_2007/system/sampleDict   |   1 -
 .../atmFlatTerrain/precursor/0.orig/T         |   2 +-
 .../atmFlatTerrain/precursor/0.orig/U         |   2 +-
 .../atmFlatTerrain/precursor/0.orig/alphat    |   2 +-
 .../atmFlatTerrain/precursor/0.orig/epsilon   |   2 +-
 .../atmFlatTerrain/precursor/0.orig/k         |   2 +-
 .../precursor/0.orig/leafAreaDensity          |   4 +-
 .../atmFlatTerrain/precursor/0.orig/nut       |   2 +-
 .../atmFlatTerrain/precursor/0.orig/omega     |   2 +-
 .../atmFlatTerrain/precursor/0.orig/p_rgh     |   2 +-
 .../atmFlatTerrain/precursor/0.orig/plantCd   |   4 +-
 .../precursor/0.orig/qPlant.template          |   4 +-
 .../precursor/constant/fvOptions.template     |   1 -
 .../atmFlatTerrain/precursor/constant/g       |   2 +-
 .../constant/turbulenceProperties.template    |   3 +-
 .../precursor/system/controlDict              |   2 +-
 .../atmFlatTerrain/precursor/system/fvSchemes |  11 +-
 .../precursor/system/fvSolution               |   1 -
 .../atmFlatTerrain/precursor/system/samples   |   2 +-
 .../atmFlatTerrain/successor/0.orig/T         |   2 +-
 .../atmFlatTerrain/successor/0.orig/U         |   2 +-
 .../atmFlatTerrain/successor/0.orig/alphat    |   2 +-
 .../atmFlatTerrain/successor/0.orig/epsilon   |   2 +-
 .../atmFlatTerrain/successor/0.orig/k         |   2 +-
 .../atmFlatTerrain/successor/0.orig/nut       |   2 +-
 .../atmFlatTerrain/successor/0.orig/omega     |   2 +-
 .../atmFlatTerrain/successor/0.orig/p_rgh     |   2 +-
 .../successor/constant/fvOptions              |   1 -
 .../atmFlatTerrain/successor/constant/g       |   1 -
 .../successor/system/controlDict              |   1 -
 .../atmFlatTerrain/successor/system/fvSchemes |  11 +-
 .../successor/system/fvSolution               |   1 -
 .../successor/system/mapFieldsDict            |   1 +
 .../atmFlatTerrain/successor/system/samples   |   2 +-
 .../atmForestStability/0.orig/T               |   2 +-
 .../atmForestStability/0.orig/U               |   2 +-
 .../atmForestStability/0.orig/alphat          |   2 +-
 .../atmForestStability/0.orig/epsilon         |   2 +-
 .../atmForestStability/0.orig/k               |   2 +-
 .../atmForestStability/0.orig/leafAreaDensity |   2 +-
 .../atmForestStability/0.orig/nut             |   2 +-
 .../atmForestStability/0.orig/omega           |   2 +-
 .../atmForestStability/0.orig/p_rgh           |   5 +-
 .../atmForestStability/0.orig/plantCd         |   4 +-
 .../atmForestStability/0.orig/qPlant.template |   4 +-
 .../{Allrun => Allrun-parallel}               |   0
 .../constant/fvOptions.template               |   5 +-
 .../atmForestStability/constant/g             |   1 -
 .../constant/turbulenceProperties.template    |   2 +-
 .../atmForestStability/system/controlDict     |   1 -
 .../atmForestStability/system/fvSchemes       |  11 +-
 .../atmForestStability/system/fvSolution      |   1 -
 .../atmForestStability/system/samples         |   2 +-
 .../system/turbulenceFields                   |   2 +-
 80 files changed, 323 insertions(+), 311 deletions(-)
 rename tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/{Allrun.par => Allrun-parallel} (92%)
 delete mode 100644 tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/README
 create mode 100644 tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/README.md
 rename tutorials/verificationAndValidation/atmosphericModels/atmForestStability/{Allrun => Allrun-parallel} (100%)

diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C
index 0826e3ef91..ed933c37a4 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C
@@ -61,7 +61,7 @@ atmBoundaryLayer::atmBoundaryLayer
     const dictionary& dict
 )
 :
-    initABL_(dict.getOrDefault<Switch>("initABL", true)),
+    initABL_(dict.getOrDefault<bool>("initABL", true)),
     kappa_
     (
         dict.getCheckOrDefault<scalar>("kappa", 0.41, scalarMinMax::ge(SMALL))
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H
index c792a6a2bc..bc3e84dcfb 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H
@@ -156,7 +156,7 @@ Usage
 
     where the entries mean:
     \table
-      Property  | Description              | Type          | Req'd | Deflt
+      Property  | Description              | Type          | Reqd  | Deflt
       flowDir   | Flow direction           | TimeFunction1<vector> | yes | -
       zDir      | Ground-normal direction  | TimeFunction1<vector> | yes | -
       Uref      | Reference mean streamwise flow speed being used in <!--
@@ -237,7 +237,7 @@ protected:
 
         //- Flag to initialise profiles with the theoretical ABL expressions,
         //- otherwise initialises by using "value" entry content
-        Switch initABL_;
+        bool initABL_;
 
 
 private:
@@ -325,7 +325,7 @@ public:
             tmp<scalarField> Ustar(const scalarField& z0) const;
 
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             void autoMap(const fvPatchFieldMapper&);
@@ -334,7 +334,7 @@ public:
             void rmap(const atmBoundaryLayer&, const labelList&);
 
 
-        // Evaluate functions
+        // Evaluation
 
             //- Return the velocity distribution for the ATM
             tmp<vectorField> U(const vectorField& pCf) const;
@@ -348,6 +348,7 @@ public:
             //- Return the specific dissipation rate distribution for the ATM
             tmp<scalarField> omega(const vectorField& pCf) const;
 
+
         //- Write
         void write(Ostream&) const;
 };
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
index d1be83af20..2641003b1f 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
@@ -32,8 +32,8 @@ Group
 
 Description
     This boundary condition provides a log-law type ground-normal inlet
-    boundary condition for the turbulent kinetic energy dissipation rate,
-    i.e. \c epsilon, for homogeneous, two-dimensional, dry-air, equilibrium
+    boundary condition for the turbulent kinetic energy dissipation rate
+    (i.e. \c epsilon) for homogeneous, two-dimensional, dry-air, equilibrium
     and neutral atmospheric boundary layer modelling.
 
     The ground-normal \c epsilon profile expression:
@@ -75,8 +75,8 @@ Usage
 
     where the entries mean:
     \table
-      Property  | Description                             | Type | Req'd | Deflt
-      type      | Type name: atmBoundaryLayerInletEpsilon | word | yes   | -
+      Property  | Description                             | Type | Reqd | Deflt
+      type      | Type name: atmBoundaryLayerInletEpsilon | word | yes  | -
     \endtable
 
     The inherited entries are elaborated in:
@@ -183,7 +183,7 @@ public:
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.H
index 3d774fa9f4..e2d8d70213 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.H
@@ -32,7 +32,7 @@ Group
 
 Description
     This boundary condition provides a log-law type ground-normal inlet
-    boundary condition for the turbulent kinetic energy, i.e. \c k,
+    boundary condition for the turbulent kinetic energy (i.e. \c k)
     for homogeneous, two-dimensional, dry-air, equilibrium and neutral
     atmospheric boundary layer modelling.
 
@@ -74,8 +74,8 @@ Usage
 
     where the entries mean:
     \table
-      Property  | Description                         | Type | Req'd | Deflt
-      type      | Type name: atmBoundaryLayerInletK   | word | yes   | -
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: atmBoundaryLayerInletK   | word | yes  | -
     \endtable
 
     The inherited entries are elaborated in:
@@ -183,7 +183,7 @@ public:
         virtual void updateCoeffs();
 
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletOmega/atmBoundaryLayerInletOmegaFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletOmega/atmBoundaryLayerInletOmegaFvPatchScalarField.H
index 7423a0f80c..ffe51f9d37 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletOmega/atmBoundaryLayerInletOmegaFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletOmega/atmBoundaryLayerInletOmegaFvPatchScalarField.H
@@ -30,7 +30,7 @@ Group
 
 Description
     This boundary condition provides a log-law type ground-normal inlet
-    boundary condition for the specific dissipation rate, i.e. \c omega,
+    boundary condition for the specific dissipation rate (i.e. \c omega)
     for homogeneous, two-dimensional, dry-air, equilibrium and neutral
     atmospheric boundary layer modelling.
 
@@ -71,8 +71,8 @@ Usage
 
     where the entries mean:
     \table
-      Property  | Description                            | Type | Req'd | Deflt
-      type      | Type name: atmBoundaryLayerInletOmega  | word | yes   | -
+      Property  | Description                            | Type | Reqd | Deflt
+      type      | Type name: atmBoundaryLayerInletOmega  | word | yes  | -
     \endtable
 
     The inherited entries are elaborated in:
@@ -180,7 +180,7 @@ public:
         virtual void updateCoeffs();
 
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
index 3afc5aad1d..ddbbd8179e 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
@@ -32,8 +32,8 @@ Group
 
 Description
     This boundary condition provides a log-law type ground-normal inlet
-    boundary condition for the streamwise component of wind velocity,
-    i.e. \c u, for homogeneous, two-dimensional, dry-air, equilibrium
+    boundary condition for the streamwise component of wind velocity
+    (i.e. \c u) for homogeneous, two-dimensional, dry-air, equilibrium
     and neutral atmospheric boundary layer modelling.
 
     The ground-normal streamwise flow speed profile expression:
@@ -78,8 +78,8 @@ Usage
 
     where the entries mean:
     \table
-      Property  | Description                            | Type | Req'd | Deflt
-      type      | Type name: atmBoundaryLayerInletVelocity | word | yes   | -
+      Property  | Description                              | Type | Reqd | Deflt
+      type      | Type name: atmBoundaryLayerInletVelocity | word | yes  | -
     \endtable
 
     The inherited entries are elaborated in:
@@ -186,7 +186,8 @@ public:
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
 
-        // Mapping functions
+
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmTurbulentHeatFluxTemperature/atmTurbulentHeatFluxTemperatureFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/atmTurbulentHeatFluxTemperature/atmTurbulentHeatFluxTemperatureFvPatchScalarField.C
index 9f2729a0a7..bd376e47d0 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmTurbulentHeatFluxTemperature/atmTurbulentHeatFluxTemperatureFvPatchScalarField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/atmTurbulentHeatFluxTemperature/atmTurbulentHeatFluxTemperatureFvPatchScalarField.C
@@ -195,7 +195,7 @@ void atmTurbulentHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
     {
         FatalErrorInFunction
             << "Cp0 = " << Cp0 << " must be positive."
-            << exit(FatalIOError);
+            << exit(FatalError);
     }
 
     const scalarField q(q_->value(t));
@@ -220,7 +220,7 @@ void atmTurbulentHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
             FatalErrorInFunction
                 << "Unknown heat source type. Valid types are: "
                 << heatSourceTypeNames << nl
-                << exit(FatalIOError);
+                << exit(FatalError);
         }
     }
 
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmTurbulentHeatFluxTemperature/atmTurbulentHeatFluxTemperatureFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/atmTurbulentHeatFluxTemperature/atmTurbulentHeatFluxTemperatureFvPatchScalarField.H
index a5afba860a..f58563d8cb 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmTurbulentHeatFluxTemperature/atmTurbulentHeatFluxTemperatureFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/atmTurbulentHeatFluxTemperature/atmTurbulentHeatFluxTemperatureFvPatchScalarField.H
@@ -31,8 +31,8 @@ Group
     grpAtmWallFunctions
 
 Description
-    This boundary condition provides a fixed heat constraint on temperature,
-    i.e. \c T, to specify temperature gradient through an input heat source
+    This boundary condition provides a fixed heat constraint on temperature
+    (i.e. \c T) to specify temperature gradient through an input heat source
     which can either be specified as absolute power [W], or as flux [W/m2].
 
     Required fields:
@@ -62,7 +62,7 @@ Usage
 
     where the entries mean:
     \table
-      Property     | Description                         | Type | Req'd | Dflt
+      Property     | Description                         | Type | Reqd  | Dflt
       heatSource   | Heat source type                    | word |  yes  | -
       alphaEff     | Name of turbulent thermal diff. field [kg/m/s] <!--
                                                      --> | word | yes   | -
@@ -80,7 +80,7 @@ Usage
     Options for the \c heatSource entry:
     \verbatim
       power  | Absolute power heat source [W]
-      flux   | Flux heat source [W/m2]
+      flux   | Flux heat source           [W/m2]
     \endverbatim
 
 SourceFiles
@@ -212,7 +212,7 @@ public:
 
     // Member Functions
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap(const fvPatchFieldMapper&);
@@ -225,7 +225,7 @@ public:
             );
 
 
-        // Evaluation functions
+        // Evaluation
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmAlphatkWallFunction/atmAlphatkWallFunctionFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmAlphatkWallFunction/atmAlphatkWallFunctionFvPatchScalarField.C
index 3c9f4790de..844dcbc379 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmAlphatkWallFunction/atmAlphatkWallFunctionFvPatchScalarField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmAlphatkWallFunction/atmAlphatkWallFunctionFvPatchScalarField.C
@@ -116,7 +116,7 @@ atmAlphatkWallFunctionFvPatchScalarField
         (
             "Cmu",
             0.09,
-            scalarMinMax::ge(0)
+            scalarMinMax::ge(SMALL)
         )
     ),
     kappa_
@@ -125,7 +125,7 @@ atmAlphatkWallFunctionFvPatchScalarField
         (
             "kappa",
             0.41,
-            scalarMinMax::ge(0)
+            scalarMinMax::ge(SMALL)
         )
     ),
     Pr_(TimeFunction1<scalar>(db().time(), "Pr", dict)),
@@ -183,14 +183,15 @@ void atmAlphatkWallFunctionFvPatchScalarField::updateCoeffs()
     const label patchi = patch().index();
 
     // Retrieve turbulence properties from model
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
-    (
-        IOobject::groupName
+    const auto& turbModel =
+        db().lookupObject<turbulenceModel>
         (
-            turbulenceModel::propertiesName,
-            internalField().group()
-        )
-    );
+            IOobject::groupName
+            (
+                turbulenceModel::propertiesName,
+                internalField().group()
+            )
+        );
 
     const scalarField& y = turbModel.y()[patchi];
 
@@ -211,7 +212,7 @@ void atmAlphatkWallFunctionFvPatchScalarField::updateCoeffs()
         FatalErrorInFunction
             << "Pr cannot be negative or zero. "
             << "Please check input Pr = " << Pr
-            << exit(FatalIOError);
+            << exit(FatalError);
     }
     #endif
 
@@ -226,7 +227,7 @@ void atmAlphatkWallFunctionFvPatchScalarField::updateCoeffs()
             FatalErrorInFunction
                 << "Elements of input surface fields can only be positive. "
                 << "Please check input fields z0 and Prt."
-                << exit(FatalIOError);
+                << exit(FatalError);
         }
     }
     #endif
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmAlphatkWallFunction/atmAlphatkWallFunctionFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmAlphatkWallFunction/atmAlphatkWallFunctionFvPatchScalarField.H
index 5838e78967..598e1f4934 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmAlphatkWallFunction/atmAlphatkWallFunctionFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmAlphatkWallFunction/atmAlphatkWallFunctionFvPatchScalarField.H
@@ -32,7 +32,7 @@ Group
 
 Description
     This boundary condition provides a wall constraint on the kinematic
-    turbulent thermal conductivity, i.e. \c alphat, for atmospheric boundary
+    turbulent thermal conductivity (i.e. \c alphat) for atmospheric boundary
     layer modelling.  It assumes a logarithmic distribution of the potential
     temperature within the first cell.
 
@@ -65,8 +65,8 @@ Usage
 
     where the entries mean:
     \table
-      Property | Description                      | Type   | Req'd     | Dflt
-      type     | Type name: atmAlphatkWallFunction | word   | yes       | -
+      Property | Description                       | Type   | Reqd     | Dflt
+      type     | Type name: atmAlphatkWallFunction | word   | yes      | -
       Pr       | Molecular Prandtl number     | TimeFunction1<scalar>  | yes | -
       Prt      | Turbulent Prandtl number     | PatchFunction1<scalar> | yes | -
       z0       | Surface roughness length [m] | PatchFunction1<scalar> | yes | -
@@ -212,13 +212,13 @@ public:
 
     // Member Functions
 
-        // Evaluation functions
+        // Evaluation
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap(const fvPatchFieldMapper&);
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.C
index 16177e1a89..c032167279 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.C
@@ -74,7 +74,7 @@ void Foam::atmEpsilonWallFunctionFvPatchScalarField::calculate
             FatalErrorInFunction
                 << "z0 field can only contain positive values. "
                 << "Please check input field z0."
-                << exit(FatalIOError);
+                << exit(FatalError);
         }
     }
     #endif
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.H
index cd0b3243f9..655c772ceb 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.H
@@ -32,8 +32,8 @@ Group
 
 Description
     This boundary condition provides a wall constraint on the turbulent kinetic
-    energy dissipation rate, i.e. \c epsilon, and the turbulent kinetic energy
-    production contribution, i.e. \c G, for atmospheric boundary layer
+    energy dissipation rate (i.e. \c epsilon) and the turbulent kinetic energy
+    production contribution (i.e. \c G) for atmospheric boundary layer
     modelling.
 
     References:
@@ -75,7 +75,7 @@ Usage
 
     where the entries mean:
     \table
-      Property | Description                       | Type   | Req'd    | Dflt
+      Property | Description                       | Type   | Reqd     | Dflt
       type     | Type name: atmEpsilonWallFunction | word   | yes      | -
       z0       | Surface roughness length [m] | PatchFunction1<scalar> | yes | -
     \endtable
@@ -206,7 +206,7 @@ public:
 
     // Member Functions
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap(const fvPatchFieldMapper&);
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.C
index 2bbfcfb541..2035a29061 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.C
@@ -43,14 +43,15 @@ tmp<scalarField> atmNutUWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
-    (
-        IOobject::groupName
+    const auto& turbModel =
+        db().lookupObject<turbulenceModel>
         (
-            turbulenceModel::propertiesName,
-            internalField().group()
-        )
-    );
+            IOobject::groupName
+            (
+                turbulenceModel::propertiesName,
+                internalField().group()
+            )
+        );
 
     const scalarField& y = turbModel.y()[patchi];
 
@@ -65,20 +66,20 @@ tmp<scalarField> atmNutUWallFunctionFvPatchScalarField::calcNut() const
     const scalarField z0(z0_->value(t));
 
     #ifdef FULLDEBUG
-    for (const auto& z : z0)
+    for (const scalar z : z0)
     {
         if (z < VSMALL)
         {
             FatalErrorInFunction
                 << "z0 field can only contain positive values. "
                 << "Please check input field z0."
-                << exit(FatalIOError);
+                << exit(FatalError);
         }
     }
     #endif
 
-    tmp<scalarField> tnutw(new scalarField(*this));
-    scalarField& nutw = tnutw.ref();
+    auto tnutw = tmp<scalarField>::New(*this);
+    auto& nutw = tnutw.ref();
 
     forAll(nutw, facei)
     {
@@ -90,7 +91,7 @@ tmp<scalarField> atmNutUWallFunctionFvPatchScalarField::calcNut() const
 
     if (boundNut_)
     {
-        nutw = max(nutw, scalar(0.0));
+        nutw = max(nutw, scalar(0));
     }
 
     return tnutw;
@@ -133,7 +134,7 @@ atmNutUWallFunctionFvPatchScalarField::atmNutUWallFunctionFvPatchScalarField
 )
 :
     nutUWallFunctionFvPatchScalarField(p, iF, dict),
-    boundNut_(dict.getOrDefault<Switch>("boundNut", true)),
+    boundNut_(dict.getOrDefault<bool>("boundNut", true)),
     z0_(PatchFunction1<scalar>::New(p.patch(), "z0", dict))
 {}
 
@@ -190,7 +191,8 @@ void atmNutUWallFunctionFvPatchScalarField::rmap
 
 void atmNutUWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
-    nutUWallFunctionFvPatchScalarField::write(os);
+    fvPatchField<scalar>::write(os);
+    nutWallFunctionFvPatchScalarField::writeLocalEntries(os);
     os.writeEntry("boundNut", boundNut_);
     z0_->writeData(os);
     writeEntry("value", os);
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.H
index 1ba8d848fc..fdd5e5d5ff 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.H
@@ -32,7 +32,7 @@ Group
 
 Description
     This boundary condition provides a wall constraint on the turbulent
-    viscosity, i.e. \c nut, based on velocity, i.e. \c U, for atmospheric
+    viscosity (i.e. \c nut) based on velocity (i.e. \c U) for atmospheric
     boundary layer modelling.  It is designed to be used in conjunction
     with the \c atmBoundaryLayerInletVelocity boundary condition.
 
@@ -52,8 +52,8 @@ Description
 
     Required fields:
     \verbatim
-      nut     | Turbulent viscosity    [m2/s]
-      U       | Velocity               [m/s]
+      nut     | Turbulent viscosity         [m2/s]
+      U       | Velocity                    [m/s]
     \endverbatim
 
 Usage
@@ -77,10 +77,10 @@ Usage
 
     where the entries mean:
     \table
-      Property | Description                      | Type   | Req'd  | Dflt
-      type     | Type name: atmNutUWallFunction   | word   | yes    | -
+      Property | Description                      | Type   | Reqd  | Dflt
+      type     | Type name: atmNutUWallFunction   | word   | yes   | -
       z0       | Surface roughness length [m] | PatchFunction1<scalar> | yes | -
-      boundNut | Flag: zero-bound nut near wall   | bool   | no     | true
+      boundNut | Flag: zero-bound nut near wall   | bool   | no    | true
     \endtable
 
     The inherited entries are elaborated in:
@@ -119,7 +119,7 @@ class atmNutUWallFunctionFvPatchScalarField
 
         //- Flag to zero-bound nut to prevent negative nut
         //- at the wall arising from negative heat fluxes
-        const Switch boundNut_;
+        const bool boundNut_;
 
         //- Surface roughness length field [m]
         autoPtr<PatchFunction1<scalar>> z0_;
@@ -204,7 +204,7 @@ public:
 
     // Member Functions
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap(const fvPatchFieldMapper&);
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.C
index e9e31e756f..5fe0246a8f 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.C
@@ -44,14 +44,15 @@ tmp<scalarField> atmNutWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
-    (
-        IOobject::groupName
+    const auto& turbModel =
+        db().lookupObject<turbulenceModel>
         (
-            turbulenceModel::propertiesName,
-            internalField().group()
-        )
-    );
+            IOobject::groupName
+            (
+                turbulenceModel::propertiesName,
+                internalField().group()
+            )
+        );
     const scalarField& y = turbModel.y()[patchi];
 
     const tmp<volScalarField> tk = turbModel.k();
@@ -60,8 +61,8 @@ tmp<scalarField> atmNutWallFunctionFvPatchScalarField::calcNut() const
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
-    tmp<scalarField> tnutw(new scalarField(*this));
-    scalarField& nutw = tnutw.ref();
+    auto tnutw = tmp<scalarField>::New(*this);
+    auto& nutw = tnutw.ref();
 
     const scalar Cmu25 = pow025(Cmu_);
 
@@ -72,14 +73,14 @@ tmp<scalarField> atmNutWallFunctionFvPatchScalarField::calcNut() const
     const scalarField z0(z0_->value(t));
 
     #ifdef FULLDEBUG
-    for (const auto& z : z0)
+    for (const scalar z : z0)
     {
         if (z < VSMALL)
         {
             FatalErrorInFunction
                 << "z0 field can only contain positive values. "
                 << "Please check input field z0."
-                << exit(FatalIOError);
+                << exit(FatalError);
         }
     }
     #endif
@@ -211,7 +212,8 @@ void atmNutWallFunctionFvPatchScalarField::rmap
 
 void atmNutWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
-    nutkWallFunctionFvPatchScalarField::write(os);
+    fvPatchField<scalar>::write(os);
+    nutWallFunctionFvPatchScalarField::writeLocalEntries(os);
     os.writeEntry("z0Min", z0Min_);
     z0_->writeData(os);
     writeEntry("value", os);
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.H
index ed898c7f80..0aaf4555b5 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.H
@@ -32,8 +32,8 @@ Group
 
 Description
     This boundary condition provides a wall constraint on the turbulent
-    viscosity, i.e. \c nut, based on the turbulent kinetic energy, \c k,
-    and velocity, i.e. \c U, for atmospheric boundary layer modelling.
+    viscosity (i.e. \c nut) based on the turbulent kinetic energy (i.e. \c k)
+    and velocity (i.e. \c U) for atmospheric boundary layer modelling.
 
     The governing equation of the boundary condition:
 
@@ -114,9 +114,9 @@ Usage
 
     where the entries mean:
     \table
-      Property | Description                       | Type   | Req'd  | Dflt
-      type     | Type name: nutAtmWallFunction     | word   | yes    | -
-      z0Min | Minimum surface roughness length [m] | scalar | yes    | -
+      Property | Description                       | Type   | Reqd  | Dflt
+      type     | Type name: nutAtmWallFunction     | word   | yes   | -
+      z0Min | Minimum surface roughness length [m] | scalar | yes   | -
       z0    | Surface roughness length [m] | PatchFunction1<scalar> | yes | -
     \endtable
 
@@ -242,7 +242,7 @@ public:
 
     // Member Functions
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap(const fvPatchFieldMapper&);
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.C
index 79d9c57a22..712da459ca 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.C
@@ -45,14 +45,15 @@ tmp<scalarField> atmNutkWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
-    (
-        IOobject::groupName
+    const auto& turbModel =
+        db().lookupObject<turbulenceModel>
         (
-            turbulenceModel::propertiesName,
-            internalField().group()
-        )
-    );
+            IOobject::groupName
+            (
+                turbulenceModel::propertiesName,
+                internalField().group()
+            )
+        );
     const scalarField& y = turbModel.y()[patchi];
 
     const tmp<volScalarField> tk = turbModel.k();
@@ -61,8 +62,8 @@ tmp<scalarField> atmNutkWallFunctionFvPatchScalarField::calcNut() const
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
-    tmp<scalarField> tnutw(new scalarField(*this));
-    scalarField& nutw = tnutw.ref();
+    auto tnutw = tmp<scalarField>::New(*this);
+    auto& nutw = tnutw.ref();
 
     const scalar Cmu25 = pow025(Cmu_);
 
@@ -70,14 +71,14 @@ tmp<scalarField> atmNutkWallFunctionFvPatchScalarField::calcNut() const
     const scalarField z0(z0_->value(t));
 
     #ifdef FULLDEBUG
-    for (const auto& z : z0)
+    for (const scalar z : z0)
     {
         if (z < VSMALL)
         {
             FatalErrorInFunction
                 << "z0 field can only contain positive values. "
                 << "Please check input field z0."
-                << exit(FatalIOError);
+                << exit(FatalError);
         }
     }
     #endif
@@ -98,7 +99,7 @@ tmp<scalarField> atmNutkWallFunctionFvPatchScalarField::calcNut() const
 
     if (boundNut_)
     {
-        nutw = max(nutw, scalar(0.0));
+        nutw = max(nutw, scalar(0));
     }
 
     return tnutw;
@@ -141,7 +142,7 @@ atmNutkWallFunctionFvPatchScalarField::atmNutkWallFunctionFvPatchScalarField
 )
 :
     nutkWallFunctionFvPatchScalarField(p, iF, dict),
-    boundNut_(dict.getOrDefault<Switch>("boundNut", false)),
+    boundNut_(dict.getOrDefault<bool>("boundNut", false)),
     z0_(PatchFunction1<scalar>::New(p.patch(), "z0", dict))
 {}
 
@@ -198,7 +199,8 @@ void atmNutkWallFunctionFvPatchScalarField::rmap
 
 void atmNutkWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
-    nutkWallFunctionFvPatchScalarField::write(os);
+    fvPatchField<scalar>::write(os);
+    nutWallFunctionFvPatchScalarField::writeLocalEntries(os);
     os.writeEntry("boundNut", boundNut_);
     z0_->writeData(os);
     writeEntry("value", os);
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.H
index aa8e5fb3dc..7d49e8a62c 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.H
@@ -33,7 +33,7 @@ Group
 
 Description
     This boundary condition provides a wall constraint on the turbulent
-    viscosity, i.e. \c nut, based on the turbulent kinetic energy, i.e. \c k,
+    viscosity (i.e. \c nut) based on the turbulent kinetic energy (i.e. \c k)
     for atmospheric boundary layer modelling.  It is designed to be used in
     conjunction with the \c atmBoundaryLayerInletVelocity boundary condition.
 
@@ -47,7 +47,7 @@ Description
     \vartable
       \u^*     | Friction velocity
       \kappa   | von Kármán constant
-      z_0      | Surface roughness length [m]
+      z_0      | Surface roughness length       [m]
       z        | Ground-normal coordinate
     \endvartable
 
@@ -94,10 +94,10 @@ Usage
 
     where the entries mean:
     \table
-      Property | Description                      | Type   | Req'd    | Dflt
-      type     | Type name: atmNutkWallFunction   | word   | yes      | -
+      Property | Description                      | Type   | Reqd    | Dflt
+      type     | Type name: atmNutkWallFunction   | word   | yes     | -
       z0       | Surface roughness length [m] | PatchFunction1<scalar> | yes | -
-      boundNut | Flag: zero-bound nut near wall   | bool   | no       | false
+      boundNut | Flag to zero-bound nut near wall | bool   | no      | false
     \endtable
 
     The inherited entries are elaborated in:
@@ -139,7 +139,7 @@ class atmNutkWallFunctionFvPatchScalarField
 
         //- Flag to zero-bound nut to prevent negative nut
         //- at the wall arising from negative heat fluxes
-        const Switch boundNut_;
+        const bool boundNut_;
 
         //- Surface roughness length field [m]
         autoPtr<PatchFunction1<scalar>> z0_;
@@ -224,7 +224,7 @@ public:
 
     // Member Functions
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap(const fvPatchFieldMapper&);
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.C
index 5a2d8409a4..6605729c98 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.C
@@ -66,14 +66,14 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::calculate
     const scalarField z0(z0_->value(t));
 
     #ifdef FULLDEBUG
-    for (const auto& z : z0)
+    for (const scalar z : z0)
     {
         if (z < VSMALL)
         {
             FatalErrorInFunction
                 << "z0 field can only contain positive values. "
                 << "Please check input field z0."
-                << exit(FatalIOError);
+                << exit(FatalError);
         }
     }
     #endif
@@ -198,7 +198,6 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::write
 {
     omegaWallFunctionFvPatchScalarField::write(os);
     z0_->writeData(os);
-    writeEntry("value", os);
 }
 
 
diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.H
index 96ad16e0f1..82553be234 100644
--- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.H
@@ -32,8 +32,8 @@ Group
 
 Description
     This boundary condition provides a wall constraint on the specific
-    dissipation rate, i.e. \c omega, and the turbulent kinetic energy
-    production contribution, i.e. \c G, for atmospheric boundary
+    dissipation rate (i.e. \c omega) and the turbulent kinetic energy
+    production contribution (i.e. \c G) for atmospheric boundary
     layer modelling.
 
     References:
@@ -74,7 +74,7 @@ Usage
 
     where the entries mean:
     \table
-      Property | Description                      | Type   | Req'd  | Dflt
+      Property | Description                      | Type   | Reqd   | Dflt
       type     | Type name: atmOmegaWallFunction  | word   | yes    | -
       z0    | Surface roughness length [m] | PatchFunction1<scalar> | yes | -
     \endtable
@@ -206,7 +206,7 @@ public:
 
     // Member Functions
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap(const fvPatchFieldMapper&);
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun.par b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun-parallel
similarity index 92%
rename from tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun.par
rename to tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun-parallel
index e4725835dc..02ee5086d3 100755
--- a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun.par
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun-parallel
@@ -18,4 +18,6 @@ runParallel renumberMesh -overwrite
 
 runParallel $(getApplication)
 
+runParallel redistributePar -reconstruct -latestTime
+
 #------------------------------------------------------------------------------
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/README b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/README
deleted file mode 100644
index f5aa98553e..0000000000
--- a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/README
+++ /dev/null
@@ -1,117 +0,0 @@
-#------------------------------------------------------------------------------
-
-Overview
-
-    "By setting appropriate profiles for wind velocity and the turbulence
-    quantities at the inlet, it is often assumed that the boundary layer will
-    be maintained up to the buildings or obstructions in the flow." (HW:p. 355).
-    However, it was quantified by (HW:p. 355) that "even in the absence of
-    obstructions, ..., the velocity and turbulence profiles decay along the
-    fetch" (HW:p. 355). It was shown by (HW:p. 355) that a set of modifications
-    were required to maintain a neutral atmospheric boundary layer throughout
-    an empty and long computational domain of a RANS computation.
-
-    Aim:
-
-        Verification of the following boundary conditions in terms of the
-        maintenance of inlet quantities downstream within a RANS computation:
-
-        - atmBoundaryLayerInletVelocity
-        - atmBoundaryLayerInletK
-        - atmBoundaryLayerInletEpsilon
-        - atmBoundaryLayerInletOmega
-
-    Benchmark (Physical phenomenon):
-
-        The benchmark is an empty fetch computational domain, steady-state
-        RANS simulation involving the following traits:
-
-        - External flow
-        - The surface layer portion of the neutral-stratified equilibrium
-        atmospheric boundary layer (no Ekman layer)
-        - Dry air
-        - Homogeneous, smooth terrain
-        - Spatiotemporal-invariant aerodynamic roughness length
-        - No displacement height
-        - Newtonian, single-phase, incompressible, non-reacting
-
-    Benchmark scenario:
-
-        - Computational domain: (HW:Fig. 1)
-        - Benchmark dataset: (HW:Fig. 6) (Obtained by the WebPlotDigitizer-4.2
-        (Rohatgi, 2019))
-
-    Resources:
-
-        Computational study (tag:HW):
-            Hargreaves, D. M., & Wright, N. G. (2007).
-            On the use of the k–ε model in commercial CFD software
-            to model the neutral atmospheric boundary layer.
-            Journal of wind engineering and
-            industrial aerodynamics, 95(5), 355-369.
-            DOI:10.1016/j.jweia.2006.08.002
-
-        Wind profile (tag:RQP):
-            Richards, P. J., Quinn, A. D., & Parker, S. (2002).
-            A 6 m cube in an atmospheric boundary layer flow-Part 2.
-            Computational solutions. Wind and structures, 5(2_3_4), 177-192.
-            DOI:10.12989/was.2002.5.2_3_4.177
-
-    Physical modelling:
-
-        - The governing equations for:
-            - Steady-state, Newtonian, single-phase, incompressible fluid flows,
-            excluding any thermal chemical, electromagnetic and scalar
-            interactions
-        - Mathematical approach for the turbulence modelling:
-            - Reynolds-averaged Navier-Stokes simulation (RANS)
-        - Turbulence closure model:
-            - kEpsilon and kOmegaSST linear eddy viscosity closure models
-        - The sets of input (HW:Table 1):
-            - Reference height, Zref = 6 [m]
-            - Aerodynamic roughness height, z0 = 0.01 [m]
-            - Displacement height, d = 0 [m]
-            - Reference mean wind speed, Uref = 10 [m/s]
-
-    Computational domain modelling:
-
-        - Rectangular prism
-        - (x1, x2, x3)
-            = (5000, 100, 500) [m]
-            = (streamwise, spanwise, ground-normal) directions
-
-    Computational domain discretisation:
-
-        - Spatial resolution:
-            - (x1, x2, x3) = (500, 5, 50) [cells]
-            - Refer to the `system/blockMeshDict` for the grading details
-        - Temporal resolution: Steady state
-
-    Equation discretisation:
-
-        - Spatial derivatives and variables:
-            - Convection: Second order
-            - Others: Second order with various limiters
-
-        - Temporal derivatives and variables: First order
-
-    Numerical boundary/initial conditions:
-
-        - Refer to `0.orig`
-
-    Pressure-velocity coupling algorithm:
-
-        - SIMPLEC
-
-    Linear solvers:
-
-        - Refer to `system/fvSolution`
-
-    Initialisation and sampling:
-
-        - No initialisation/averaging
-        - Sampling at the end of the simulation via `system/sampleDict`
-        - Refer to `system/controlDict` for further details
-
-
-#------------------------------------------------------------------------------
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/README.md b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/README.md
new file mode 100644
index 0000000000..1aa1a797b9
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/README.md
@@ -0,0 +1,112 @@
+<!------------------------------------------------------------------------- -->
+
+# Overview
+
+"By setting appropriate profiles for wind velocity and the turbulence
+quantities at the inlet, it is often assumed that the boundary layer will
+be maintained up to the buildings or obstructions in the flow." (HW:p. 355).
+However, it was quantified by (HW:p. 355) that "even in the absence of
+obstructions, ..., the velocity and turbulence profiles decay along the
+fetch" (HW:p. 355). It was shown by (HW:p. 355) that a set of modifications
+were required to maintain a neutral atmospheric boundary layer throughout
+an empty and long computational domain of a RANS computation.
+
+## Aim
+
+- Verification of the atmospheric boundary-layer boundary conditions in terms
+of the maintenance of inlet quantities downstream within a RANS computation:
+  - atmBoundaryLayerInletVelocity
+  - atmBoundaryLayerInletK
+  - atmBoundaryLayerInletEpsilon
+  - atmBoundaryLayerInletOmega
+
+## Benchmark (Physical phenomenon)
+
+- The benchmark is an empty fetch computational
+domain, steady-state RANS simulation.
+- Flow characteristics:
+  - External flow
+  - The surface layer portion of the neutral-stratified
+  equilibrium atmospheric boundary layer (no Ekman layer)
+  - Dry air
+  - Homogeneous, smooth terrain
+  - Spatiotemporal-invariant aerodynamic roughness length
+  - No displacement height
+  - Newtonian, single-phase, incompressible, non-reacting
+- Benchmark scenario:
+  - Computational domain: (HW:Fig. 1)
+  - Benchmark dataset: (HW:Fig. 6)
+  (Obtained by the WebPlotDigitizer-4.2 (Rohatgi, 2019))
+
+## Resources
+
+    Computational study (tag:HW):
+        Hargreaves, D. M., & Wright, N. G. (2007).
+        On the use of the k–ε model in commercial CFD software
+        to model the neutral atmospheric boundary layer.
+        Journal of wind engineering and
+        industrial aerodynamics, 95(5), 355-369.
+        DOI:10.1016/j.jweia.2006.08.002
+
+    Wind profile (tag:RQP):
+        Richards, P. J., Quinn, A. D., & Parker, S. (2002).
+        A 6 m cube in an atmospheric boundary layer flow-Part 2.
+        Computational solutions. Wind and structures, 5(2_3_4), 177-192.
+        DOI:10.12989/was.2002.5.2_3_4.177
+
+
+# Numerics
+## Physical modelling:
+
+- The governing equations for:
+  - Steady-state, Newtonian, single-phase, incompressible fluid flows,
+    excluding any thermal chemical, electromagnetic and scalar interactions
+- Mathematical approach for the turbulence modelling:
+  - Reynolds-averaged Navier-Stokes simulation (RANS)
+- Turbulence closure model:
+  - kEpsilon and kOmegaSST linear eddy viscosity closure models
+- The sets of input (HW:Table 1):
+  - Reference height, Zref = 6 [m]
+  - Aerodynamic roughness height, z0 = 0.01 [m]
+  - Displacement height, d = 0 [m]
+  - Reference mean wind speed, Uref = 10 [m/s]
+
+## Computational domain modelling:
+
+- Rectangular prism
+- (x1, x2, x3) = (5000, 100, 500) [m] = (streamwise, spanwise, ground-normal) directions
+
+## Computational domain discretisation:
+
+- Spatial resolution:
+  - (x1, x2, x3) = (500, 5, 50) [cells]
+  - Refer to the `system/blockMeshDict` for the grading details
+- Temporal resolution: Steady state
+
+## Equation discretisation:
+
+- Spatial derivatives and variables:
+  - Convection: Second order
+  - Others: Second order with various limiters
+- Temporal derivatives and variables: First order
+
+## Numerical boundary/initial conditions:
+
+- Refer to `0.orig`
+
+## Pressure-velocity coupling algorithm:
+
+- SIMPLEC
+
+## Linear solvers:
+
+- Refer to `system/fvSolution`
+
+## Initialisation and sampling:
+
+- No initialisation/averaging
+- Sampling at the end of the simulation via `system/sampleDict`
+- Refer to `system/controlDict` for further details
+
+
+<!------------------------------------------------------------------------- -->
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/controlDict b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/controlDict
index f161c28683..85c578e863 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/controlDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/controlDict
@@ -47,7 +47,6 @@ timePrecision   6;
 
 runTimeModifiable yes;
 
-
 functions
 {
     #include "sampleDict"
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/decomposeParDict b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/decomposeParDict
index 8c168f1475..a2ba30f384 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/decomposeParDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/decomposeParDict
@@ -16,11 +16,12 @@ FoamFile
 
 numberOfSubdomains 8;
 
-method          hierarchical;
+method             hierarchical;
 
 coeffs
 {
-    n           (8 1 1);
+    n              (8 1 1);
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/fvSchemes b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/fvSchemes
index 65911f9975..4d7c420932 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/fvSchemes
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/fvSchemes
@@ -27,10 +27,14 @@ gradSchemes
 divSchemes
 {
     default             none;
+
     div(phi,U)          bounded Gauss linear;
-    div(phi,epsilon)    bounded Gauss limitedLinear 1;
-    div(phi,omega)      bounded Gauss limitedLinear 1;
-    div(phi,k)          bounded Gauss limitedLinear 1;
+
+    turbulence          bounded Gauss limitedLinear 1;
+    div(phi,epsilon)    $turbulence;
+    div(phi,omega)      $turbulence;
+    div(phi,k)          $turbulence;
+
     div((nuEff*dev2(T(grad(U)))))    Gauss linear;
 }
 
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/sampleDict b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/sampleDict
index 93a335937c..a84fb6af7c 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/sampleDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/sampleDict
@@ -10,7 +10,6 @@ FoamFile
     version         2.0;
     format          ascii;
     class           dictionary;
-    location        system;
     object          sampleDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/T b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/T
index 7eff4a6bef..0c6aeb4f4b 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/T
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/T
@@ -36,7 +36,7 @@ boundaryField
         value          uniform 300;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/U b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/U
index 88057c922b..4152f9deee 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/U
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/U
@@ -31,7 +31,7 @@ boundaryField
         type            slip;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/alphat b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/alphat
index b65c235454..89992ade62 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/alphat
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/alphat
@@ -37,7 +37,7 @@ boundaryField
         value           uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/epsilon b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/epsilon
index 477e99b7d5..1c3b68a5a8 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/epsilon
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/epsilon
@@ -34,7 +34,7 @@ boundaryField
         type            slip;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/k b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/k
index 30b68562af..92b83b2515 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/k
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/k
@@ -31,7 +31,7 @@ boundaryField
         type            slip;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/leafAreaDensity b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/leafAreaDensity
index 2542580d22..2947c12847 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/leafAreaDensity
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/leafAreaDensity
@@ -20,13 +20,13 @@ internalField   uniform 0.0;
 
 boundaryField
 {
-    "bottom|top"
+    "(bottom|top)"
     {
         type            fixedValue;
         value           uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/nut b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/nut
index 9411874bbb..9f32db8e73 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/nut
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/nut
@@ -35,7 +35,7 @@ boundaryField
         value           uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/omega b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/omega
index 33528c09fa..26f8ee089a 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/omega
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/omega
@@ -34,7 +34,7 @@ boundaryField
         type            slip;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/p_rgh b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/p_rgh
index 7042c43d7d..9df254bcec 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/p_rgh
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/p_rgh
@@ -27,7 +27,7 @@ boundaryField
        value            uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/plantCd b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/plantCd
index 68243f1aa0..70fb206c50 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/plantCd
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/plantCd
@@ -20,13 +20,13 @@ internalField   uniform 0.0;
 
 boundaryField
 {
-    "bottom|top"
+    "(bottom|top)"
     {
         type            fixedValue;
         value           uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/qPlant.template b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/qPlant.template
index 7a57452fa9..f8ebebffe9 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/qPlant.template
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/0.orig/qPlant.template
@@ -20,13 +20,13 @@ internalField  uniform  Q_PLANT;
 
 boundaryField
 {
-    "bottom|top"
+    "(bottom|top)"
     {
         type            fixedValue;
         value           uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/fvOptions.template b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/fvOptions.template
index 9f9643c411..a7929ff194 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/fvOptions.template
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/fvOptions.template
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/g b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/g
index 97b92bccc3..5d28f81891 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/g
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/g
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       uniformDimensionedVectorField;
-    location    "constant";
     object      g;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -18,4 +17,5 @@ FoamFile
 dimensions      [0 1 -2 0 0 0 0];
 value           (0 0 -9.81);
 
+
 // ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/turbulenceProperties.template b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/turbulenceProperties.template
index 9c5280d457..912cc2b71d 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/turbulenceProperties.template
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/constant/turbulenceProperties.template
@@ -14,7 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType  RAS;
+simulationType      RAS;
 
 RAS
 {
@@ -33,4 +33,5 @@ RAS
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/controlDict b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/controlDict
index 9e8967ad79..bba7e64f48 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/controlDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/controlDict
@@ -47,10 +47,10 @@ timePrecision   6;
 
 runTimeModifiable false;
 
-
 functions
 {
     #includeFunc "samples"
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/fvSchemes b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/fvSchemes
index f6e43c35c6..12faad3064 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/fvSchemes
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/fvSchemes
@@ -28,11 +28,14 @@ gradSchemes
 divSchemes
 {
     default             none;
-    div(phi,T)          bounded Gauss upwind;
     div(phi,U)          bounded Gauss upwind;
-    div(phi,k)          bounded Gauss upwind;
-    div(phi,epsilon)    bounded Gauss upwind;
-    div(phi,omega)      bounded Gauss upwind;
+
+    turbulence          bounded Gauss upwind;
+    div(phi,k)          $turbulence;
+    div(phi,epsilon)    $turbulence;
+    div(phi,omega)      $turbulence;
+
+    div(phi,T)          bounded Gauss upwind;
     div((nuEff*dev(T(grad(U)))))  Gauss linear;
     div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/fvSolution b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/fvSolution
index b94f3bbf0d..4afa515f47 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/fvSolution
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/fvSolution
@@ -69,7 +69,6 @@ relaxationFactors
 cache
 {
     grad(U);
-    grad(T);
 }
 
 
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/samples b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/samples
index b872c40b86..b7310e1760 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/samples
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/system/samples
@@ -1,4 +1,4 @@
-/*--------------------------------*- C++ -*----------------------------------*/
+// -*- C++ -*-
 
 type                sets;
 libs                (sampling);
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/T b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/T
index 76d7a7dc9d..8a09df123a 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/T
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/T
@@ -49,7 +49,7 @@ boundaryField
         type            slip;
     }
 
-    "upperInterface|lowerInterface"
+    "(upperInterface|lowerInterface)"
     {
         type            cyclicAMI;
         value           $internalField;
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/U b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/U
index e586bf5ba1..79f4c5cf15 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/U
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/U
@@ -45,7 +45,7 @@ boundaryField
         type            slip;
     }
 
-    "upperInterface|lowerInterface"
+    "(upperInterface|lowerInterface)"
     {
         type            cyclicAMI;
         value           $internalField;
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/alphat b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/alphat
index 4ce6c66d33..f846a50c5e 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/alphat
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/alphat
@@ -50,7 +50,7 @@ boundaryField
         type            slip;
     }
 
-    "upperInterface|lowerInterface"
+    "(upperInterface|lowerInterface)"
     {
         type            cyclicAMI;
         value           $internalField;
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/epsilon b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/epsilon
index 25060dad77..31ccb4127c 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/epsilon
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/epsilon
@@ -48,7 +48,7 @@ boundaryField
         type            slip;
     }
 
-    "upperInterface|lowerInterface"
+    "(upperInterface|lowerInterface)"
     {
         type            cyclicAMI;
         value           $internalField;
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/k b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/k
index d6ac7f4e8a..9b8037e8cd 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/k
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/k
@@ -45,7 +45,7 @@ boundaryField
         type            slip;
     }
 
-    "upperInterface|lowerInterface"
+    "(upperInterface|lowerInterface)"
     {
         type            cyclicAMI;
         value           $internalField;
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/nut b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/nut
index 1e2131aeeb..4ef7575c0c 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/nut
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/nut
@@ -48,7 +48,7 @@ boundaryField
         type            slip;
     }
 
-    "upperInterface|lowerInterface"
+    "(upperInterface|lowerInterface)"
     {
         type            cyclicAMI;
         value           $internalField;
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/omega b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/omega
index d6ef661764..173e484f60 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/omega
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/omega
@@ -40,7 +40,7 @@ boundaryField
         type            slip;
     }
 
-    "upperInterface|lowerInterface"
+    "(upperInterface|lowerInterface)"
     {
         type            cyclicAMI;
         value           $internalField;
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/p_rgh b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/p_rgh
index ebb3cb2f2a..4400c8dfc1 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/p_rgh
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/0.orig/p_rgh
@@ -38,7 +38,7 @@ boundaryField
         value           uniform 10.123;
     }
 
-    "upperInterface|lowerInterface"
+    "(upperInterface|lowerInterface)"
     {
         type            cyclicAMI;
         value           $internalField;
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/constant/fvOptions b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/constant/fvOptions
index 5e5197a3f5..c4193a5d55 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/constant/fvOptions
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/constant/fvOptions
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/constant/g b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/constant/g
index 07cc25626c..5d28f81891 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/constant/g
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/constant/g
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       uniformDimensionedVectorField;
-    location    "constant";
     object      g;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/controlDict b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/controlDict
index 8499961fbb..206e3f2844 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/controlDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/controlDict
@@ -47,7 +47,6 @@ timePrecision   6;
 
 runTimeModifiable false;
 
-
 functions
 {
     fieldAverage1
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/fvSchemes b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/fvSchemes
index ebd32d1712..e4c9deaf5a 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/fvSchemes
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/fvSchemes
@@ -28,11 +28,14 @@ gradSchemes
 divSchemes
 {
     default                       none;
-    div(phi,T)                    bounded Gauss upwind;
     div(phi,U)                    bounded Gauss upwind;
-    div(phi,k)                    bounded Gauss upwind;
-    div(phi,epsilon)              bounded Gauss upwind;
-    div(phi,omega)                bounded Gauss upwind;
+
+    turbulence                    bounded Gauss upwind;
+    div(phi,k)                    $turbulence;
+    div(phi,epsilon)              $turbulence;
+    div(phi,omega)                $turbulence;
+
+    div(phi,T)                    bounded Gauss upwind;
     div((nuEff*dev(T(grad(U)))))  Gauss linear;
     div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/fvSolution b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/fvSolution
index 1ecd83519e..7de8b41635 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/fvSolution
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/fvSolution
@@ -70,7 +70,6 @@ relaxationFactors
 cache
 {
     grad(U);
-    grad(T);
 }
 
 
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/mapFieldsDict b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/mapFieldsDict
index 1ba407e72c..dbc7e46af0 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/mapFieldsDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/mapFieldsDict
@@ -21,4 +21,5 @@ cuttingPatches
     terrain top
 );
 
+
 // ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/samples b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/samples
index d4af390008..13ec1014e9 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/samples
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/system/samples
@@ -1,4 +1,4 @@
-/*--------------------------------*- C++ -*----------------------------------*/
+// -*- C++ -*-
 
 type                sets;
 libs                (sampling);
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/T b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/T
index b49b4daa7a..c292e67f27 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/T
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/T
@@ -36,7 +36,7 @@ boundaryField
         value          uniform 300;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/U b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/U
index e53a393ada..be1a3d7d72 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/U
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/U
@@ -31,7 +31,7 @@ boundaryField
         type            slip;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/alphat b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/alphat
index 089c249d88..61630db1e3 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/alphat
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/alphat
@@ -37,7 +37,7 @@ boundaryField
         value           uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/epsilon b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/epsilon
index 9286e72f9c..29a6bb537c 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/epsilon
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/epsilon
@@ -34,7 +34,7 @@ boundaryField
         type            slip;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/k b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/k
index 30b68562af..92b83b2515 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/k
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/k
@@ -31,7 +31,7 @@ boundaryField
         type            slip;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/leafAreaDensity b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/leafAreaDensity
index 50ac7e229e..9a1e326947 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/leafAreaDensity
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/leafAreaDensity
@@ -32,7 +32,7 @@ boundaryField
         value           uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/nut b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/nut
index eddbcff03b..75eef89bf5 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/nut
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/nut
@@ -36,7 +36,7 @@ boundaryField
         value           uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/omega b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/omega
index 060bff4ef6..be97e78321 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/omega
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/omega
@@ -34,7 +34,7 @@ boundaryField
         type            slip;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/p_rgh b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/p_rgh
index 931bf6fde2..b723813d7e 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/p_rgh
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/p_rgh
@@ -20,18 +20,17 @@ internalField   uniform 0.0;
 
 boundaryField
 {
-    "bottom|top"
+    "(bottom|top)"
     {
        type             fixedFluxPressure;
        rho              rhok;
        value            uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
-
 }
 
 
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/plantCd b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/plantCd
index 5fe266d6e7..25c4871a43 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/plantCd
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/plantCd
@@ -20,13 +20,13 @@ internalField  uniform  0.2;
 
 boundaryField
 {
-    "bottom|top"
+    "(bottom|top)"
     {
         type            fixedValue;
         value           uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/qPlant.template b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/qPlant.template
index 7a57452fa9..f8ebebffe9 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/qPlant.template
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/0.orig/qPlant.template
@@ -20,13 +20,13 @@ internalField  uniform  Q_PLANT;
 
 boundaryField
 {
-    "bottom|top"
+    "(bottom|top)"
     {
         type            fixedValue;
         value           uniform 0;
     }
 
-    "inlet|outlet|left|right"
+    "(inlet|outlet|left|right)"
     {
         type            cyclic;
     }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/Allrun b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/Allrun-parallel
similarity index 100%
rename from tutorials/verificationAndValidation/atmosphericModels/atmForestStability/Allrun
rename to tutorials/verificationAndValidation/atmosphericModels/atmForestStability/Allrun-parallel
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/fvOptions.template b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/fvOptions.template
index 960ad8edff..58307189e3 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/fvOptions.template
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/fvOptions.template
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
     object      fvOptions;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -26,6 +25,7 @@ pressureGradient
     }
 }
 
+
 atmCoriolisUSource1
 {
     type                atmCoriolisUSource;
@@ -36,6 +36,7 @@ atmCoriolisUSource1
     }
 }
 
+
 atmAmbientTurbSource1
 {
     type                atmAmbientTurbSource;
@@ -47,6 +48,7 @@ atmAmbientTurbSource1
     }
 }
 
+
 atmBuoyancyTurbSource1
 {
     type                atmBuoyancyTurbSource;
@@ -59,6 +61,7 @@ atmBuoyancyTurbSource1
     }
 }
 
+
 atmLengthScaleTurbSource1
 {
     type                atmLengthScaleTurbSource;
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/g b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/g
index 07cc25626c..5d28f81891 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/g
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/g
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       uniformDimensionedVectorField;
-    location    "constant";
     object      g;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/turbulenceProperties.template b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/turbulenceProperties.template
index 09705fd041..912cc2b71d 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/turbulenceProperties.template
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/constant/turbulenceProperties.template
@@ -14,7 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType  RAS;
+simulationType      RAS;
 
 RAS
 {
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/controlDict b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/controlDict
index 7e5b0b5b8e..400c2f0a2e 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/controlDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/controlDict
@@ -47,7 +47,6 @@ timePrecision   6;
 
 runTimeModifiable false;
 
-
 functions
 {
     ObukhovLength1
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/fvSchemes b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/fvSchemes
index 309091fa81..b8aa5c1731 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/fvSchemes
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/fvSchemes
@@ -28,11 +28,14 @@ gradSchemes
 divSchemes
 {
     default             none;
-    div(phi,T)          bounded Gauss upwind;
     div(phi,U)          bounded Gauss upwind;
-    div(phi,k)          bounded Gauss upwind;
-    div(phi,epsilon)    bounded Gauss upwind;
-    div(phi,omega)      bounded Gauss upwind;
+
+    turbulence          bounded Gauss upwind;
+    div(phi,k)          $turbulence;
+    div(phi,epsilon)    $turbulence;
+    div(phi,omega)      $turbulence;
+
+    div(phi,T)          bounded Gauss upwind;
     div((nuEff*dev(T(grad(U)))))  Gauss linear;
     div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/fvSolution b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/fvSolution
index b94f3bbf0d..4afa515f47 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/fvSolution
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/fvSolution
@@ -69,7 +69,6 @@ relaxationFactors
 cache
 {
     grad(U);
-    grad(T);
 }
 
 
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/samples b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/samples
index d9a88977bc..ddbdb8ed88 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/samples
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/samples
@@ -1,4 +1,4 @@
-/*--------------------------------*- C++ -*----------------------------------*/
+// -*- C++ -*-
 
 type                sets;
 libs                (sampling);
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/turbulenceFields b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/turbulenceFields
index f76ff843cf..b5444da37d 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/turbulenceFields
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/system/turbulenceFields
@@ -1,4 +1,4 @@
-/*--------------------------------*- C++ -*----------------------------------*/
+// -*- C++ -*-
 
 type            turbulenceFields;
 libs            (fieldFunctionObjects);