diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C
index 1bc235dbeb8694fb9d321c24cde99fc404599397..36a995056504bc68d8a56565a2b58f10a9b0765f 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2014-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,12 +37,13 @@ namespace Foam
 
 atmBoundaryLayer::atmBoundaryLayer(const Time& time, const polyPatch& pp)
 :
+    initABL_(false),
+    kappa_(0.41),
+    Cmu_(0.09),
     time_(time),
     patch_(pp),
     flowDir_(time, "flowDir"),
     zDir_(time, "zDir"),
-    kappa_(0.41),
-    Cmu_(0.09),
     Uref_(time, "Uref"),
     Zref_(time, "Zref"),
     z0_(),
@@ -57,12 +58,13 @@ atmBoundaryLayer::atmBoundaryLayer
     const dictionary& dict
 )
 :
+    initABL_(dict.getOrDefault<Switch>("initABL", true)),
+    kappa_(dict.getOrDefault<scalar>("kappa", 0.41)),
+    Cmu_(dict.getOrDefault<scalar>("Cmu", 0.09)),
     time_(time),
     patch_(pp),
     flowDir_(TimeFunction1<vector>(time, "flowDir", dict)),
     zDir_(TimeFunction1<vector>(time, "zDir", dict)),
-    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
-    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
     Uref_(TimeFunction1<scalar>(time, "Uref", dict)),
     Zref_(TimeFunction1<scalar>(time, "Zref", dict)),
     z0_(PatchFunction1<scalar>::New(pp, "z0", dict)),
@@ -77,12 +79,13 @@ atmBoundaryLayer::atmBoundaryLayer
     const fvPatchFieldMapper& mapper
 )
 :
+    initABL_(abl.initABL_),
+    kappa_(abl.kappa_),
+    Cmu_(abl.Cmu_),
     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_)),
@@ -92,12 +95,13 @@ atmBoundaryLayer::atmBoundaryLayer
 
 atmBoundaryLayer::atmBoundaryLayer(const atmBoundaryLayer& abl)
 :
+    initABL_(abl.initABL_),
+    kappa_(abl.kappa_),
+    Cmu_(abl.Cmu_),
     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_)),
@@ -110,7 +114,7 @@ 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)
@@ -128,7 +132,7 @@ 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)
@@ -149,6 +153,7 @@ tmp<scalarField> atmBoundaryLayer::Ustar(const scalarField& z0) const
     const scalar Uref = Uref_.value(t);
     const scalar Zref = Zref_.value(t);
 
+    // (derived from RH:Eq. 6, HW:Eq. 5)
     return kappa_*Uref/(log((Zref + z0)/z0));
 }
 
@@ -177,6 +182,7 @@ tmp<vectorField> atmBoundaryLayer::U(const vectorField& p) const
     const scalarField zGround(zGround_->value(t));
     const scalarField z0(max(z0_->value(t), ROOTVSMALL));
 
+    // (RH:Eq. 6, HW:Eq. 5)
     scalarField Un((Ustar(z0)/kappa_)*log(((zDir() & p) - zGround + z0)/z0));
 
     return flowDir()*Un;
@@ -188,6 +194,7 @@ tmp<scalarField> atmBoundaryLayer::k(const vectorField& p) const
     const scalar t = time_.timeOutputValue();
     const scalarField z0(max(z0_->value(t), ROOTVSMALL));
 
+    // (RH:Eq. 7, HW:Eq. 6)
     return sqr(Ustar(z0))/sqrt(Cmu_);
 }
 
@@ -198,19 +205,21 @@ tmp<scalarField> atmBoundaryLayer::epsilon(const vectorField& p) const
     const scalarField zGround(zGround_->value(t));
     const scalarField z0(max(z0_->value(t), ROOTVSMALL));
 
+    // (RH:Eq. 8, HW:Eq. 7)
     return pow3(Ustar(z0))/(kappa_*((zDir() & p) - zGround + 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_);
+    flowDir_.writeData(os);
+    zDir_.writeData(os);
     Uref_.writeData(os);
     Zref_.writeData(os);
+    z0_->writeData(os) ;
     zGround_->writeData(os);
 }
 
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H
index 888801476f99f22a61f31c93b5d596a6faa1e97f..812b9f1ef83541a0ec00c49e4271e9ff1c7fd618 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,101 @@ 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 appropriate ground-normal inlet profiles for wind velocity
+    and turbulence quantities for homogeneous two-dimensional atmospheric
+    boundary layer (ABL) modelling.
 
