diff --git a/src/atmosphericModels/Make/files b/src/atmosphericModels/Make/files
index 644149f5553786301485f97d36296cc4c66358b3..887e1e60596c2febb73f3cce83da8d626289c7fd 100644
--- a/src/atmosphericModels/Make/files
+++ b/src/atmosphericModels/Make/files
@@ -1,10 +1,22 @@
+/* Models */
+
+atmosphericTurbulentTransportModels.C
+porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C
+
+
+/* Boundary conditions */
+
 derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C
 derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
 derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C
 derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
+derivedFvPatchFields/atmBoundaryLayerInletOmega/atmBoundaryLayerInletOmegaFvPatchScalarField.C
 derivedFvPatchFields/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.C
 
-atmosphericTurbulentTransportModels.C
-porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C
+
+/* Wall function BCs */
+
+
+/* fvOptions */
 
 LIB = $(FOAM_LIBBIN)/libatmosphericModels
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C
index d80ccaac8dbc6abf5ce30ae4361b27f976b7aa0a..dca3c7b8d046b54f006818b1239a47e8400c74a6 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C
@@ -37,16 +37,20 @@ namespace Foam
 
 atmBoundaryLayer::atmBoundaryLayer(const Time& time, const polyPatch& pp)
 :
+    initABL_(false),
+    kappa_(0.41),
+    Cmu_(0.09),
+    C1_(0.0),
+    C2_(1.0),
+    ppMin_((boundBox(pp.points())).min()),
     time_(time),
     patch_(pp),
     flowDir_(time, "flowDir"),
     zDir_(time, "zDir"),
-    kappa_(0.41),
-    Cmu_(0.09),
     Uref_(time, "Uref"),
     Zref_(time, "Zref"),
     z0_(),
-    zGround_()
+    d_()
 {}
 
 
@@ -57,6 +61,15 @@ atmBoundaryLayer::atmBoundaryLayer
     const dictionary& dict
 )
 :
+    initABL_(dict.getOrDefault<Switch>("initABL", true)),
+    kappa_
+    (
+        dict.getCheckOrDefault<scalar>("kappa", 0.41, scalarMinMax::ge(SMALL))
+    ),
+    Cmu_(dict.getCheckOrDefault<scalar>("Cmu", 0.09, scalarMinMax::ge(SMALL))),
+    C1_(dict.getOrDefault("C1", 0.0)),
+    C2_(dict.getOrDefault("C2", 1.0)),
+    ppMin_((boundBox(pp.points())).min()),
     time_(time),
     patch_(pp),
     flowDir_(TimeFunction1<vector>(time, "flowDir", dict)),
@@ -66,7 +79,7 @@ atmBoundaryLayer::atmBoundaryLayer
     Uref_(TimeFunction1<scalar>(time, "Uref", dict)),
     Zref_(TimeFunction1<scalar>(time, "Zref", dict)),
     z0_(PatchFunction1<scalar>::New(pp, "z0", dict)),
-    zGround_(PatchFunction1<scalar>::New(pp, "zGround", dict))
+    d_(PatchFunction1<scalar>::New(pp, "d", dict))
 {}
 
 
@@ -77,31 +90,39 @@ atmBoundaryLayer::atmBoundaryLayer
     const fvPatchFieldMapper& mapper
 )
 :
+    initABL_(abl.initABL_),
+    kappa_(abl.kappa_),
+    Cmu_(abl.Cmu_),
+    C1_(abl.C1_),
+    C2_(abl.C2_),
+    ppMin_(abl.ppMin_),
     time_(abl.time_),
     patch_(patch.patch()),
     flowDir_(abl.flowDir_),
     zDir_(abl.zDir_),
-    kappa_(abl.kappa_),
-    Cmu_(abl.Cmu_),
     Uref_(abl.Uref_),
     Zref_(abl.Zref_),
     z0_(abl.z0_.clone(patch_)),
-    zGround_(abl.zGround_.clone(patch_))
+    d_(abl.d_.clone(patch_))
 {}
 
 
 atmBoundaryLayer::atmBoundaryLayer(const atmBoundaryLayer& abl)
 :
+    initABL_(abl.initABL_),
+    kappa_(abl.kappa_),
+    Cmu_(abl.Cmu_),
+    C1_(abl.C1_),
+    C2_(abl.C2_),
+    ppMin_(abl.ppMin_),
     time_(abl.time_),
     patch_(abl.patch_),
     flowDir_(abl.flowDir_),
     zDir_(abl.zDir_),
-    kappa_(abl.kappa_),
-    Cmu_(abl.Cmu_),
     Uref_(abl.Uref_),
     Zref_(abl.Zref_),
     z0_(abl.z0_.clone(patch_)),
-    zGround_(abl.zGround_.clone(patch_))
+    d_(abl.d_.clone(patch_))
 {}
 
 
@@ -110,13 +131,13 @@ atmBoundaryLayer::atmBoundaryLayer(const atmBoundaryLayer& abl)
 vector atmBoundaryLayer::flowDir() const
 {
     const scalar t = time_.timeOutputValue();
-    vector dir(flowDir_.value(t));
+    const vector dir(flowDir_.value(t));
     const scalar magDir = mag(dir);
 
     if (magDir < SMALL)
     {
         FatalErrorInFunction
-            << "magnitude of " << flowDir_.name()
+            << "magnitude of " << flowDir_.name() << " = " << magDir
             << " vector must be greater than zero"
             << abort(FatalError);
     }
@@ -128,13 +149,13 @@ vector atmBoundaryLayer::flowDir() const
 vector atmBoundaryLayer::zDir() const
 {
     const scalar t = time_.timeOutputValue();
-    vector dir(zDir_.value(t));
+    const vector dir(zDir_.value(t));
     const scalar magDir = mag(dir);
 
     if (magDir < SMALL)
     {
         FatalErrorInFunction
-            << "magnitude of " << zDir_.name()
+            << "magnitude of " << zDir_.name() << " = " << magDir
             << " vector must be greater than zero"
             << abort(FatalError);
     }
@@ -149,14 +170,22 @@ tmp<scalarField> atmBoundaryLayer::Ustar(const scalarField& z0) const
     const scalar Uref = Uref_.value(t);
     const scalar Zref = Zref_.value(t);
 
-    return kappa_*Uref/(log((Zref + z0)/z0));
+    if (Zref < 0)
+    {
+        FatalErrorInFunction
+            << "Negative entry in " << Zref_.name() << " = " << Zref
+            << abort(FatalError);
+    }
+
+    // (derived from RH:Eq. 6, HW:Eq. 5)
+    return kappa_*Uref/log((Zref + z0)/z0);
 }
 
 
 void atmBoundaryLayer::autoMap(const fvPatchFieldMapper& mapper)
 {
     z0_->autoMap(mapper);
-    zGround_->autoMap(mapper);
+    d_->autoMap(mapper);
 }
 
 
@@ -167,51 +196,80 @@ void atmBoundaryLayer::rmap
 )
 {
     z0_->rmap(abl.z0_(), addr);
-    zGround_->rmap(abl.zGround_(), addr);
+    d_->rmap(abl.d_(), addr);
 }
 
 
-tmp<vectorField> atmBoundaryLayer::U(const vectorField& p) const
+tmp<vectorField> atmBoundaryLayer::U(const vectorField& pCf) const
 {
     const scalar t = time_.timeOutputValue();
-    const scalarField zGround(zGround_->value(t));
+    const scalarField d(d_->value(t));
     const scalarField z0(max(z0_->value(t), ROOTVSMALL));
+    const scalar groundMin = zDir() & ppMin_;
 
-    scalarField Un((Ustar(z0)/kappa_)*log(((zDir() & p) - zGround + z0)/z0));
+    // (YGCJ:Table 1, RH:Eq. 6, HW:Eq. 5)
+    scalarField Un
+    (
+        (Ustar(z0)/kappa_)*log(((zDir() & pCf) - groundMin - d + z0)/z0)
+    );
 
     return flowDir()*Un;
 }
 
 
-tmp<scalarField> atmBoundaryLayer::k(const vectorField& p) const
+tmp<scalarField> atmBoundaryLayer::k(const vectorField& pCf) const
 {
     const scalar t = time_.timeOutputValue();
+    const scalarField d(d_->value(t));
     const scalarField z0(max(z0_->value(t), ROOTVSMALL));
+    const scalar groundMin = zDir() & ppMin_;
 
-    return sqr(Ustar(z0))/sqrt(Cmu_);
+    // (YGCJ:Eq. 21; RH:Eq. 7, HW:Eq. 6 when C1=0 and C2=1)
+    return
+        sqr(Ustar(z0))/sqrt(Cmu_)
+       *sqrt(C1_*log(((zDir() & pCf) - groundMin - d + z0)/z0) + C2_);
 }
 
 
-tmp<scalarField> atmBoundaryLayer::epsilon(const vectorField& p) const
+tmp<scalarField> atmBoundaryLayer::epsilon(const vectorField& pCf) const
 {
     const scalar t = time_.timeOutputValue();
-    const scalarField zGround(zGround_->value(t));
+    const scalarField d(d_->value(t));
     const scalarField z0(max(z0_->value(t), ROOTVSMALL));
+    const scalar groundMin = zDir() & ppMin_;
 
-    return pow3(Ustar(z0))/(kappa_*((zDir() & p) - zGround + z0));
+    // (YGCJ:Eq. 22; RH:Eq. 8, HW:Eq. 7 when C1=0 and C2=1)
+    return
+        pow3(Ustar(z0))/(kappa_*((zDir() & pCf) - groundMin - d + z0))
+       *sqrt(C1_*log(((zDir() & pCf) - groundMin - d + z0)/z0) + C2_);
+}
+
+
+tmp<scalarField> atmBoundaryLayer::omega(const vectorField& pCf) const
+{
+    const scalar t = time_.timeOutputValue();
+    const scalarField d(d_->value(t));
+    const scalarField z0(max(z0_->value(t), ROOTVSMALL));
+    const scalar groundMin = zDir() & ppMin_;
+
+    // (YGJ:Eq. 13)
+    return Ustar(z0)/(kappa_*sqrt(Cmu_)*((zDir() & pCf) - groundMin - d + z0));
 }
 
 
 void atmBoundaryLayer::write(Ostream& os) const
 {
-    z0_->writeData(os) ;
-    flowDir_.writeData(os);
-    zDir_.writeData(os);
+    os.writeEntry("initABL", initABL_);
     os.writeEntry("kappa", kappa_);
     os.writeEntry("Cmu", Cmu_);
+    os.writeEntry("C1", C1_);
+    os.writeEntry("C2", C2_);
+    flowDir_.writeData(os);
+    zDir_.writeData(os);
     Uref_.writeData(os);
     Zref_.writeData(os);
-    zGround_->writeData(os);
+    z0_->writeData(os) ;
+    d_->writeData(os);
 }
 
 
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H
index 888801476f99f22a61f31c93b5d596a6faa1e97f..c030615facd4a56492859e9ddda7b94c684e5f7e 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2014-2018 OpenFOAM Foundation
-    Copyright (C) 2015 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,84 +31,156 @@ Group
     grpRASBoundaryConditions grpInletBoundaryConditions
 
 Description
-    This class provides functions to evaluate the velocity and turbulence
-    distributions appropriate for atmospheric boundary layers (ABL).
+    Base class to set log-law type ground-normal inflow boundary conditions for
+    wind velocity and turbulence quantities for homogeneous, two-dimensional,
+    dry-air, equilibrium and neutral atmospheric boundary layer (ABL) modelling.
 
-    The profile is derived from the friction velocity, flow direction and
-    "vertical" direction:
+    The ground-normal profile expressions are due to \c YGCJ
+    (refer to references below) whereat \c RH expressions were generalised:
 
-        \f[
-            U = \frac{U^*}{\kappa} ln\left(\frac{z - z_g + z_0}{z_0}\right)
-        \f]
+    \f[
+        u = \frac{u^*}{\kappa} \ln \left( \frac{z - d + z_0}{z_0} \right)
+    \f]
 