-    The profile is derived from the friction velocity, flow direction and
-    "vertical" direction:
+    The ground-normal profile expressions:
 
         \f[
-            U = \frac{U^*}{\kappa} ln\left(\frac{z - z_g + z_0}{z_0}\right)
+            u = \frac{u^*}{\kappa} \ln \left( \frac{z - z_g + z_0}{z_0} \right)
         \f]
 
         \f[
-            k = \frac{(U^*)^2}{\sqrt{C_{\mu}}}
+            k = \frac{(u^*)^2}{\sqrt{C_\mu}}
         \f]
 
         \f[
-            \epsilon = \frac{(U^*)^3}{\kappa(z - z_g + z_0)}
+            \epsilon = \frac{(u^*)^3}{\kappa (z - z_g + z_0)}
         \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)}
+            u^* =
+                \frac{U_{ref} \kappa}{\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]
+      k        | Ground-normal turbulent kinetic energy (TKE) profile [m^2/s^2]
+      \epsilon | Ground-normal TKE 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]
+      z_g      | Minimum ground-normal coordinate of the domain       [m]
+      z_0      | Surface roughness length                             [m]
+      U_{ref}  | Reference mean wind speed at Z_{ref}                 [m/s]
+      Z_{ref}  | Reference height                                     [m]
     \endvartable
-
-    Use in the atmBoundaryLayerInletVelocity, atmBoundaryLayerInletK and
-    atmBoundaryLayerInletEpsilon boundary conditions.
+    and
 
     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
+        Homogeneous two-dimensional ABL 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 ABL 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
+    \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 \\
+                   | TimeFunction1<scalar>             | yes | -
+      Zref         | Reference height                [m]  \\
+                   | TimeFunction1<scalar>             | yes | -
+      z0           | Surface roughness length        [m]  \\
+                   | PatchFunction1<scalar>            | yes | -
+      zGround      | Minimum ground-normal coordinate  [m]  \\
+                   | 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 | -
     \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 (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 side 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 without the need for re-meshing.
+
+See also
+    Foam::atmBoundaryLayer::atmBoundaryLayerInletVelocity
+    Foam::atmBoundaryLayer::atmBoundaryLayerInletK
+    Foam::atmBoundaryLayer::atmBoundaryLayerInletEpsilon
 
 SourceFiles
     atmBoundaryLayer.C
@@ -133,7 +150,24 @@ 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_;
 
         //- Reference to the time database
         const Time& time_;
@@ -144,22 +178,16 @@ class atmBoundaryLayer
         //- 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
         TimeFunction1<scalar> Uref_;
 
         //- Reference height
         TimeFunction1<scalar> Zref_;
 
-        //- Surface roughness height
+        //- Surface roughness length
         autoPtr<PatchFunction1<scalar>> z0_;
 
         //- Minimum coordinate value in z direction
@@ -193,14 +221,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
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
index 9063fa309abcb62d02e40509068005b706ad81df..da71ad401f33f190d4232d49f47fcfc4af81a774 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,19 +62,20 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
     inletOutletFvPatchScalarField(p, iF),
     atmBoundaryLayer(iF.time(), p.patch(), dict)
 {
-    phiName_ = dict.lookupOrDefault<word>("phi", "phi");
+    phiName_ = dict.getOrDefault<word>("phi", "phi");
 
     refValue() = epsilon(patch().Cf());
     refGrad() = 0;
     valueFraction() = 1;
 
-    if (dict.found("value"))
+    if (!initABL_)
     {
         scalarField::operator=(scalarField("value", dict, p.size()));
     }
     else
     {
         scalarField::operator=(refValue());
+        initABL_ = false;
     }
 }
 
@@ -147,8 +149,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..c61d035bf854366a4bf688ea85d6abf7a411bb86 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,27 +31,82 @@ 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 ground-normal inlet profile for turbulent
+    kinetic energy dissipation rate, i.e. \c epsilon, for homogeneous
+    two-dimensional atmospheric boundary layer modelling.
+    Refer to Foam::atmBoundaryLayer for further details.
+
+    The ground-normal \c epsilon profile expression:
+
+        \f[
+            \epsilon = \frac{(u^*)^3}{\kappa (z - z_g + z_0)}
+        \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      | Surface roughness length                       [m]
+        z_g      | Minimum ground-normal coordinate of the domain [m]
+    \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 entries (runtime modifiable)
+        flowDir         (1 0 0);
         z               (0 0 1);
         Uref            10.0;
         Zref            20.0;
         z0              uniform 0.1;
         zGround         uniform 0.0;
+
+        // Optional entries (runtime modifiable)
+        kappa           0.41;
+        Cmu             0.09;
+        initABL         true;
+
+        // Conditional mandatory entries (runtime modifiable)
+        value           uniform 0;    // when initABL=false
     }
     \endverbatim
 
+    \table
+      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 \\
+                   | TimeFunction1<scalar>             | yes | -
+      Zref         | Reference height                [m]  \\
+                   | TimeFunction1<scalar>             | yes | -
+      z0           | Surface roughness length        [m]  \\
+                   | PatchFunction1<scalar>            | yes | -
+      zGround      | Minimum ground-normal coordinate  [m]  \\
+                   | 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 | -
+    \endtable
+
 See also
-    Foam::atmBoundaryLayer,
-    Foam::atmBoundaryLayerInletVelocityFvPatchVectorField,
+    Foam::atmBoundaryLayer
+    Foam::atmBoundaryLayerInletVelocityFvPatchVectorField
     Foam::atmBoundaryLayerInletKFvPatchScalarField
 
 SourceFiles
@@ -104,7 +160,7 @@ public:
         );
 
         //- Construct by mapping given