-        \f[
-            k = \frac{(U^*)^2}{\sqrt{C_{\mu}}}
-        \f]
+    \f[
+        v = w = 0
+    \f]
 
-        \f[
-            \epsilon = \frac{(U^*)^3}{\kappa(z - z_g + z_0)}
-        \f]
+    \f[
+        k = \frac{(u^*)^2}{\sqrt{C_\mu}}
+        \sqrt{C_1 \ln \left( \frac{z - d + z_0}{z_0} \right) + C_2}
+    \f]
+
+    \f[
+        \epsilon = \frac{(u^*)^3}{\kappa (z - d + z_0)}
+        \sqrt{C_1 \ln \left( \frac{z - d + z_0}{z_0} \right) + C_2}
+    \f]
+
+    \f[
+        \omega = \frac{u^*}{\kappa \sqrt{C_\mu}} \frac{1}{z - d + z_0}
+    \f]
+
+    \f[
+        u^* =
+            \frac{u_{ref} \kappa}{\ln\left(\frac{z_{ref} + z_0}{z_0}\right)}
+    \f]
 
     where
     \vartable
-        U^*     | Friction velocity
-        \kappa  | von Karman's constant
-        C_{\mu} | Turbulence viscosity coefficient
-        z       | Vertical coordinate
-        z_0     | Surface roughness height [m]
-        z_g     | Minimum z-coordinate [m]
-    \endvartable
-    and
-        \f[
-            U^* = \kappa\frac{U_{ref}}{ln\left(\frac{Z_{ref} + z_0}{z_0}\right)}
-        \f]
-    where
-    \vartable
-        U_{ref} | Reference velocity at \f$Z_{ref}\f$ [m/s]
-        Z_{ref} | Reference height [m]
+      u        | Ground-normal streamwise flow speed profile          [m/s]
+      v        | Spanwise flow speed                                  [m/s]
+      w        | Ground-normal flow speed                             [m/s]
+      k        | Ground-normal turbulent kinetic energy (TKE) profile [m^2/s^2]
+      \epsilon | Ground-normal TKE dissipation rate profile           [m^2/s^3]
+      \omega   | Ground-normal specific dissipation rate profile      [m^2/s^3]
+      u^*      | Friction velocity                                    [m/s]
+      \kappa   | von Kármán constant                                  [-]
+      C_\mu    | Empirical model constant                             [-]
+      z        | Ground-normal coordinate component                   [m]
+      d        | Ground-normal displacement height                    [m]
+      z_0      | Aerodynamic roughness length                         [m]
+      u_{ref}  | Reference mean streamwise wind speed at \f$z_{ref}\f$ [m/s]
+      z_{ref}  | Reference height being used in \f$u^*\f$ estimations [m]
+      C_1      | Curve-fitting coefficient for \c YGCJ profiles       [-]
+      C_2      | Curve-fitting coefficient for \c YGCJ profiles       [-]
     \endvartable
 
-    Use in the atmBoundaryLayerInletVelocity, atmBoundaryLayerInletK and
-    atmBoundaryLayerInletEpsilon boundary conditions.
+Note
+    - The \c RH expressions are special cases of those in \c YGCJ when \c C1=0
+    and \c C2=1. Both \c C1 and \c C2 can be determined by nonlinear fitting
+    of (\c YGCJ:Eq. 19) with an experimental dataset for \c k. By default,
+    \c atmBoundaryLayerInlet boundary conditions compute \c RH expressions.
+    - \c z is the ground-normal height relative to the global minimum height
+    of the inlet patch; therefore, the minimum of \c z is always zero
+    irrespective of the absolute z-coordinate of the computational patch.
 
     Reference:
-        D.M. Hargreaves and N.G. Wright,  "On the use of the k-epsilon model
-        in commercial CFD software to model the neutral atmospheric boundary
-        layer", Journal of Wind Engineering and Industrial Aerodynamics
-        95(2007), pp 355-369.
+    \verbatim
+        The ground-normal profile expressions (tag:RH):
+            Richards, P. J., & Hoxey, R. P. (1993).
+            Appropriate boundary conditions for computational wind
+            engineering models using the k-ε turbulence model.
+            In Computational Wind Engineering 1 (pp. 145-153).
+            DOI:10.1016/B978-0-444-81688-7.50018-8
+
+        Modifications to preserve the profiles downstream (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
+
+        Expression generalisations to allow height
+        variation for turbulence quantities (tag:YGCJ):
+            Yang, Y., Gu, M., Chen, S., & Jin, X. (2009).
+            New inflow boundary conditions for modelling the neutral equilibrium
+            atmospheric boundary layer in computational wind engineering.
+            J. of Wind Engineering and Industrial Aerodynamics, 97(2), 88-95.
+            DOI:10.1016/j.jweia.2008.12.001
+
+        The generalised ground-normal profile expression for omega (tag:YGJ):
+            Yang, Y., Gu, M., & Jin, X., (2009).
+            New inflow boundary conditions for modelling the
+            neutral equilibrium atmospheric boundary layer in SST k-ω model.
+            In: The Seventh Asia-Pacific Conference on Wind Engineering,
+            November 8-12, Taipei, Taiwan.
+
+        Theoretical remarks (tag:E):
+            Emeis, S. (2013).
+            Wind Energy Meteorology: Atmospheric
+            Physics for Wind Power Generation.
+            Springer-Verlag Berlin Heidelberg.
+            DOI:10.1007/978-3-642-30523-8
+    \endverbatim
 
 Usage
     \table
-        Property     | Description                      | Required  | Default
-        flowDir      | Flow direction                   | yes       |
-        zDir         | Vertical direction               | yes       |
-        kappa        | von Karman's constant            | no        | 0.41
-        Cmu          | Turbulence viscosity coefficient | no        | 0.09
-        Uref         | Reference velocity [m/s]         | yes       |
-        Zref         | Reference height [m]             | yes       |
-        z0           | Surface roughness height [m]     | yes       |
-        zGround      | Minimum z-coordinate [m]         | yes       |
+      Property  | Description              | Type          | Req'd | Deflt
+      flowDir   | Flow direction           | TimeFunction1<vector> | yes | -
+      zDir      | Ground-normal direction  | TimeFunction1<vector> | yes | -
+      Uref      | Reference mean streamwise flow speed being used in <!--
+             --> \f$u^*\f$ estimations [m/s] | TimeFunction1<scalar> | yes | -
+      Zref      | Reference height being used in \f$u^*\f$ estimations [m]  <!--
+                                       --> | TimeFunction1<scalar> | yes | -
+      z0        | Surface roughness length [m] <!--
+                                       --> | PatchFunction1<scalar> | yes | -
+      d         | Displacement height [m] - see Notes <!--
+                                       --> | PatchFunction1<scalar> | yes | -
+      kappa     | von Kármán constant      | scalar | no  | 0.41
+      Cmu       | Empirical model constant | scalar | no  | 0.09
+      initABL   | Flag to initialise profiles with the theoretical <!--
+            --> ABL expressions, otherwise use "value" list        <!--
+                                       --> | bool   | no  | true
+      value     | ABL profile content when initABL=false <!--
+                                       --> | scalarList | conditional | -
+      phi       | Name of the flux field   | word         | no        | phi
+      C1        | Curve-fitting coefficient \c YGCJ profiles | scalar | no | 0.0
+      C2        | Curve-fitting coefficient \c YGCJ profiles | scalar | no | 1.0
     \endtable
 
-    Example of the boundary condition specification:
-    \verbatim
-    ground
-    {
-        type            atmBoundaryLayerInletVelocity;
-        flowDir         (1 0 0);
-        zDir            (0 0 1);
-        Uref            10.0;
-        Zref            20.0;
-        z0              uniform 0.1;
-        zGround         uniform 0.0;
-    }
-    \endverbatim
-
 Note
-    D.M. Hargreaves and N.G. Wright recommend Gamma epsilon in the
-    k-epsilon model should be changed from 1.3 to 1.11 for consistency.
-    The roughness height (Er) is given by Er = 20 z0 following the same
-    reference.
+    - The derived ABL expressions automatically satisfy the simplified transport
+    equation for \c k. Yet the same expressions only satisfy the simplified
+    transport equation for \c epsilon when the model constants \c sigmaEpsilon
+    is 1.11 with \c kappa=0.4 (\c HW:p. 358).
+    - \c atmBoundaryLayerInlet boundary conditions inherit \c inletOutlet
+    traits, so that a given inlet condition can be supplied from all sides of
+    the domain, e.g. a ground-normal cylinder domain having a single
+    inlet/outlet boundary where the changes between inlet and outlet depend
+    on the wind direction and patch normals, so that any change in inflow
+    orientation can be handled with the same mesh.
+    - \c d is the displacement height, and "is relevant for flows over forests
+    and cities" (E:p. 28). "The displacement height gives the vertical
+    displacement of the entire flow regime over areas which are densely covered
+    with obstacles such as trees or buildings" (E:p. 28).
+
+See also
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletVelocityFvPatchVectorField
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletKFvPatchScalarField
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletEpsilonFvPatchScalarField
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletOmegaFvPatchScalarField
+    - ExtendedCodeGuide::atmBoundaryLayer
 
 SourceFiles
     atmBoundaryLayer.C
@@ -133,7 +205,33 @@ namespace Foam
 
 class atmBoundaryLayer
 {
-    // Private data
+protected:
+
+    // Protected Data
+
+        //- Flag to initialise profiles with the theoretical ABL expressions,
+        //- otherwise initialises by using "value" entry content
+        Switch initABL_;
+
+
+private:
+
+    // Private Data
+
+        //- von Kármán constant
+        const scalar kappa_;
+
+        //- Empirical model constant
+        const scalar Cmu_;
+
+        //- Curve-fitting coefficient
+        const scalar C1_;
+
+        //- Curve-fitting coefficient
+        const scalar C2_;
+
+        //- Minimum coordinate vector of this patch
+        const vector ppMin_;
 
         //- Reference to the time database
         const Time& time_;
@@ -141,29 +239,23 @@ class atmBoundaryLayer
         //- Reference to the patch
         const polyPatch& patch_;
 
-        //- Flow direction
+        //- Streamwise flow direction
         TimeFunction1<vector> flowDir_;
 
-        //- Direction of the z-coordinate
+        //- Direction of the ground-normal coordinate
         TimeFunction1<vector> zDir_;
 
-        //- Von Karman constant
-        const scalar kappa_;
-
-        //- Turbulent viscosity coefficient
-        const scalar Cmu_;
-
-        //- Reference velocity
+        //- Reference mean streamwise flow speed being used in Ustar estimations
         TimeFunction1<scalar> Uref_;
 
-        //- Reference height
+        //- Reference height being used in Ustar estimations
         TimeFunction1<scalar> Zref_;
 
-        //- Surface roughness height
+        //- Surface roughness length
         autoPtr<PatchFunction1<scalar>> z0_;
 
-        //- Minimum coordinate value in z direction
-        autoPtr<PatchFunction1<scalar>> zGround_;
+        //- Displacement height
+        autoPtr<PatchFunction1<scalar>> d_;
 
 
 public:
@@ -193,14 +285,14 @@ public:
         atmBoundaryLayer(const atmBoundaryLayer&);
 
 
-    // Member functions
+    // Member Functions
 
         // Access
 
             //- Return flow direction
             vector flowDir() const;
 
-            //- Return z-direction
+            //- Return the ground-normal direction
             vector zDir() const;
 
             //- Return friction velocity
@@ -219,14 +311,16 @@ public:
         // Evaluate functions
 
             //- Return the velocity distribution for the ATM
-            tmp<vectorField> U(const vectorField& p) const;
+            tmp<vectorField> U(const vectorField& pCf) const;
 
             //- Return the turbulent kinetic energy distribution for the ATM
-            tmp<scalarField> k(const vectorField& p) const;
+            tmp<scalarField> k(const vectorField& pCf) const;
 
             //- Return the turbulent dissipation rate distribution for the ATM
-            tmp<scalarField> epsilon(const vectorField& p) const;
+            tmp<scalarField> epsilon(const vectorField& pCf) const;
 
+            //- 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.C b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
index 1a78cf0601bf60f2b2580ef7969ab70787ed9b07..aff00e89bae2604df1acf9abe4697951a267b9fb 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
@@ -28,6 +28,7 @@ License
 
 #include "atmBoundaryLayerInletEpsilonFvPatchScalarField.H"
 #include "addToRunTimeSelectionTable.H"
+#include "turbulenceModel.H"
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
 #include "surfaceFields.H"
@@ -68,13 +69,14 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
     refGrad() = 0;
     valueFraction() = 1;
 
-    if (dict.found("value"))
+    if (!initABL_)
     {
         scalarField::operator=(scalarField("value", dict, p.size()));
     }
     else
     {
         scalarField::operator=(refValue());
+        initABL_ = false;
     }
 }
 
@@ -148,8 +150,8 @@ void atmBoundaryLayerInletEpsilonFvPatchScalarField::rmap
 void atmBoundaryLayerInletEpsilonFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
-    atmBoundaryLayer::write(os);
     os.writeEntryIfDifferent<word>("phi", "phi", phiName_);
+    atmBoundaryLayer::write(os);
     writeEntry("value", os);
 }
 
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
index 30dea1d0704b4d8e7b7986c135bc075d2d245505..e6c8318612043582734d273ab0e76fab3d39586b 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,28 +31,74 @@ Group
     grpRASBoundaryConditions grpInletBoundaryConditions
 
 Description