-        //  atmBoundaryLayerInletEpsilonFvPatchScalarField onto a new patch
+        //- atmBoundaryLayerInletEpsilonFvPatchScalarField onto a new patch
         atmBoundaryLayerInletEpsilonFvPatchScalarField
         (
             const atmBoundaryLayerInletEpsilonFvPatchScalarField&,
@@ -142,7 +198,7 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C
index 4c04f7c4a82688b9ee190ed4edd857fc346b1bce..7fc3d768e3f0d38032aa3925db54a86c7327d149 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,19 +62,20 @@ atmBoundaryLayerInletKFvPatchScalarField
     inletOutletFvPatchScalarField(p, iF),
     atmBoundaryLayer(iF.time(), p.patch(), dict)
 {
-    phiName_ = dict.lookupOrDefault<word>("phi", "phi");
+    phiName_ = dict.getOrDefault<word>("phi", "phi");
 
     refValue() = k(patch().Cf());
     refGrad() = 0;
     valueFraction() = 1;
 
-    if (dict.found("value"))
+    if (!initABL_)
     {
         scalarField::operator=(scalarField("value", dict, p.size()));
     }
     else
     {
         scalarField::operator=(refValue());
+        initABL_ = false;
     }
 }
 
@@ -147,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..1c47b927b5a7b10c19d79a255b211713f2dde7cc 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,27 +31,79 @@ 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 ground-normal inlet profile for turbulent
+    kinetic energy, i.e. \c k, for homogeneous two-dimensional atmospheric
+    boundary layer modelling.
+    Refer to Foam::atmBoundaryLayer for further details.
+
+    The ground-normal \c k profile expression:
+
+        \f[
+            k = \frac{(u^*)^2}{\sqrt{C_\mu}}
+        \f]
+
+    where
+    \vartable
+        k        | Ground-normal turbulent kinetic energy profile [m^2/s^3]
+        u^*      | Friction velocity                              [m/s]
+        C_\mu    | Empirical model constant                       [-]
+    \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;
+
+        // Mandatory entries (runtime modifiable)
+        flowDir         (1 0 0);
         z               (0 0 1);
         Uref            10.0;
         Zref            20.0;
         z0              uniform 0.1;
         zGround         uniform 0.0;
+
+        // Optional entries (runtime modifiable)
+        kappa           0.41;
+        Cmu             0.09;
+        initABL         true;
+
+        // Conditional mandatory entries (runtime modifiable)
+        value           uniform 0;    // when initABL=false
     }
     \endverbatim
 
+    \table
+      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 \\
+                   | TimeFunction1<scalar>             | yes | -
+      Zref         | Reference height                [m]  \\
+                   | TimeFunction1<scalar>             | yes | -
+      z0           | Surface roughness length        [m]  \\
+                   | PatchFunction1<scalar>            | yes | -
+      zGround      | Minimum ground-normal coordinate  [m]  \\
+                   | 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 | -
+    \endtable
+
 See also
-    Foam::atmBoundaryLayer,
-    Foam::atmBoundaryLayerInletVelocityFvPatchVectorField,
+    Foam::atmBoundaryLayer
+    Foam::atmBoundaryLayerInletVelocityFvPatchVectorField
     Foam::atmBoundaryLayerInletEpsilonFvPatchScalarField
 
 SourceFiles
@@ -104,7 +157,7 @@ public:
         );
 
         //- Construct by mapping given
-        //  atmBoundaryLayerInletKFvPatchScalarField onto a new patch
+        //- atmBoundaryLayerInletKFvPatchScalarField onto a new patch
         atmBoundaryLayerInletKFvPatchScalarField
         (
             const atmBoundaryLayerInletKFvPatchScalarField&,
@@ -142,7 +195,7 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
index c57ef0c13c642831ee930845709a497238db1880..077a47391b1d05112207ad2a8796d484c99e013f 100644
--- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
+++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,19 +62,20 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
     inletOutletFvPatchVectorField(p, iF),
     atmBoundaryLayer(iF.time(), p.patch(), dict)
 {
-    phiName_ = dict.lookupOrDefault<word>("phi", "phi");
+    phiName_ = dict.getOrDefault<word>("phi", "phi");
 
     refValue() = U(patch().Cf());
     refGrad() = Zero;
     valueFraction() = 1;
 
-    if (dict.found("value"))
+    if (!initABL_)
     {
         vectorField::operator=(vectorField("value", dict, p.size()));
     }
     else
     {
         vectorField::operator=(refValue());
+        initABL_ = false;
     }
 }
 
@@ -147,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..4a826562dd836c893c979363e4c70104c6bd28cd 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,25 +31,79 @@ 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 ground-normal inlet profile for
+    streamwise component of flow velocity, i.e. \c U, for homogeneous
+    two-dimensional atmospheric boundary layer modelling.
+    Refer to Foam::atmBoundaryLayer for further details.
+
+    The ground-normal \c U profile expression:
+
+        \f[
+            u = \frac{u^*}{\kappa} \ln \left( \frac{z - z_g + z_0}{z_0} \right)
+        \f]
+
+    where
+    \vartable
+        u        | Ground-normal streamwise flow speed profile    [m/s]
+        u^*      | Friction velocity                              [m/s]
+        \kappa   | von Kármán constant                            [-]
+        z        | Ground-normal coordinate component             [m]
+        z_g      | Minimum ground-normal coordinate of the domain [m]
+        z_0      | Surface 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 entries (runtime modifiable)
+        flowDir         (1 0 0);
         z               (0 0 1);
         Uref            10.0;
         Zref            20.0;
         z0              uniform 0.1;
         zGround         uniform 0.0;
+
+        // Optional entries (runtime modifiable)
+        kappa           0.41;
+        Cmu             0.09;
+        initABL         true;
+
+        // Conditional mandatory entries (runtime modifiable)
+        value           uniform (0 0 0);    // when initABL=false
     }
     \endverbatim
 
+    \table
+      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 \\
+                   | TimeFunction1<scalar>             | yes | -
+      Zref         | Reference height                [m]  \\
+                   | TimeFunction1<scalar>             | yes | -
+      z0           | Surface roughness length        [m]  \\
+                   | PatchFunction1<scalar>            | yes | -
+      zGround      | Minimum ground-normal coordinate  [m]  \\
+                   | 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 | -
+    \endtable
+
 See also
     Foam::atmBoundaryLayer,
     Foam::atmBoundaryLayerInletKFvPatchScalarField,
@@ -105,7 +160,7 @@ public:
         );
 
         //- Construct by mapping given
-        // atmBoundaryLayerInletVelocityFvPatchVectorField onto a new patch
+        //- atmBoundaryLayerInletVelocityFvPatchVectorField onto a new patch
         atmBoundaryLayerInletVelocityFvPatchVectorField
         (
             const atmBoundaryLayerInletVelocityFvPatchVectorField&,
@@ -143,7 +198,7 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();