-    This boundary condition specifies an inlet value for the turbulence
-    dissipation, \f$\epsilon\f$, appropriate for atmospheric boundary layers.
-
-    See Foam::atmBoundaryLayer for details.
+    This boundary condition provides a log-law type ground-normal inflow
+    boundary condition for 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:
+
+        \f[
+            \epsilon = \frac{(u^*)^3}{\kappa (z - d + z_0)}
+            \sqrt{C_1 \ln \left( \frac{z - d + z_0}{z_0} \right) + C_2}
+        \f]
+
+    where
+    \vartable
+      \epsilon | Ground-normal TKE dissipation rate profile     [m^2/s^3]
+      u^*      | Friction velocity                              [m/s]
+      \kappa   | von Kármán constant                            [-]
+      z        | Ground-normal coordinate component             [m]
+      z_0      | Aerodynamic roughness length                   [m]
+      d        | Ground-normal displacement height              [m]
+      C_1      | Curve-fitting coefficient for \c YGCJ profiles [-]
+      C_2      | Curve-fitting coefficient for \c YGCJ profiles [-]
+ \endvartable
+
+Usage
+    \heading Required fields
+    \plaintable
+        epsilon   | Turbulent kinetic energy dissipation rate   [m^2/s^3]
+    \endplaintable
 
     Example of the boundary condition specification:
     \verbatim
-    ground
+    inlet
     {
+        // Mandatory entries (unmodifiable)
         type            atmBoundaryLayerInletEpsilon;
+
+        // Mandatory (inherited) entries (unmodifiable)
+        flowDir         (1 0 0);
         z               (0 0 1);
         Uref            10.0;
-        Zref            20.0;
+        Zref            0.0;
         z0              uniform 0.1;
-        zGround         uniform 0.0;
+        d               uniform 0.0;
+
+        // Optional (inherited) entries (unmodifiable)
+        kappa           0.41;
+        Cmu             0.09;
+        initABL         true;
+        phi             phi;
+        C1              0.0;
+        C2              1.0;
+
+        // Conditional mandatory (inherited) entries (runtime modifiable)
+        value           uniform 0;    // when initABL=false
     }
     \endverbatim
 
+    The inherited entries are elaborated in:
+     - \link atmBoundaryLayer.H \endlink
+     - \link inletOutletFvPatchField.H \endlink
+
 See also
-    Foam::atmBoundaryLayer,
-    Foam::atmBoundaryLayerInletVelocityFvPatchVectorField,
-    Foam::atmBoundaryLayerInletKFvPatchScalarField
+    - Foam::atmBoundaryLayer
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletKFvPatchScalarField
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletOmegaFvPatchScalarField
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletVelocityFvPatchVectorField
+    - ExtendedCodeGuide::atmBoundaryLayer::atmBoundaryLayerInletEpsilon
 
 SourceFiles
     atmBoundaryLayerInletEpsilonFvPatchScalarField.C
@@ -104,7 +151,7 @@ public:
         );
 
         //- Construct by mapping given
-        //  atmBoundaryLayerInletEpsilonFvPatchScalarField onto a new patch
+        //- atmBoundaryLayerInletEpsilonFvPatchScalarField onto a new patch
         atmBoundaryLayerInletEpsilonFvPatchScalarField
         (
             const atmBoundaryLayerInletEpsilonFvPatchScalarField&,
@@ -142,12 +189,11 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
 
-
         // Mapping functions
 
             //- Map (and resize as needed) from self given a mapping object
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C
index b297eadfeae977a58b1f7e3fd02f5481b76c0527..7fc3d768e3f0d38032aa3925db54a86c7327d149 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C
@@ -68,13 +68,14 @@ atmBoundaryLayerInletKFvPatchScalarField
     refGrad() = 0;
     valueFraction() = 1;
 
-    if (dict.found("value"))
+    if (!initABL_)
     {
         scalarField::operator=(scalarField("value", dict, p.size()));
     }
     else
     {
         scalarField::operator=(refValue());
+        initABL_ = false;
     }
 }
 
@@ -148,8 +149,8 @@ void atmBoundaryLayerInletKFvPatchScalarField::rmap
 void atmBoundaryLayerInletKFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
-    atmBoundaryLayer::write(os);
     os.writeEntryIfDifferent<word>("phi", "phi", phiName_);
+    atmBoundaryLayer::write(os);
     writeEntry("value", os);
 }
 
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.H
index f4f98cce1fced409089a806d4c6f1feffdc80b22..a5a2e47f086860b6d280f265959700db1de4b7f6 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2014-2018 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,28 +31,71 @@ Group
     grpRASBoundaryConditions grpInletBoundaryConditions
 
 Description
-    This boundary condition specifies an inlet value for the turbulence
-    kinetic energy, \f$k\f$, appropriate for atmospheric boundary layers.
-
-    See Foam::atmBoundaryLayer for details.
+    This boundary condition provides a log-law type ground-normal inflow
+    boundary condition for turbulent kinetic energy, i.e. \c k,
+    for homogeneous, two-dimensional, dry-air, equilibrium and neutral
+    atmospheric boundary layer modelling.
+
+    The ground-normal \c k profile expression:
+
+        \f[
+            k = \frac{(u^*)^2}{\sqrt{C_\mu}}
+            \sqrt{C_1 \ln \left( \frac{z - d + z_0}{z_0} \right) + C_2}
+        \f]
+
+    where
+    \vartable
+      k        | Ground-normal turbulent kinetic energy profile [m^2/s^3]
+      u^*      | Friction velocity                              [m/s]
+      C_\mu    | Empirical model constant                       [-]
+      C_1      | Curve-fitting coefficient for \c YGCJ profiles [-]
+      C_2      | Curve-fitting coefficient for \c YGCJ profiles [-]
+    \endvartable
+
+Usage
+    \heading Required fields
+    \plaintable
+        k    | Turbulent kinetic energy                         [m^2/s^2]
+    \endplaintable
 
     Example of the boundary condition specification:
     \verbatim
-    ground
+    inlet
     {
+        // Mandatory entries (unmodifiable)
         type            atmBoundaryLayerInletK;
-        z               (0 0 1);
+
+        // Mandatory (inherited) entries (unmodifiable)
         Uref            10.0;
-        Zref            20.0;
+        Zref            0.0;
         z0              uniform 0.1;
-        zGround         uniform 0.0;
+        d               uniform 0.0;
+        flowDir         (1 0 0);    // not used
+        z               (0 0 1);    // not used
+
+        // Optional (inherited) entries (unmodifiable)
+        kappa           0.41;
+        Cmu             0.09;
+        initABL         true;
+        phi             phi;
+        C1              0.0;
+        C2              1.0;
+
+        // Conditional mandatory (inherited) entries (runtime modifiable)
+        value           uniform 0;    // when initABL=false
     }
     \endverbatim
 
+    The inherited entries are elaborated in:
+     - \link atmBoundaryLayer.H \endlink
+     - \link inletOutletFvPatchField.H \endlink
+
 See also
-    Foam::atmBoundaryLayer,
-    Foam::atmBoundaryLayerInletVelocityFvPatchVectorField,
-    Foam::atmBoundaryLayerInletEpsilonFvPatchScalarField
+    - Foam::atmBoundaryLayer
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletVelocityFvPatchVectorField
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletEpsilonFvPatchScalarField
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletOmegaFvPatchScalarField
+    - ExtendedCodeGuide::atmBoundaryLayer::atmBoundaryLayerInletK
 
 SourceFiles
     atmBoundaryLayerInletKFvPatchScalarField.C
@@ -104,7 +148,7 @@ public:
         );
 
         //- Construct by mapping given
-        //  atmBoundaryLayerInletKFvPatchScalarField onto a new patch
+        //- atmBoundaryLayerInletKFvPatchScalarField onto a new patch
         atmBoundaryLayerInletKFvPatchScalarField
         (
             const atmBoundaryLayerInletKFvPatchScalarField&,
@@ -142,7 +186,7 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletOmega/atmBoundaryLayerInletOmegaFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletOmega/atmBoundaryLayerInletOmegaFvPatchScalarField.C
new file mode 100644
index 0000000000000000000000000000000000000000..7cc4e9490d762f9a28ae5157b52082b2f33fd7e5
--- /dev/null
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletOmega/atmBoundaryLayerInletOmegaFvPatchScalarField.C
@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "atmBoundaryLayerInletOmegaFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+atmBoundaryLayerInletOmegaFvPatchScalarField::
+atmBoundaryLayerInletOmegaFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    inletOutletFvPatchScalarField(p, iF),
+    atmBoundaryLayer(iF.time(), p.patch())
+{}
+
+
+atmBoundaryLayerInletOmegaFvPatchScalarField::
+atmBoundaryLayerInletOmegaFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    inletOutletFvPatchScalarField(p, iF),
+    atmBoundaryLayer(iF.time(), p.patch(), dict)
+{
+    phiName_ = dict.getOrDefault<word>("phi", "phi");
+
+    refValue() = omega(patch().Cf());
+    refGrad() = 0;
+    valueFraction() = 1;
+
+    if (!initABL_)
+    {
+        scalarField::operator=(scalarField("value", dict, p.size()));
+    }
+    else
+    {
+        scalarField::operator=(refValue());
+        initABL_ = false;
+    }
+}
+
+
+atmBoundaryLayerInletOmegaFvPatchScalarField::
+atmBoundaryLayerInletOmegaFvPatchScalarField
+(
+    const atmBoundaryLayerInletOmegaFvPatchScalarField& psf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    inletOutletFvPatchScalarField(psf, p, iF, mapper),
+    atmBoundaryLayer(psf, p, mapper)
+{}
+
+
+atmBoundaryLayerInletOmegaFvPatchScalarField::
+atmBoundaryLayerInletOmegaFvPatchScalarField
+(
+    const atmBoundaryLayerInletOmegaFvPatchScalarField& psf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    inletOutletFvPatchScalarField(psf, iF),
+    atmBoundaryLayer(psf)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void atmBoundaryLayerInletOmegaFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    refValue() = omega(patch().Cf());
+
+    inletOutletFvPatchScalarField::updateCoeffs();
+}
+
+
+void atmBoundaryLayerInletOmegaFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    inletOutletFvPatchScalarField::autoMap(m);
+    atmBoundaryLayer::autoMap(m);
+}
+
+
+void atmBoundaryLayerInletOmegaFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& psf,
+    const labelList& addr
+)
+{
+    inletOutletFvPatchScalarField::rmap(psf, addr);
+
+    const atmBoundaryLayerInletOmegaFvPatchScalarField& blpsf =
+        refCast<const atmBoundaryLayerInletOmegaFvPatchScalarField>(psf);
+
+    atmBoundaryLayer::rmap(blpsf, addr);
+}
+
+
+void atmBoundaryLayerInletOmegaFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeEntryIfDifferent<word>("phi", "phi", phiName_);
+    atmBoundaryLayer::write(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    atmBoundaryLayerInletOmegaFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletOmega/atmBoundaryLayerInletOmegaFvPatchScalarField.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletOmega/atmBoundaryLayerInletOmegaFvPatchScalarField.H
new file mode 100644
index 0000000000000000000000000000000000000000..5b755ae223ddcc6a17c07867bc01b7ee29a31681
--- /dev/null
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletOmega/atmBoundaryLayerInletOmegaFvPatchScalarField.H
@@ -0,0 +1,223 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::atmBoundaryLayerInletOmegaFvPatchScalarField
+
+Group
+    grpRASBoundaryConditions grpInletBoundaryConditions
+
+Description
+    This boundary condition provides a log-law type ground-normal inflow
+    boundary condition for specific dissipation rate, i.e. \c omega,
+    for homogeneous, two-dimensional, dry-air, equilibrium and neutral
+    atmospheric boundary layer modelling.
+
+    The ground-normal \c omega profile expression:
+
+        \f[
+            \omega = \frac{u^*}{\kappa \sqrt{C_\mu}} \frac{1}{z - d + z_0}
+        \f]
+
+    where
+    \vartable
+      \omega | Ground-normal specific dissipation rate profile [m^2/s^3]
+      u^*    | Friction velocity                               [m/s]
+      \kappa | von Kármán constant                             [-]
+      C_\mu  | Empirical model constant                        [-]
+      z      | Ground-normal coordinate component              [m]
+      d      | Ground-normal displacement height               [m]
+      z_0    | Aerodynamic roughness length                    [m]
+    \endvartable
+
+Usage
+    \heading Required fields
+    \plaintable
+        omega   | Specific dissipation rate                    [1/s]
+    \endplaintable
+
+    Example of the boundary condition specification:
+    \verbatim
+    inlet
+    {
+        // Mandatory entries (unmodifiable)
+        type            atmBoundaryLayerInletOmega;
+
+        // Mandatory (inherited) entries (unmodifiable)
+        flowDir         (1 0 0);
+        z               (0 0 1);
+        Uref            10.0;
+        Zref            0.0;
+        z0              uniform 0.1;
+        d               uniform 0.0;
+
+        // Optional (inherited) entries (unmodifiable)
+        kappa           0.41;
+        Cmu             0.09;
+        initABL         true;
+        phi             phi;
+        C1              0.0;
+        C2              1.0;
+
+        // Conditional mandatory (inherited) entries (runtime modifiable)
+        value           uniform 0;    // when initABL=false
+    }
+    \endverbatim
+
+    The inherited entries are elaborated in:
+     - \link atmBoundaryLayer.H \endlink
+     - \link inletOutletFvPatchField.H \endlink
+
+See also
+    - Foam::atmBoundaryLayer
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletVelocityFvPatchVelocityField
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletEpsilonFvPatchScalarField
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletKFvPatchScalarField
+    - ExtendedCodeGuide::atmBoundaryLayer::atmBoundaryLayerInletOmega
+
+SourceFiles
+    atmBoundaryLayerInletOmegaFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef atmBoundaryLayerInletOmegaFvPatchScalarField_H
+#define atmBoundaryLayerInletOmegaFvPatchScalarField_H
+
+#include "fvPatchFields.H"
+#include "inletOutletFvPatchFields.H"
+#include "atmBoundaryLayer.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+        Class atmBoundaryLayerInletOmegaFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class atmBoundaryLayerInletOmegaFvPatchScalarField
+:
+    public inletOutletFvPatchScalarField,
+    public atmBoundaryLayer
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("atmBoundaryLayerInletOmega");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        atmBoundaryLayerInletOmegaFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        atmBoundaryLayerInletOmegaFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //- atmBoundaryLayerInletOmegaFvPatchScalarField onto a new patch
+        atmBoundaryLayerInletOmegaFvPatchScalarField
+        (
+            const atmBoundaryLayerInletOmegaFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new atmBoundaryLayerInletOmegaFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        atmBoundaryLayerInletOmegaFvPatchScalarField
+        (
+            const atmBoundaryLayerInletOmegaFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new atmBoundaryLayerInletOmegaFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member Functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
index 23e2e21ac8989e4c16149e3a8a7f99b52a617583..077a47391b1d05112207ad2a8796d484c99e013f 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
@@ -68,13 +68,14 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
     refGrad() = Zero;
     valueFraction() = 1;
 
-    if (dict.found("value"))
+    if (!initABL_)
     {
         vectorField::operator=(vectorField("value", dict, p.size()));
     }
     else
     {
         vectorField::operator=(refValue());
+        initABL_ = false;
     }
 }
 
@@ -148,8 +149,8 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::rmap
 void atmBoundaryLayerInletVelocityFvPatchVectorField::write(Ostream& os) const
 {
     fvPatchVectorField::write(os);
-    atmBoundaryLayer::write(os);
     os.writeEntryIfDifferent<word>("phi", "phi", phiName_);
+    atmBoundaryLayer::write(os);
     writeEntry("value", os);
 }
 
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
index 4a5a68ffaf184dfb759260cc17c211ecd3a620b3..e6fa274cc8fb0d8723526224e1d25c699b951cbc 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,29 +31,77 @@ Group
     grpRASBoundaryConditions grpInletBoundaryConditions
 
 Description
-    This boundary condition specifies a velocity inlet profile appropriate
-    for atmospheric boundary layers (ABL).
-
-    See Foam::atmBoundaryLayer for details.
+    This boundary condition provides a log-law type ground-normal inflow
+    boundary condition for 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:
+
+        \f[
+            u = \frac{u^*}{\kappa} \ln \left( \frac{z - d + z_0}{z_0} \right)
+        \f]
+
+        \f[
+            v = w = 0
+        \f]
+
+    where
+    \vartable
+      u        | Ground-normal streamwise flow speed profile    [m/s]
+      v        | Spanwise flow speed                            [m/s]
+      w        | Ground-normal flow speed                       [m/s]
+      u^*      | Friction velocity                              [m/s]
+      \kappa   | von Kármán constant                            [-]
+      z        | Ground-normal coordinate component             [m]
+      d        | Ground-normal displacement height              [m]
+      z_0      | Aerodynamic roughness length                   [m]
+    \endvartable
+
+Usage
+    \heading Required fields
+    \plaintable
+        U   | Velocity                                          [m/s]
+    \endplaintable
 
     Example of the boundary condition specification:
     \verbatim
-    ground
+    inlet
     {
+        // Mandatory entries (unmodifiable)
         type            atmBoundaryLayerInletVelocity;
-        n               (1 0 0);
+
+        // Mandatory (inherited) entries (unmodifiable)
+        flowDir         (1 0 0);
         z               (0 0 1);
         Uref            10.0;
-        Zref            20.0;
+        Zref            0.0;
         z0              uniform 0.1;
-        zGround         uniform 0.0;
+        d               uniform 0.0;
+
+        // Optional (inherited) entries (unmodifiable)
+        kappa           0.41;
+        Cmu             0.09;
+        initABL         true;
+        phi             phi;
+        C1              0.0;
+        C2              1.0;
+
+        // Conditional mandatory (inherited) entries (runtime modifiable)
+        value           uniform (0 0 0);    // when initABL=false
     }
     \endverbatim
 
+    The inherited entries are elaborated in:
+     - \link atmBoundaryLayer.H \endlink
+     - \link inletOutletFvPatchField.H \endlink
+
 See also
-    Foam::atmBoundaryLayer,
-    Foam::atmBoundaryLayerInletKFvPatchScalarField,
-    Foam::atmBoundaryLayerInletEpsilonFvPatchScalarField
+    - Foam::atmBoundaryLayer
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletKFvPatchScalarField
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletEpsilonFvPatchScalarField
+    - Foam::atmBoundaryLayer::atmBoundaryLayerInletOmegaFvPatchScalarField
+    - ExtendedCodeGuide::atmBoundaryLayer::atmBoundaryLayerInletVelocity
 
 SourceFiles
     atmBoundaryLayerInletVelocityFvPatchVectorField.C
@@ -105,7 +154,7 @@ public:
         );
 
         //- Construct by mapping given
-        // atmBoundaryLayerInletVelocityFvPatchVectorField onto a new patch
+        //- atmBoundaryLayerInletVelocityFvPatchVectorField onto a new patch
         atmBoundaryLayerInletVelocityFvPatchVectorField
         (
             const atmBoundaryLayerInletVelocityFvPatchVectorField&,
@@ -143,7 +192,7 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/U b/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/U
index a9827c03f09e5215c8d2e40c353d6f6ef9a710f5..60959f2bbf1133206e2cd0cc0650e0acd2512af3 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/U
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/U
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volVectorField;
-    location    "0";
     object      U;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/epsilon b/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/epsilon
index af522fe72d38541920d6453ef0657700dfae28aa..1cfbc5f0b19740836d7a4287ed4f5f66476d1ab9 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/epsilon
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/epsilon
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      epsilon;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/include/ABLConditions b/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/include/ABLConditions
index 05c1a96e96cccc0955f1678dc7d17b299c695ea7..9d4e3460784e2105b7f95898ea8327dbbc35e3e2 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/include/ABLConditions
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/include/ABLConditions
@@ -11,6 +11,6 @@ Zref                 20;
 zDir                 (0 0 1);
 flowDir              (1 0 0);
 z0                   uniform 0.1;
-zGround              uniform 935.0;
+d                    uniform 0.0;
 
 // ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/nut b/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/nut
index ce7241ea478b442fbdd0f7c2495b1f7d395d08a5..aa772bc020f08c7f183e1070e84051ef0d9eb7dd 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/nut
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/0.orig/nut
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      nut;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/Allrun b/tutorials/incompressible/simpleFoam/turbineSiting/Allrun
index a4566f2fbc132f34ba02e5e57405708e7f136f7c..86a3f1fc77c957dace816317476ec0e6edba89c4 100755
--- a/tutorials/incompressible/simpleFoam/turbineSiting/Allrun
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/Allrun
@@ -3,15 +3,16 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
-# Make dummy 0 directory
-rm -rf 0
-mkdir 0
+restore0Dir
 
 runApplication blockMesh
-cp system/decomposeParDict.hierarchical system/decomposeParDict
+
+cp -f system/decomposeParDict.hierarchical system/decomposeParDict
+
 runApplication decomposePar
 
-cp system/decomposeParDict.ptscotch system/decomposeParDict
+cp -f system/decomposeParDict.ptscotch system/decomposeParDict
+
 runParallel snappyHexMesh -overwrite
 
 find . -iname '*level*' -type f -delete
@@ -20,9 +21,11 @@ find . -iname '*level*' -type f -delete
 restore0Dir -processor
 
 runParallel topoSet
+
 runParallel $(getApplication)
 
 runApplication reconstructParMesh -constant
+
 runApplication reconstructPar
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/resources/dataset/atm-HargreavesWright-2007/Ux-HW-RH-Fig6a b/tutorials/resources/dataset/atm-HargreavesWright-2007/Ux-HW-RH-Fig6a
new file mode 100644
index 0000000000000000000000000000000000000000..29f2449e1166c97a2e1c604e473ec074cc07c793
--- /dev/null
+++ b/tutorials/resources/dataset/atm-HargreavesWright-2007/Ux-HW-RH-Fig6a
@@ -0,0 +1,51 @@
+# 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
+# Fig. 6a
+5.91040498512291, 0.1664466550044139
+6.231131486544684, 0.20004461022864461
+6.510015557466532, 0.19864177284976137
+6.802843831934471, 0.19716879360193929
+7.332743206409284, 0.264925839001954
+7.806875586561379, 0.3681746700877184
+8.25311955962129, 0.4715637849113534
+8.62964251495074, 0.5753036090797181
+9.034053877372374, 0.678903149510198
+9.522159920655517, 0.8877154933568576
+9.996321760392572, 1.0965979790724703
+10.37987573666498, 1.4115700415660655
+10.833170370391153, 1.7966138311348487
+11.23767993142931, 2.2523255536648463
+11.670146638591216, 2.9543755199266144
+12.060820014561454, 3.797480784635006
+12.535286270009722, 5.097911034859102
+12.884235053819857, 6.328550125483808
+13.135525313499084, 7.3836241181413556
+13.38685485262492, 8.579542983638703
+13.582407577566388, 9.775742416611841
+13.791924145777251, 11.042294144135937
+13.931660777087739, 12.09792927174503
+14.057492484298745, 13.294479414062891
+14.18332419150975, 14.491029556380738
+14.330111483486501, 15.828319358734987
+14.448990728647765, 17.09532700840721
+14.61685144373009, 18.99589108931726
+14.756823751720233, 20.896595453965197
+14.875761916051404, 22.37487041289713
+15.001799840357107, 24.310856137623965
+15.113873921393413, 26.176489567799837
+15.239960945007372, 28.28853138357643
+15.345131663302196, 30.400678412156434
+15.457294123093371, 32.58321280622188
+15.527476674321482, 34.237789352744855
+15.611672168127257, 36.13877428486858
+15.702849583567724, 38.07493536426778
+15.766060033022793, 39.72954698172525
+15.843312884640476, 41.73620063941329
+15.906533153957195, 43.42602347508071
+15.969841802028782, 45.43274727463769
+16.03311117065377, 47.29862620135487
+16.089418257367353, 49.199751417216476
\ No newline at end of file
diff --git a/tutorials/resources/dataset/atm-HargreavesWright-2007/epsilon-HW-RH-Fig6c b/tutorials/resources/dataset/atm-HargreavesWright-2007/epsilon-HW-RH-Fig6c
new file mode 100644
index 0000000000000000000000000000000000000000..16ed7e7e210d038fbee0b1a9c830b50c2c3c175d
--- /dev/null
+++ b/tutorials/resources/dataset/atm-HargreavesWright-2007/epsilon-HW-RH-Fig6c
@@ -0,0 +1,50 @@
+# 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
+# Fig. 6c
+0.0020701794999818555, 49.592305129990216
+0.002126271773891171, 48.292814797999284
+0.002164558818180094, 47.16429009397673
+0.0022035486578201443, 45.96736595968585
+0.002263213410047251, 44.873073918499905
+0.0023244936848750017, 43.77878187731396
+0.002387462205022902, 42.54769097559135
+0.0024740215292696124, 41.179834161034535
+0.0025410406291819093, 39.948743259311925
+0.002633152079108662, 38.64928587502344
+0.0027286108031715625, 37.31562877560078
+0.0028401715363583204, 35.70839042895602
+0.0029826690252419655, 33.99858588461122
+0.003132287402529308, 32.391380485668904
+0.0033188290732853176, 30.442210883087387
+0.003595390196335023, 28.390524504359476
+0.003895115615003372, 25.99684097428992
+0.004372395839556802, 23.329691513956814
+0.004908099067927725, 20.799340914160364
+0.0053644197877401765, 18.884486343671554
+0.0057596497045315715, 17.721959610729485
+0.006351093794633881, 16.08073570901643
+0.0069412791686608594, 14.644677150405933
+0.007791105651922957, 13.003519144097766
+0.008667461744937913, 11.704325341428799
+0.009815095638299107, 10.336798003896376
+0.011164037426490142, 9.037686570483537
+0.012698101857369398, 7.977973143009834
+0.01489867869651591, 6.81577588709218
+0.017636199130990377, 5.790410439413613
+0.020876528350971663, 4.86764413713756
+0.02549155925888904, 4.013392582088372
+0.03210854870598559, 3.193456059131897
+0.04246543677826559, 2.476299893941352
+0.05793424483227349, 1.8276584759776568
+0.07527116643710868, 1.4866332818580332
+0.10592731524255083, 1.043305471657888
+0.14840638796815175, 0.805159478411511
+0.21543066611772452, 0.6697444213773878
+0.30316307258220004, 0.5000143322505579
+0.4519488475844894, 0.4330975485920945
+0.6531567674299283, 0.29766601770675294
+0.9565872888053345, 0.19648362350923776
diff --git a/tutorials/resources/dataset/atm-HargreavesWright-2007/k-HW-Fig6b-2500 b/tutorials/resources/dataset/atm-HargreavesWright-2007/k-HW-Fig6b-2500
new file mode 100644
index 0000000000000000000000000000000000000000..edf36ad4d54e26607b66f5a15501c8d8a24f1487
--- /dev/null
+++ b/tutorials/resources/dataset/atm-HargreavesWright-2007/k-HW-Fig6b-2500
@@ -0,0 +1,75 @@
+# 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
+# Fig. 6b - 2500
+1.2900874635568513, 49.577353047723314
+1.2896015549076774, 48.40818131143202
+1.2905733722060253, 47.65161419254281
+1.2905733722060253, 46.68875311964047
+1.291059280855199, 45.313220591832724
+1.292031098153547, 44.212774518335536
+1.292031098153547, 43.249913445433194
+1.292517006802721, 42.14948408131183
+1.2930029154518952, 40.98027892626888
+1.2930029154518952, 39.74231468968016
+1.2930029154518952, 38.64190203493463
+1.2930029154518952, 37.67904096203229
+1.2934888241010691, 36.372284225146146
+1.2934888241010691, 34.92799261579263
+1.2939747327502429, 33.655623774367285
+1.293488824101069, 32.52083993353679
+1.2934888241010691, 31.351651487869663
+1.293488824101069, 30.216850937663335
+1.293002915451895, 29.116454992293637
+1.2925170068027212, 27.500240615011972
+1.292031098153547, 26.124741505955893
+1.2910592808551993, 24.64609541989325
+1.290087463556851, 23.064285647448216
+1.2881438289601554, 21.344957711911654
+1.2866861030126335, 20.141431498911224
+1.2857142857142856, 18.559621726466183
+1.2837706511175897, 17.218560640998398
+1.28134110787172, 15.705576787602446
+1.2784256559766762, 13.986282270817554
+1.2759961127308066, 12.576462103803983
+1.2740524781341107, 11.132237331953803
+1.2725947521865888, 9.894323223492584
+1.2696793002915452, 8.243804497629263
+1.2696793002915452, 6.455633933667784
+1.2701652089407192, 5.01132561493845
+1.2716229348882409, 3.7045354593006437
+1.2779397473275025, 2.569517687208503
+1.2881438289601554, 1.7438572992569021
+1.303692905733722, 1.2962799582399285
+1.324101068999028, 1.0204750007686414
+1.3479105928085517, 0.8821046595097073
+1.3639455782312924, 0.8471653546464779
+1.381438289601555, 0.743400130734166
+1.4086491739552964, 0.6736886147660073
+1.4329446064139941, 0.6040773550528584
+1.4514091350826042, 0.5346666078496725
+1.4747327502429544, 0.533864557809764
+1.4941690962099123, 0.533196182776507
+1.5155490767735667, 0.5324609702399243
+1.5388726919339164, 0.49727102473924845
+1.5573372206025267, 0.46224817299683707
+1.5777453838678328, 0.4271584837511426
+1.5991253644314867, 0.3920353757537569
+1.6214771622934887, 0.39126674446553267
+1.6404275996112732, 0.39061507880812485
+1.6637512147716227, 0.38981302876820223
+1.6875607385811464, 0.38899426935246595
+1.7089407191448005, 0.35387116135510155
+1.727405247813411, 0.3876241005343033
+1.748785228377065, 0.3181130970761288
+1.7759961127308066, 0.3171773720295761
+1.8032069970845481, 0.3162416469830234
+1.8304178814382892, 0.3153059219364778
+1.8624878522837705, 0.3142031031316108
+1.899416909620991, 0.3129331905684367
+1.933430515063168, 0.31176353426025116
+1.9620991253644315, 0.2763897856253976
+1.9839650145772594, 0.27563786371298704
diff --git a/tutorials/resources/dataset/atm-HargreavesWright-2007/k-HW-Fig6b-4000 b/tutorials/resources/dataset/atm-HargreavesWright-2007/k-HW-Fig6b-4000
new file mode 100644
index 0000000000000000000000000000000000000000..e16d231f0ad7332f1389cd8017037c7c144c8421
--- /dev/null
+++ b/tutorials/resources/dataset/atm-HargreavesWright-2007/k-HW-Fig6b-4000
@@ -0,0 +1,85 @@
+# 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
+# Fig. 6b - 4000
+1.274623968947113, 49.543199450483584
+1.2755943716642406, 48.54523075114714
+1.2755943716642406, 47.306414508889745
+1.2760795730228045, 46.03317000781064
+1.2760795730228045, 44.82876532783818
+1.2760795730228045, 43.65877221015063
+1.2765647743813684, 42.38552770907153
+1.2765647743813684, 41.14671146681413
+1.2770499757399323, 39.907878528019964
+1.2770499757399323, 38.703473848047494
+1.2770499757399323, 37.499069168075025
+1.2770499757399323, 36.08819511439299
+1.2770499757399323, 34.952613558990365
+1.2775351770984957, 33.78260374476606
+1.2770499757399323, 32.578215761330355
+1.277049975739932, 31.304987956788032
+1.2765647743813684, 29.99736528649755
+1.2765647743813684, 28.861783731094935
+1.2760795730228045, 27.554161060804454
+1.2755943716642406, 26.315361515083826
+1.274623968947113, 24.973343979045183
+1.274623968947113, 23.768939299072716
+1.2736535662299855, 22.495744887603927
+1.2731683648714216, 21.39459159102301
+1.2721979621542938, 20.052574054984376
+1.2712275594371665, 18.81379120580052
+1.2702571567200387, 17.609419918901587
+1.2692867540029111, 16.301813945147877
+1.2688015526443472, 15.097425961712162
+1.2678311499272197, 13.824231550243383
+1.2663755458515285, 12.585465397596302
+1.2658903444929646, 11.346665851875684
+1.2658903444929646, 10.039026485048431
+1.2658903444929646, 8.834621805075933
+1.2658903444929646, 7.561394000533639
+1.2668607472100921, 6.391367489772541
+1.2692867540029111, 5.2901140139709995
+1.2712275594371665, 4.0856425478514495
+1.2760795730228045, 3.018717151650982
+1.2852983988355167, 2.1236992980442153
+1.2954876273653566, 1.5383521119282548
+1.3066472586123243, 1.3314987178729467
+1.3158660844250365, 1.1247121099647401
+1.328481319747695, 1.0210433131539034
+1.341581756428918, 0.9517693820912356
+1.3575934012615236, 0.8823952718079369
+1.3716642406598738, 0.8130879476717325
+1.3876758854924793, 0.7781253996733568
+1.4056283357593402, 0.7086845032429778
+1.4201843765162545, 0.6393604825699981
+1.432799611838913, 0.6389263726139447
+1.4468704512372634, 0.6384421730475864
+1.4619116933527414, 0.6035130181227615
+1.4793789422610382, 0.5340888182291508
+1.495390587093644, 0.5335378325157123
+1.5080058224163029, 0.5331037225596802
+1.5308102862688013, 0.4979074230465059
+1.549247937894226, 0.4628613923642817
+1.5657447840853955, 0.46229371011408205
+1.5822416302765647, 0.4617260278638611
+1.6026200873362444, 0.4266132110345495
+1.6200873362445416, 0.39160057342586185
+1.6424065987384764, 0.3908325327344002
+1.6618146530810287, 0.39016467126356247
+1.680737506065017, 0.389513506329493
+1.6996603590490054, 0.3888623413954235
+1.7229500242600682, 0.38806090763041823
+1.7462396894711305, 0.3528479115804899
+1.7656477438136826, 0.35218005010965925
+1.7889374090247454, 0.35137861634463263
+1.8131974769529355, 0.35054378950608367
+1.8316351285783603, 0.3499093211087967
+1.8510431829209124, 0.34924145963795894
+1.8685104318292094, 0.3486403843142014
+1.8898592916060166, 0.3134941744113604
+1.913148956817079, 0.312692740646348
+1.950024260067928, 0.31142380385176693
+1.977195536147501, 0.27607723550765684
\ No newline at end of file
diff --git a/tutorials/resources/dataset/atm-HargreavesWright-2007/k-RH-Fig6b b/tutorials/resources/dataset/atm-HargreavesWright-2007/k-RH-Fig6b
new file mode 100644
index 0000000000000000000000000000000000000000..abfb0a52384784eca56126ca448f56d231ccbbc5
--- /dev/null
+++ b/tutorials/resources/dataset/atm-HargreavesWright-2007/k-RH-Fig6b
@@ -0,0 +1,9 @@
+# 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
+# Fig. 6b
+1.3027210884353742, 0.3334523040892563
+1.3027210884353742, 49.164263858422125
diff --git a/tutorials/resources/dataset/atm-HargreavesWright-2007/mut-HW-Fig6d-2500 b/tutorials/resources/dataset/atm-HargreavesWright-2007/mut-HW-Fig6d-2500
new file mode 100644
index 0000000000000000000000000000000000000000..923ed5135b4d402ac3e0bb1663d21bbaaf3e9a64
--- /dev/null
+++ b/tutorials/resources/dataset/atm-HargreavesWright-2007/mut-HW-Fig6d-2500
@@ -0,0 +1,51 @@
+# 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
+# Fig. 6d - 2500
+1.091431306845383, 0.577445652173914
+2.4700813786500753, 1.3247282608695627
+3.848731450454764, 2.0720108695652186
+5.3997127812350385, 2.955163043478258
+6.95069411201532, 3.7703804347826093
+8.788894207754907, 4.789402173913047
+10.971756821445666, 5.944293478260867
+12.924844423168985, 7.03125
+14.935375777884154, 8.118206521739125
+17.290569650550502, 9.408967391304344
+20.105313547151752, 10.9375
+22.69028243178554, 12.398097826086953
+25.102920057443754, 13.620923913043477
+27.63044518908569, 15.08152173913043
+30.617520344662523, 16.67798913043478
+33.3748204882719, 18.172554347826086
+35.84490186692197, 19.49728260869565
+38.659645763523216, 21.0258152173913
+41.18717089516515, 22.452445652173914
+43.48492101483963, 23.743206521739125
+46.127333652465296, 25.13586956521739
+48.310196266156055, 26.324728260869563
+50.37817137386308, 27.445652173913043
+52.561033987553856, 28.702445652173914
+54.51412158927717, 29.789402173913043
+56.58209669698421, 30.978260869565215
+58.42029679272379, 31.9633152173913
+60.258496888463384, 33.05027173913044
+62.15414073719482, 34.06929347826087
+63.8774533269507, 35.05434782608695
+65.5433221637147, 36.03940217391305
+67.15174724748682, 37.02445652173913
+69.10483484921015, 38.11141304347826
+70.94303494494973, 39.19836956521739
+72.32168501675443, 40.047554347826086
+74.27477261847774, 41.202445652173914
+76.11297271421732, 42.255434782608695
+77.95117280995692, 43.34239130434783
+79.55959789372905, 44.327445652173914
+81.22546673049307, 45.346467391304344
+83.1211105792245, 46.56929347826087
+84.67209191000478, 47.452445652173914
+86.16562948779321, 48.40353260869565
+87.8314983245572, 49.49048913043478
\ No newline at end of file
diff --git a/tutorials/resources/dataset/atm-HargreavesWright-2007/mut-HW-Fig6d-4000 b/tutorials/resources/dataset/atm-HargreavesWright-2007/mut-HW-Fig6d-4000
new file mode 100644
index 0000000000000000000000000000000000000000..7c156ab2b04424ac7db09722686700530970b60f
--- /dev/null
+++ b/tutorials/resources/dataset/atm-HargreavesWright-2007/mut-HW-Fig6d-4000
@@ -0,0 +1,65 @@
+# 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
+# Fig. 6d - 4000
+0.7467687888942081, 0.5095108695652186
+1.7807563427477255, 1.188858695652172
+3.1594064145524143, 2.0040760869565233
+4.8827190043082815, 2.955163043478258
+6.376256582096698, 3.702445652173914
+8.099569171852565, 4.721467391304351
+9.593106749640981, 5.502717391304344
+11.431306845380568, 6.555706521739133
+13.039731929152708, 7.404891304347821
+14.360938247965535, 8.118206521739125
+15.797032072762086, 8.865489130434781
+17.63523216850168, 9.884510869565212
+18.898994734322642, 10.631793478260867
+20.392532312111058, 11.345108695652172
+21.886069889899474, 12.194293478260867
+23.37960746768789, 12.975543478260867
+24.873145045476306, 13.824728260869563
+26.65390138822403, 14.70788043478261
+28.20488271900431, 15.557065217391298
+29.640976543800853, 16.37228260869565
+31.134514121589277, 17.119565217391298
+32.62805169937769, 18.002717391304344
+34.29392053614169, 18.851902173913047
+35.78745811393011, 19.63315217391304
+37.28099569171853, 20.482336956521735
+38.831977022498805, 21.33152173913043
+40.44040210627094, 22.18070652173913
+41.93393968405937, 23.09782608695652
+43.48492101483963, 24.048913043478258
+45.03590234561992, 24.864130434782606
+46.471996170416475, 25.679347826086953
+48.022977501196735, 26.4945652173913
+49.516515078985165, 27.377717391304348
+51.01005265677358, 28.226902173913043
+52.618477740545714, 29.144021739130434
+53.99712781235041, 29.99320652173913
+55.605552896122546, 30.842391304347824
+57.21397797989468, 31.725543478260867
+58.59262805169938, 32.57472826086956
+60.02872187649592, 33.42391304347826
+61.57970320727621, 34.341032608695656
+63.01579703207277, 35.224184782608695
+64.6242221158449, 36.17527173913044
+66.0028721876496, 36.99048913043478
+67.55385351842988, 37.87364130434783
+68.93250359023456, 38.82472826086956
+70.31115366203926, 39.63994565217391
+71.86213499281952, 40.59103260869565
+73.29822881761608, 41.474184782608695
+74.67687888942078, 42.42527173913044
+76.11297271421732, 43.34239130434783
+77.54906653901388, 44.25951086956522
+78.92771661081856, 45.21059782608695
+80.30636668262328, 46.09375
+81.68501675442796, 47.01086956521739
+83.06366682623265, 47.961956521739125
+84.38487314504546, 48.91304347826087
+85.59119195787459, 49.76222826086956
\ No newline at end of file
diff --git a/tutorials/resources/dataset/atm-HargreavesWright-2007/mut-RH-Fig6d b/tutorials/resources/dataset/atm-HargreavesWright-2007/mut-RH-Fig6d
new file mode 100644
index 0000000000000000000000000000000000000000..9abbddd6fa713a81168835b63af861f9cdf28154
--- /dev/null
+++ b/tutorials/resources/dataset/atm-HargreavesWright-2007/mut-RH-Fig6d
@@ -0,0 +1,43 @@
+# 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
+# Fig. 6d
+1.2063188128291067, 0.6453804347826093
+3.6189564384873165, 1.936141304347828
+5.801819052178079, 3.0910326086956488
+8.099569171852565, 4.449728260869563
+10.741981809478219, 5.842391304347821
+13.097175682144567, 7.1331521739130395
+15.165150789851609, 8.254076086956516
+17.46290090952609, 9.544836956521735
+19.932982288176163, 10.835597826086953
+22.805169937769264, 12.398097826086953
+25.67735758736238, 13.960597826086953
+28.319770224988034, 15.455163043478258
+30.962182862613687, 16.813858695652172
+33.489707994255625, 18.24048913043478
+35.67257060794638, 19.39538043478261
+37.97032072762087, 20.686141304347824
+40.26807084729536, 21.976902173913043
+42.967927237912875, 23.403532608695652
+45.55289612254667, 24.728260869565215
+47.7932024892293, 26.05298913043478
+50.20584011488751, 27.377717391304348
+52.7908089995213, 28.77038043478261
+55.72044040210628, 30.4008152173913
+59.052178075634274, 32.13315217391304
+61.86692197223552, 33.661684782608695
+64.73910962182862, 35.224184782608695
+67.32407850646241, 36.71875
+70.02393489707994, 38.17934782608695
+72.43657252273816, 39.43614130434783
+74.61943513642892, 40.692934782608695
+77.14696026807086, 42.08559782608695
+79.90426041168023, 43.58016304347826
+82.77644806127334, 45.14266304347826
+85.99329822881762, 46.908967391304344
+88.75059837242699, 48.40353260869565
+91.16323599808521, 49.69429347826087
\ No newline at end of file
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/U b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/U
new file mode 100644
index 0000000000000000000000000000000000000000..47f6eb3608490ca082656c57cf8e773df9706687
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/U
@@ -0,0 +1,66 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            atmBoundaryLayerInletVelocity;
+        #include        "include/ABLConditions"
+        value           uniform (0 0 0);
+    }
+
+    ground
+    {
+        type            noSlip;
+    }
+
+    top
+    {
+    //  (HW:p. 365):
+    //  "In addition, as suggested by RH and often ignored by others, a"
+    //  "constant shear stress of rho*(u^*)^2 was applied at the top boundary."
+    //
+    //  u^* ~ Uref*kappa/ln((Zref+z0)/z0)
+    //  (HW:Table 1):
+    //  Uref = 10 m/s
+    //  Zref = 6 m
+    //  z0 = 0.01 m
+    //  tau = rho*(u^*)^2 = 0.390796574
+        type            fixedShearStress;
+        tau             (0.390796574 0 0);
+        value           uniform (0 0 0);
+    }
+
+    sides
+    {
+        type            symmetry;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/epsilon b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..2e94e69bd366c3514d57e442bd39b20b2a507f4b
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/epsilon
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+    #include            "include/ABLConditions"
+
+    inlet
+    {
+        type            atmBoundaryLayerInletEpsilon;
+        #include        "include/ABLConditions"
+        value           uniform 0;
+    }
+
+    ground
+    {
+        type            epsilonWallFunction;
+        Cmu             $Cmu;
+        kappa           $kappa;
+        value           $internalField;
+    }
+
+    top
+    {
+        type            zeroGradient;
+    }
+
+    sides
+    {
+        type            symmetry;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/include/ABLConditions b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/include/ABLConditions
new file mode 100644
index 0000000000000000000000000000000000000000..020da41e30c023fb61cf2f3f05c7b31beb9d4b3c
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/include/ABLConditions
@@ -0,0 +1,19 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+kappa                0.40;          // (HW:p. 358)
+Cmu                  0.09;          // (HW:p. 358)
+flowDir              (1 0 0);       // (HW:Fig. 1)
+zDir                 (0 0 1);       // (HW:Fig. 1)
+Uref                 10.0;          // (HW:Table 1)
+Zref                 6.0;           // (HW:Table 1)
+z0                   uniform 0.01;  // (HW:Table 1)
+d                    uniform 0.0;
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/k b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/k
new file mode 100644
index 0000000000000000000000000000000000000000..397f0d8ee7165bae8a1948c645ea4bee4b9f8bd1
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/k
@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 1.3;    // (HW:Eq. 6)
+
+boundaryField
+{
+    inlet
+    {
+        type            atmBoundaryLayerInletK;
+        #include        "include/ABLConditions"
+        value           uniform 0;
+    }
+
+    ground
+    {
+        type            kqRWallFunction;
+        value           uniform 0.0;
+    }
+
+    top
+    {
+        type            zeroGradient;
+    }
+
+    sides
+    {
+        type            symmetry;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/nut b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/nut
new file mode 100644
index 0000000000000000000000000000000000000000..36662853d88856f8d738ea543f91c4486ad59cd5
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/nut
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #include            "include/ABLConditions"
+
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    ground
+    {
+        type            nutkAtmRoughWallFunction;
+        kappa           $kappa;
+        Cmu             $Cmu;
+        z0              $z0;
+        value           uniform 0.0;
+    }
+
+    top
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    sides
+    {
+        type            symmetry;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/omega b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/omega
new file mode 100644
index 0000000000000000000000000000000000000000..eeddbff51492633817539367663f3df8b420d675
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/omega
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      omega;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 -1 0 0 0 0];
+
+internalField   uniform 0.0007;
+
+boundaryField
+{
+    #include            "include/ABLConditions"
+
+    inlet
+    {
+        type            atmBoundaryLayerInletOmega;
+        #include        "include/ABLConditions"
+        value           uniform 0;
+    }
+
+    ground
+    {
+        type            omegaWallFunction;
+        Cmu             $Cmu;
+        kappa           $kappa;
+        value           $internalField;
+    }
+
+    top
+    {
+        type            zeroGradient;
+    }
+
+    sides
+    {
+        type            symmetry;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/p b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/p
new file mode 100644
index 0000000000000000000000000000000000000000..2ee320f1673158a27023de0cb5d424eda49c166a
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/0.orig/p
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    ground
+    {
+        type            zeroGradient;
+    }
+
+    top
+    {
+        type            zeroGradient;
+    }
+
+    sides
+    {
+        type            symmetry;
+    }
+
+    outlet
+    {
+        type            uniformFixedValue;
+        uniformValue    constant 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allclean b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..ce6bcd045ea01ff12d72409fa2e67d15a4473946
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allclean
@@ -0,0 +1,12 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase0
+
+rm -rf plots
+rm -f constant/turbulenceProperties
+rm -rf system/atm-HargreavesWright-2007
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..fdca8561bb31a55c3a445cda7851e9dd2eaff5d9
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun
@@ -0,0 +1,21 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+restore0Dir
+
+cp -rf $FOAM_TUTORIALS/resources/dataset/atm-HargreavesWright-2007 system/.
+
+rasModel="kEpsilon"    # Tested options="kOmegaSST","kEpsilon"
+
+sed "s|RAS_MODEL|$rasModel|g" constant/turbulenceProperties.template > \
+    constant/turbulenceProperties
+
+runApplication blockMesh
+
+runApplication renumberMesh -overwrite
+
+runApplication $(getApplication)
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun.par b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun.par
new file mode 100755
index 0000000000000000000000000000000000000000..e4725835dce17412310116ffdda3ae0fe21d944c
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/Allrun.par
@@ -0,0 +1,21 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+restore0Dir
+
+rasModel="kEpsilon"    # Tested options="kOmegaSST","kEpsilon"
+
+sed "s|RAS_MODEL|$rasModel|g" constant/turbulenceProperties.template > \
+    constant/turbulenceProperties
+
+runApplication blockMesh
+
+runApplication decomposePar
+
+runParallel renumberMesh -overwrite
+
+runParallel $(getApplication)
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/README b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/README
new file mode 100644
index 0000000000000000000000000000000000000000..f5aa98553e95a9b5fdc4f2c6061caaa8086109c5
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/README
@@ -0,0 +1,117 @@
+#------------------------------------------------------------------------------
+
+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/constant/transportProperties b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..da462617e0f4626d6d5a9631a12e3b72a073e3aa
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/constant/transportProperties
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              1.5e-05;
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/constant/turbulenceProperties.template b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/constant/turbulenceProperties.template
new file mode 100644
index 0000000000000000000000000000000000000000..14fb1087fa60023658279424aa4c731ef2e996d8
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/constant/turbulenceProperties.template
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType          RAS;
+
+RAS
+{
+    RASModel            RAS_MODEL;
+    turbulence          on;
+    printCoeffs         on;
+
+    RAS_MODELCoeffs     // Valid only for epsilon-based models
+    {
+        Cmu         0.09;
+        C1          1.44;
+        C2          1.92;
+        sigmaEps    1.11; //Original value:1.44
+        // See p. 358:
+        //  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
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/plot b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/plot
new file mode 100755
index 0000000000000000000000000000000000000000..bcf70141f94cf7c42568777ba7be82353ef4c504
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/plot
@@ -0,0 +1,339 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+# Postprocess according to the existence of "epsilon" or "omega"
+baseEpsilonOrOmega="epsilon"    # options: "epsilon", "omega".
+
+# Note: Benchmark data is available for the standard k-epsilon model from:
+#    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
+#    Figure 6.
+#------------------------------------------------------------------------------
+
+plotUxUpstream() {
+    timeDir=$1
+    zMin=$2
+    echo "    Plotting the ground-normal flow speed profile (upstream)."
+
+    outName="plots/Ux_upstream.png"
+    gnuplot<<PLT_UX_UPSTREAM
+    set terminal pngcairo font "helvetica,20" size 1000, 800
+    set xrange [4:18]
+    set yrange [0:50]
+    set grid
+    set key top left
+    set xlabel "Ux [m s^{-1}]"
+    set ylabel "Non-dimensionalised height, z/z_{ref}"
+    set offset .05, .05
+    set output "$outName"
+
+    zRef = 6
+    inp0="$timeDir/x_0mCell_U.xy"
+    inp1="$timeDir/x_0mPatch_U.xy"
+
+    plot \
+        "system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
+            u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
+        "system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
+            u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
+        "system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
+            u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
+        inp0 u 2:((\$1-$zMin)/zRef) t "OF, x=0m (Patch)" w l lw 2 lc rgb "#009E73", \
+        inp1 u 2:((\$1-$zMin)/zRef) t "OF, x=0m (Cell)" w l lw 2 lc rgb "#F0E440"
+PLT_UX_UPSTREAM
+}
+
+
+plotUxMid() {
+    timeDir=$1
+    zMin=$2
+    echo "    Plotting the ground-normal flow speed profile (mid-range)."
+
+    outName="plots/Ux_mid.png"
+    gnuplot<<PLT_UX_MID
+    set terminal pngcairo font "helvetica,20" size 1000, 800
+    set xrange [4:18]
+    set yrange [0:50]
+    set grid
+    set key top left
+    set xlabel "Ux [m s^{-1}]"
+    set ylabel "Non-dimensionalised height, z/z_{ref}"
+    set offset .05, .05
+    set output "$outName"
+
+    zRef = 6
+    inp2="$timeDir/x_2500m_U.xy"
+    inp3="$timeDir/x_4000m_U.xy"
+
+    plot \
+        "system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
+            u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
+        "system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
+            u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
+        "system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
+            u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
+        inp2 u 2:((\$1-$zMin)/zRef) t "OF, x=2500m" w l lw 2 lc rgb "#0072B2", \
+        inp3 u 2:((\$1-$zMin)/zRef) t "OF, x=4000m" w l lw 2 lc rgb "#D55E00"
+PLT_UX_MID
+}
+
+
+plotUxDownstream() {
+    timeDir=$1
+    zMin=$2
+    echo "    Plotting the ground-normal flow speed profile (downstream)."
+
+    outName="plots/Ux_downstream.png"
+    gnuplot<<PLT_UX_DOWNSTREAM
+    set terminal pngcairo font "helvetica,20" size 1000, 800
+    set xrange [4:18]
+    set yrange [0:50]
+    set grid
+    set key top left
+    set xlabel "Ux [m s^{-1}]"
+    set ylabel "Non-dimensionalised height, z/z_{ref}"
+    set offset .05, .05
+    set output "$outName"
+
+    zRef = 6
+    inp4="$timeDir/x_5000mCell_U.xy"
+    inp5="$timeDir/x_5000mPatch_U.xy"
+
+    plot \
+        "system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
+            u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
+        "system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
+            u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
+        "system/benchmark-data-HargreavesWright2007/Ux-HW-RH-Fig6a" \
+            u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
+        inp4 u 2:((\$1-$zMin)/zRef) t "OF, x=5000m (Cell)" w l lw 2 lc rgb "#CC79A7", \
+        inp5 u 2:((\$1-$zMin)/zRef) t "OF, x=5000m (Patch)" w l lw 2 lc rgb "#440154"
+PLT_UX_DOWNSTREAM
+}
+
+
+plotK() {
+    timeDir=$1
+    items=$2
+    seq=$3
+    zMin=$4
+    echo "    Plotting the ground-normal turbulent kinetic energy profile."
+
+    outName="plots/k.png"
+    gnuplot<<PLT_K
+    set terminal pngcairo font "helvetica,20" size 1000, 800
+    set xrange [1:2]
+    set yrange [0:50]
+    set grid
+    set key top right
+    set xlabel "k [m^2 s^{-2}]"
+    set ylabel "Non-dimensionalised height, z/z_{ref}"
+    set offset .05, .05
+    set output "$outName"
+
+    zRef = 6
+    inp0="$timeDir/x_0mCell_$items.xy"
+    inp1="$timeDir/x_0mPatch_$items.xy"
+    inp2="$timeDir/x_2500m_$items.xy"
+    inp3="$timeDir/x_4000m_$items.xy"
+    inp4="$timeDir/x_5000mCell_$items.xy"
+    inp5="$timeDir/x_5000mPatch_$items.xy"
+
+    plot \
+        "system/benchmark-data-HargreavesWright2007/k-RH-Fig6b" \
+            u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
+        "system/benchmark-data-HargreavesWright2007/k-HW-Fig6b-2500" \
+            u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
+        "system/benchmark-data-HargreavesWright2007/k-HW-Fig6b-4000" \
+            u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
+        inp0 u $seq:((\$1-$zMin)/zRef) t "OF, x=0m (Patch)" w l lw 2 lc rgb "#009E73", \
+        inp1 u $seq:((\$1-$zMin)/zRef) t "OF, x=0m (Cell)" w l lw 2 lc rgb "#F0E440", \
+        inp2 u $seq:((\$1-$zMin)/zRef) t "OF, x=2500m" w l lw 2 lc rgb "#0072B2", \
+        inp3 u $seq:((\$1-$zMin)/zRef) t "OF, x=4000m" w l lw 2 lc rgb "#D55E00", \
+        inp4 u $seq:((\$1-$zMin)/zRef) t "OF, x=5000m (Cell)" w l lw 2 lc rgb "#CC79A7", \
+        inp5 u $seq:((\$1-$zMin)/zRef) t "OF, x=5000m (Patch)" w l lw 2 lc rgb "#440154"
+PLT_K
+}
+
+
+plotEpsilon() {
+    timeDir=$1
+    items=$2
+    zMin=$3
+    echo "    Plotting the ground-normal turbulent kinetic"\
+              "energy dissipation rate profile."
+
+    outName="plots/epsilon.png"
+    gnuplot<<PLT_EPSILON
+    set terminal pngcairo font "helvetica,20" size 1000, 800
+    set xrange [0.001:10]
+    set yrange [0:50]
+    set grid
+    set key top right
+    set xlabel "epsilon [m^2 s^{-3}]"
+    set ylabel "Non-dimensionalised height, z/z_{ref}"
+    set offset .05, .05
+    set logscale x
+    set output "$outName"
+
+    zRef = 6
+    inp0="$timeDir/x_0mCell_$items.xy"
+    inp1="$timeDir/x_0mPatch_$items.xy"
+    inp2="$timeDir/x_2500m_$items.xy"
+    inp3="$timeDir/x_4000m_$items.xy"
+    inp4="$timeDir/x_5000mCell_$items.xy"
+    inp5="$timeDir/x_5000mPatch_$items.xy"
+
+    plot \
+        "system/benchmark-data-HargreavesWright2007/epsilon-HW-RH-Fig6c" \
+            u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
+        "system/benchmark-data-HargreavesWright2007/epsilon-HW-RH-Fig6c" \
+            u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
+        "system/benchmark-data-HargreavesWright2007/epsilon-HW-RH-Fig6c" \
+            u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
+        inp0 u 2:((\$1-$zMin)/zRef) t "OF, x=0m (Patch)" w l lw 2 lc rgb "#009E73", \
+        inp1 u 2:((\$1-$zMin)/zRef) t "OF, x=0m (Cell)" w l lw 2 lc rgb "#F0E440", \
+        inp2 u 2:((\$1-$zMin)/zRef) t "OF, x=2500m" w l lw 2 lc rgb "#0072B2", \
+        inp3 u 2:((\$1-$zMin)/zRef) t "OF, x=4000m" w l lw 2 lc rgb "#D55E00", \
+        inp4 u 2:((\$1-$zMin)/zRef) t "OF, x=5000m (Cell)" w l lw 2 lc rgb "#CC79A7", \
+        inp5 u 2:((\$1-$zMin)/zRef) t "OF, x=5000m (Patch)" w l lw 2 lc rgb "#440154"
+PLT_EPSILON
+}
+
+
+plotOmega() {
+    timeDir=$1
+    items=$2
+    zMin=$3
+    echo "    Plotting the ground-normal specific dissipation rate profile."
+
+    outName="plots/omega.png"
+    gnuplot<<PLT_OMEGA
+    set terminal pngcairo font "helvetica,20" size 1000, 800
+    set xrange [0.001:10]
+    set yrange [0:50]
+    set grid
+    set key top right
+    set xlabel "omega [s^{-1}]"
+    set ylabel "Non-dimensionalised height, z/z_{ref}"
+    set offset .05, .05
+    set logscale x
+    set output "$outName"
+
+    zRef = 6
+    inp0="$timeDir/x_0mCell_$items.xy"
+    inp1="$timeDir/x_0mPatch_$items.xy"
+    inp2="$timeDir/x_2500m_$items.xy"
+    inp3="$timeDir/x_4000m_$items.xy"
+    inp4="$timeDir/x_5000mCell_$items.xy"
+    inp5="$timeDir/x_5000mPatch_$items.xy"
+
+    plot \
+        inp0 u 4:((\$1-$zMin)/zRef) t "OF, x=0m (Patch)" w l lw 2 lc rgb "#009E73", \
+        inp1 u 4:((\$1-$zMin)/zRef) t "OF, x=0m (Cell)" w l lw 2 lc rgb "#F0E440", \
+        inp2 u 4:((\$1-$zMin)/zRef) t "OF, x=2500m" w l lw 2 lc rgb "#0072B2", \
+        inp3 u 4:((\$1-$zMin)/zRef) t "OF, x=4000m" w l lw 2 lc rgb "#D55E00", \
+        inp4 u 4:((\$1-$zMin)/zRef) t "OF, x=5000m (Cell)" w l lw 2 lc rgb "#CC79A7", \
+        inp5 u 4:((\$1-$zMin)/zRef) t "OF, x=5000m (Patch)" w l lw 2 lc rgb "#440154"
+PLT_OMEGA
+}
+
+
+plotMut() {
+    timeDir=$1
+    items=$2
+    seq=$3
+    zMin=$4
+    echo "    Plotting the ground-normal turbulent viscosity profile."
+
+    outName="plots/mut.png"
+    gnuplot<<PLT_MUT
+    set terminal pngcairo font "helvetica,20" size 1000, 800
+    set xrange [0:120]
+    set yrange [0:50]
+    set grid
+    set key bottom right
+    set xlabel "mu_t [Pa.s]"
+    set ylabel "Non-dimensionalised height, z/z_{ref}"
+    set offset .05, .05
+    set output "$outName"
+
+    zRef = 6
+    inp0="$timeDir/x_0mCell_$items.xy"
+    inp1="$timeDir/x_0mPatch_$items.xy"
+    inp2="$timeDir/x_2500m_$items.xy"
+    inp3="$timeDir/x_4000m_$items.xy"
+    inp4="$timeDir/x_5000mCell_$items.xy"
+    inp5="$timeDir/x_5000mPatch_$items.xy"
+
+    plot \
+        "system/benchmark-data-HargreavesWright2007/mut-RH-Fig6d" \
+            u 1:2 t "RH" w p ps 2 pt 6 lc rgb "#000000", \
+        "system/benchmark-data-HargreavesWright2007/mut-HW-Fig6d-2500" \
+            u 1:2 t "HW, x=2500m" w p ps 1 pt 5 lc rgb "#E69F00", \
+        "system/benchmark-data-HargreavesWright2007/mut-HW-Fig6d-4000" \
+            u 1:2 t "HW, x=4000m" w p ps 0.5 pt 4 lc rgb "#56B4E9", \
+        inp0 u $seq:((\$1-$zMin)/zRef) t "OF, x=0m (Patch)" w l lw 2 lc rgb "#009E73", \
+        inp1 u $seq:((\$1-$zMin)/zRef) t "OF, x=0m (Cell)" w l lw 2 lc rgb "#F0E440", \
+        inp2 u $seq:((\$1-$zMin)/zRef) t "OF, x=2500m" w l lw 2 lc rgb "#0072B2", \
+        inp3 u $seq:((\$1-$zMin)/zRef) t "OF, x=4000m" w l lw 2 lc rgb "#D55E00", \
+        inp4 u $seq:((\$1-$zMin)/zRef) t "OF, x=5000m (Cell)" w l lw 2 lc rgb "#CC79A7", \
+        inp5 u $seq:((\$1-$zMin)/zRef) t "OF, x=5000m (Patch)" w l lw 2 lc rgb "#440154"
+PLT_MUT
+}
+
+
+#------------------------------------------------------------------------------
+
+# Require gnuplot
+command -v gnuplot >/dev/null || {
+    echo "gnuplot not found - skipping graph creation" 1>&2
+    exit 1
+}
+
+# The latestTime in postProcessing/samples
+timeDir=$(foamListTimes -case postProcessing/samples -latestTime 2>/dev/null)
+[ -n "$timeDir" ] || {
+    echo "No postProcessing/samples found - skipping graph creation" 1>&2
+    exit 2
+}
+timeDir="postProcessing/samples/$timeDir"
+zMin=0
+mkdir -p plots
+
+
+# Postprocess flow speed
+plotUxUpstream $timeDir $zMin
+plotUxMid $timeDir $zMin
+plotUxDownstream $timeDir $zMin
+
+# Postprocess turbulence quantities
+if [ $baseEpsilonOrOmega == "epsilon" ]
+then
+    items="epsilon_k_nut"
+
+    plotEpsilon $timeDir $items $zMin
+    plotK $timeDir $items 3 $zMin
+    plotMut $timeDir $items 4 $zMin
+
+elif [ $baseEpsilonOrOmega == "omega" ]
+then
+    items="k_nut_omega"
+
+    plotK $timeDir $items 2 $zMin
+    plotMut $timeDir $items 3 $zMin
+    plotOmega $timeDir $items $zMin
+
+else
+    echo "Chosen turbulence model is neither epsilon nor omega based." 1>&2
+    exit 2
+fi
+
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/blockMeshDict b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..5d9e682019cfa507153f28442bf10a07f5243631
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/blockMeshDict
@@ -0,0 +1,111 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+    class           dictionary;
+    object          blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale               1;
+
+// x = streamwise
+// y = spanwise
+// z = wall-normal
+
+nx                  500;
+ny                  5;
+nz                  50;
+xMin                0;
+xMax                5000.0;
+yMin                0;
+yMax                100.0;
+zMin                0.0;
+zMax                500.0;
+// blockMesh calculator input:
+// width of start cell = 1.0 (HW:p. 359)
+// number of cells = 50
+// total length = 500
+// blockMesh calculator output:
+// cell-to-cell expansion ratio = 1.076030437 (consistent with 1.076 (HW:Fig.1))
+zTotalExpansion     36.25795062;
+
+vertices
+(
+    ($xMin  $yMin  $zMin)
+    ($xMax  $yMin  $zMin)
+    ($xMax  $yMax  $zMin)
+    ($xMin  $yMax  $zMin)
+    ($xMin  $yMin  $zMax)
+    ($xMax  $yMin  $zMax)
+    ($xMax  $yMax  $zMax)
+    ($xMin  $yMax  $zMax)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) ($nx $ny $nz)  simpleGrading (1 1 $zTotalExpansion)
+);
+
+edges
+(
+);
+
+boundary
+(
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    ground
+    {
+        type wall;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+    top
+    {
+        type patch;
+        faces
+        (
+            (4 5 6 7)
+        );
+    }
+    sides
+    {
+        type symmetry;
+        faces
+        (
+            (1 5 4 0)
+            (3 7 6 2)
+        );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/controlDict b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..f6b7e96e99fe14ecd44c43758543601113c7e473
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/controlDict
@@ -0,0 +1,62 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     simpleFoam;
+
+startFrom       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         5000;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   500;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  12;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+
+functions
+{
+    #include "sampleDict"
+
+    minMax
+    {
+      type          fieldMinMax;
+      libs          (fieldFunctionObjects);
+      writeControl  writeTime;
+      fields        (U);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/decomposeParDict b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..9928e41702395fa91ee9855c780dbfe535420b16
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/decomposeParDict
@@ -0,0 +1,26 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method          hierarchical;
+
+coeffs
+{
+    n           (8 1 1);
+}
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/fvSchemes b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..5e4c45d4b975af167f02015d664be35cd518f3c6
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/fvSchemes
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default             steadyState;
+}
+
+gradSchemes
+{
+    default             Gauss linear;
+}
+
+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;
+    div((nuEff*dev2(T(grad(U)))))    Gauss linear;
+}
+
+laplacianSchemes
+{
+    default             Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default             linear;
+}
+
+snGradSchemes
+{
+    default             uncorrected;
+}
+
+wallDist
+{
+    method              meshWave;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/fvSolution b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..1c87392bce6cb087a93160f88e32899deac0ac0c
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/fvSolution
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver           GAMG;
+        smoother         GaussSeidel;
+        tolerance        1e-8;
+        relTol           0.01;
+    }
+
+    "(U|k|epsilon|omega)"
+    {
+        solver           smoothSolver;
+        smoother         GaussSeidel;
+        tolerance        1e-8;
+        relTol           0.01;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    consistent               true;
+    pRefCell                 0;
+    pRefValue                0;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        U                    0.9;
+        "(k|epsilon|omega)"  0.7;
+    }
+}
+
+cache
+{
+    grad(U);
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/sampleDict b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/sampleDict
new file mode 100644
index 0000000000000000000000000000000000000000..fd8ff84622a8de41344f0d1e5fc5498e51976d50
--- /dev/null
+++ b/tutorials/verificationAndValidation/atmosphericModels/HargreavesWright_2007/system/sampleDict
@@ -0,0 +1,79 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+    class           dictionary;
+    location        system;
+    object          sampleDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// The locations of the sample profiles correspond to:
+// Hargreaves-Wright (2007), Fig.6
+// DOI:10.1016/j.jweia.2006.08.002
+
+samples
+{
+    type                    sets;
+    libs                    (sampling);
+    setFormat               raw;
+    interpolationScheme     cell;
+    fields                  (U k epsilon nut omega);
+    writeControl            writeTime;
+
+    sets
+    (
+        x_0mPatch    // inlet patch face centres
+        {
+            type        face;
+            axis        z;
+            start       (0 50 0);
+            end         (0 50 500);
+        }
+        x_0mCell    // inlet-first cell centres
+        {
+            type        midPoint;
+            axis        z;
+            start       (5.0 50 0);
+            end         (5.0 50 500);
+        }
+        x_2500m
+        {
+            type        face;
+            axis        z;
+            start       (2500 50 0);
+            end         (2500 50 500);
+        }
+        x_4000m
+        {
+            type        face;
+            axis        z;
+            start       (4000 50 0);
+            end         (4000 50 500);
+        }
+        x_5000mCell    // outlet patch face centres
+        {
+            type        face;
+            axis        z;
+            start       (4995 50 0);
+            end         (4995 50 500);
+        }
+        x_5000mPatch    // outlet-first cell centres
+        {
+            type        face;
+            axis        z;
+            start       (5000 50 0);
+            end         (5000 50 500);
+        }
+    );
+}
+
+
+// *********************************************************************** //