WIP: ENH: evapotranspirationHeatTransfer: New atmospheric-model fvOption

WIP: The units of the main governing eq are inconsistent - see EP1950

Applies sources on temperature ('T' - incompressible)
or energy ('h'/'e' - compressible) equation to incorporate
evapotranspiration heat-transfer effects from the specified
plant canopy. Heat transfer is usually calculated based on
empirical relations between plants and solar radiation.

Two submodels to incorporate heat transfer effects
- 'tree': specified tree canopy - uses empirical relations between solar
radiation and evapotranspiration.
- 'grass': specified grass canopy - uses Pemnan-Monteith Equation model.

src/atmosphericModels/Make/files

@@ -33,7 +33,15 @@ fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSource.C
 fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.C
 fvOptions/atmPlantCanopyTSource/atmPlantCanopyTSource.C
 fvOptions/atmNutSource/atmNutSource.C
+fvOptions/treeTurbulence/treeTurbulence.C
 
+etht = fvOptions/evapotranspirationHeatTransfer
+ethtModels = $(etht)/evapotranspirationHeatTransferModels
+$(etht)/evapotranspirationHeatTransfer.C
+$(ethtModels)/evapotranspirationHeatTransferModel/evapotranspirationHeatTransferModel.C
+$(ethtModels)/evapotranspirationHeatTransferModel/evapotranspirationHeatTransferModelNew.C
+$(ethtModels)/tree/tree.C
+$(ethtModels)/grass/grass.C
 
 /* functionObjects */
 functionObjects/ObukhovLength/ObukhovLength.C

src/atmosphericModels/Make/options

@@ -8,6 +8,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
@@ -18,6 +19,7 @@ LIB_LIBS = \
     -lfvOptions \
     -lmeshTools \
     -lsurfMesh \
+    -lradiationModels \
     -lfluidThermophysicalModels \
     -lsolidThermo \
     -lturbulenceModels \

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransfer.C

@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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 "evapotranspirationHeatTransfer.H"
+#include "evapotranspirationHeatTransferModel.H"
+#include "fvMesh.H"
+#include "fvMatrix.H"
+#include "basicThermo.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(evapotranspirationHeatTransfer, 0);
+    addToRunTimeSelectionTable
+    (
+        option,
+        evapotranspirationHeatTransfer,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::evapotranspirationHeatTransfer::evapotranspirationHeatTransfer
+(
+    const word& sourceName,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    fv::cellSetOption(sourceName, modelType, dict, mesh),
+    ethtModelPtr_()
+{
+    // Set the field name to that of the energy
+    // field from which the temperature is obtained
+
+    const auto& thermo = mesh_.lookupObject<basicThermo>(basicThermo::dictName);
+
+    fieldNames_.resize(1, thermo.he().name());
+
+    fv::option::resetApplied();
+
+    Info<< "    Applying evapotranspirationHeatTransfer to: " << fieldNames_[0]
+        << endl;
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::fv::evapotranspirationHeatTransfer::~evapotranspirationHeatTransfer()
+{}  // evapotranspirationHeatTransferModel was forward declared
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fv::evapotranspirationHeatTransfer::addSup
+(
+    fvScalarMatrix& eqn,
+    const label fieldi
+)
+{
+    if (this->V() < VSMALL)
+    {
+        return;
+    }
+
+    const scalar V = this->V();
+    const scalarField& Vcells = mesh_.V();
+
+    const volScalarField& he = eqn.psi();
+    scalarField& heSource = eqn.source();
+
+    // Calculate evapotranspiration heat transfer rate per volume [J/s/m^3]
+    const scalarField Q(ethtModelPtr_->Q(cells_)/V);
+
+    if (he.dimensions() == dimEnergy/dimMass)
+    {
+        forAll(cells_, i)
+        {
+            const label celli = cells_[i];
+
+            heSource[celli] += Q[i]*Vcells[celli];
+        }
+    }
+    else if (he.dimensions() == dimTemperature)
+    {
+        const auto& thermo =
+            mesh_.lookupObject<basicThermo>(basicThermo::dictName);
+
+        // Calculate density*heat capacity at constant pressure/volume
+        const volScalarField rhoCpv(thermo.rho()*thermo.Cpv());
+
+        // heSource [K m^3/s] = [J/s/m^3] * m^3 / [kg/m^3] / [J/kg/K]
+        forAll(cells_, i)
+        {
+            const label celli = cells_[i];
+
+            heSource[celli] += Q[i]*Vcells[celli]*rhoCpv[i];
+        }
+    }
+}
+
+
+void Foam::fv::evapotranspirationHeatTransfer::addSup
+(
+    const volScalarField& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+)
+{
+    addSup(eqn, fieldi);
+}
+
+
+bool Foam::fv::evapotranspirationHeatTransfer::read(const dictionary& dict)
+{
+    if (!fv::cellSetOption::read(dict))
+    {
+        return false;
+    }
+
+    if (selectionMode_ != selectionModeType::smCellZone)
+    {
+        FatalIOErrorInFunction(dict)
+            << "evapotranspirationHeatTransfer requires "
+            << selectionModeTypeNames_[selectionModeType::smCellZone]
+            << exit(FatalIOError);
+    }
+
+    ethtModelPtr_.reset
+    (
+        evapotranspirationHeatTransferModel::New(dict, mesh_)
+    );
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransfer.H

@@ -0,0 +1,183 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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/>.
+
+Class
+    Foam::fv::evapotranspirationHeatTransfer
+
+Description
+    Applies sources on temperature (\c T - incompressible)
+    or energy (\c h/e - compressible) equation to incorporate
+    evapotranspiration heat-transfer effects from the specified
+    plant canopy. Heat transfer is usually calculated based on
+    empirical relations between plants and solar radiation.
+
+Usage
+    Example by using \c constant/fvOptions:
+    \verbatim
+    evapotranspirationHeatTransfer1
+    {
+        // Mandatory entries
+        type                evapotranspirationHeatTransfer;
+        model               <word>;
+
+        // Model-specific entries
+        ...
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property | Description                         | Type | Reqd | Deflt
+      type     | Type name: evapotranspirationHeatTransfer  | word | yes  | -
+      model    | Name of heat-transfer model        | word  | yes  | -
+    \endtable
+
+    Options for the \c model entry:
+    \verbatim
+      tree   | Regression-based heat-transfer model for trees
+      grass  | Heat-transfer model for grass
+    \endverbatim
+
+    The inherited entries are elaborated in:
+      - \link fvOption.H \endlink
+      - \link cellSetOption.H \endlink
+
+Note
+  - Evapotranspiration mass transfer is not modelled.
+
+SourceFiles
+    evapotranspirationHeatTransfer.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_fv_evapotranspirationHeatTransfer_H
+#define Foam_fv_evapotranspirationHeatTransfer_H
+
+#include "cellSetOption.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class evapotranspirationHeatTransferModel;
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                Class evapotranspirationHeatTransfer Declaration
+\*---------------------------------------------------------------------------*/
+
+class evapotranspirationHeatTransfer
+:
+    public fv::cellSetOption
+{
+    // Private Data
+
+        //- Runtime-selectable evapotranspiration heat transfer model
+        autoPtr<evapotranspirationHeatTransferModel> ethtModelPtr_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("evapotranspirationHeatTransfer");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        evapotranspirationHeatTransfer
+        (
+            const word& sourceName,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- No copy construct
+        evapotranspirationHeatTransfer(const evapotranspirationHeatTransfer&) =
+            delete;
+
+        //- No copy assignment
+        void operator=(const evapotranspirationHeatTransfer&) = delete;
+
+
+    //- Destructor
+    virtual ~evapotranspirationHeatTransfer();
+
+
+    // Member Functions
+
+        //- Add explicit contribution to energy equation
+        //- for incompressible flow computations
+        virtual void addSup
+        (
+            fvScalarMatrix& eqn,
+            const label fieldi
+        );
+
+        //- Add explicit contribution to energy equation
+        //- for compressible flow computations
+        virtual void addSup
+        (
+            const volScalarField& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        );
+
+        //- Add explicit contribution to energy equation
+        //- for multiphase flow computations
+        virtual void addSup
+        (
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        )
+        {
+            NotImplemented;
+        }
+
+        //- Read dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/evapotranspirationHeatTransferModel/evapotranspirationHeatTransferModel.C

@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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 "evapotranspirationHeatTransferModel.H"
+#include "radiationModel.H"
+#include "solarLoadBase.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(evapotranspirationHeatTransferModel, 0);
+    defineRunTimeSelectionTable
+    (
+        evapotranspirationHeatTransferModel,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::volScalarField& Foam::evapotranspirationHeatTransferModel::getOrReadField
+(
+    const word& fieldName
+) const
+{
+    auto* ptr = mesh_.getObjectPtr<volScalarField>(fieldName);
+
+    if (!ptr)
+    {
+        ptr = new volScalarField
+        (
+            IOobject
+            (
+                fieldName,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh_
+        );
+        mesh_.objectRegistry::store(ptr);
+    }
+
+    return *ptr;
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModel::q() const
+{
+    // Retrieve solar-load information
+    const auto& base =
+        mesh().lookupObject<radiation::solarLoadBase>("solarLoadBase");
+    const solarCalculator& solar = base.solarCalculatorRef();
+
+    // Retrieve internal & boundary faces directly hit by solar rays
+    const faceShading& shade = base.faceShadingRef();
+    const labelList& hitFacesId = shade.rayStartFaces();
+
+    // Retrieve face zone information
+    const faceZone& zone = mesh_.faceZones()[faceZoneId_];
+    const vectorField& faceNormals = zone().faceNormals();
+    const scalarField& faceAreas = zone().magFaceAreas();
+
+    // Retrieve direct solar load [W/m^2]
+    const vector directSolarRad(solar.directSolarRad()*solar.direction());
+
+    // Identify zone faces within mesh
+    bitSet isZoneFace(mesh_.nFaces());
+    isZoneFace.set(zone);
+
+    // Calculate area-averaged incident solar load
+    // Assuming hit faces are updated by solarLoad
+    scalar q = 0;
+    scalar totalFaceArea = 0;
+
+    for (const label facei : hitFacesId)
+    {
+        if (isZoneFace[facei])
+        {
+            const label localFacei = zone.whichFace(facei);
+
+            const vector& faceNormal = faceNormals[localFacei];
+            const scalar faceArea = faceAreas[localFacei];
+
+            const scalar qIncident = directSolarRad & faceNormal;
+
+            q += qIncident*faceArea;
+            totalFaceArea += faceArea;
+        }
+    }
+    reduce(q, sumOp<scalar>());
+    reduce(totalFaceArea, sumOp<scalar>());
+    q /= totalFaceArea;
+
+
+    // Sum diffusive solar loads
+    switch(solar.sunLoadModel())
+    {
+        case solarCalculator::mSunLoadFairWeatherConditions:
+        case solarCalculator::mSunLoadTheoreticalMaximum:
+        {
+            // Calculate diffusive radiance
+            // contribution from sky and ground
+            tmp<scalarField> tdiffuseSolarRad =
+                solar.diffuseSolarRad(faceNormals);
+            const scalarField& diffuseSolarRad = tdiffuseSolarRad.cref();
+
+            // Calculate area-averaged diffusive solar load
+            scalar meanDiffuseSolarRad = 0;
+            forAll(faceAreas, i)
+            {
+                meanDiffuseSolarRad += diffuseSolarRad[i]*faceAreas[i];
+            }
+            meanDiffuseSolarRad /= totalFaceArea;
+
+            q += meanDiffuseSolarRad;
+        }
+        break;
+
+        case solarCalculator::mSunLoadConstant:
+        case solarCalculator::mSunLoadTimeDependent:
+        {
+            q += solar.diffuseSolarRad();
+        }
+        break;
+    }
+
+    return q;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::evapotranspirationHeatTransferModel::evapotranspirationHeatTransferModel
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::evapotranspirationHeatTransferModel::read(const dictionary& dict)
+{
+    faceZoneId_ = mesh_.faceZones().findZoneID(dict.get<word>("faceZone"));
+
+    if (faceZoneId_ < 0)
+    {
+        const word faceZoneName(mesh_.faceZones()[faceZoneId_].name());
+
+        FatalIOErrorInFunction(dict)
+            << type() << ' ' << typeName << ": "
+            << "    No matching face zone: " << faceZoneName  << nl
+            << "    Known face zones: "
+            << flatOutput(mesh_.faceZones().names()) << nl
+            << exit(FatalIOError);
+
+        return false;
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/evapotranspirationHeatTransferModel/evapotranspirationHeatTransferModel.H

@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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/>.
+
+Class
+    Foam::evapotranspirationHeatTransferModel
+
+Description
+    Base class for evapotranspiration models
+    to handle general evapotranspiration characteristics.
+
+SourceFiles
+    evapotranspirationHeatTransferModel.C
+    evapotranspirationHeatTransferModelNew.C
+    evapotranspirationHeatTransferModelTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_evapotranspirationHeatTransferModel_H
+#define Foam_evapotranspirationHeatTransferModel_H
+
+#include "dictionary.H"
+#include "HashSet.H"
+#include "volFields.H"
+#include "runTimeSelectionTables.H"
+#include "OFstream.H"
+#include "coordinateSystem.H"
+#include "writeFile.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+            Class evapotranspirationHeatTransferModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class evapotranspirationHeatTransferModel
+{
+    // Private Data
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+        //- Identifier of specified face zone
+        label faceZoneId_;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Return requested field from the object registry
+        //- or read+register the field to the object registry
+        volScalarField& getOrReadField(const word& fieldName) const;
+
+        //- Area-averaged heat flux due to short-wave
+        //- solar rays on face zone [W/m^2]
+        //  Short-wave: direct and diffusive solar rays
+        scalar q() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("evapotranspirationHeatTransferModel");
+
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            evapotranspirationHeatTransferModel,
+            dictionary,
+            (
+                const dictionary& dict,
+                const fvMesh& mesh
+            ),
+            (dict, mesh)
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected evapotranspiration model
+        static autoPtr<evapotranspirationHeatTransferModel> New
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        evapotranspirationHeatTransferModel
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- No copy construct
+        evapotranspirationHeatTransferModel
+        (
+            const evapotranspirationHeatTransferModel&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const evapotranspirationHeatTransferModel&) = delete;
+
+
+    //- Destructor
+    virtual ~evapotranspirationHeatTransferModel() = default;
+
+
+    // Member Functions
+
+    // Access
+
+        //- Return const reference to the mesh
+        const fvMesh& mesh() const noexcept
+        {
+            return mesh_;
+        }
+
+        //- Return the face-zone identifier
+        label faceZoneId() const noexcept
+        {
+            return faceZoneId_;
+        }
+
+
+    // Evaluation
+
+        //- Return heat-transfer rate for speficied cells [J/s]
+        virtual tmp<scalarField> Q(const labelList& cells) const = 0;
+
+
+    // I-O
+
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/evapotranspirationHeatTransferModel/evapotranspirationHeatTransferModelNew.C

@@ -0,0 +1,62 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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 "evapotranspirationHeatTransferModel.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::evapotranspirationHeatTransferModel>
+Foam::evapotranspirationHeatTransferModel::New
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+{
+    word modelType(dict.get<word>("model"));
+
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+
+    if (!cstrIter.found())
+    {
+        FatalIOErrorInLookup
+        (
+            dict,
+            "model",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
+    }
+
+    return autoPtr<evapotranspirationHeatTransferModel>
+    (
+        cstrIter()(dict, mesh)
+    );
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/grass/grass.C

@@ -0,0 +1,418 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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 "grass.H"
+#include "basicThermo.H"
+#include "fluidThermo.H"
+#include "turbulentTransportModel.H"
+#include "turbulentFluidThermoModel.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace evapotranspirationHeatTransferModels
+{
+    defineTypeNameAndDebug(grass, 0);
+    addToRunTimeSelectionTable
+    (
+        evapotranspirationHeatTransferModel,
+        grass,
+        dictionary
+    );
+}
+}
+
+
+const Foam::Enum
+<
+    Foam::evapotranspirationHeatTransferModels::grass::soilHeatFluxType
+>
+Foam::evapotranspirationHeatTransferModels::grass::soilHeatFluxTypeNames
+({
+    {
+        soilHeatFluxType::PROPORTIONAL_TO_SOLAR_RADIATION,
+        "proportionalToSolarRadiation"
+    },
+    { soilHeatFluxType::BOUNDARY, "boundary" }
+});
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::evapotranspirationHeatTransferModels::grass::E(const labelList& cells)
+const
+{
+    const scalar q = this->q();
+    const scalar G = this->G(cells);
+    const scalar Delta = this->Delta();
+    const scalar D = this->D();
+    const scalar ra = this->ra();
+    const scalar rs = this->rs();
+    const scalar gamma = this->gamma();
+
+    // (BSG:Eq. 11)
+    const scalar E =
+        (Delta*(q - G) + rho_*Cp_*D/ra)/(Delta + gamma*(1 + rs/ra));
+
+    return tmp<scalarField>::New(cells.size(), E);
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::Delta() const
+{
+    // (BSG:Eq. 15), (APR:Eq. 13) - note 0.6108 -> 610.8 due to kPa -> Pa
+    const scalar Ta = Tref_ + 237.3;
+
+    return 4098*(610.8*exp(17.27*Tref_/Ta))/sqr(Ta);
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::D() const
+{
+    // (BSG:Eq. 16)
+    const scalar RHref = RHptr_->value(mesh().time().timeOutputValue());
+
+    return (1 - RHref)*pSat();
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::pSat() const
+{
+    // (BSG:Eq. 17; Arden Buck equation)
+    const scalar p1 = (1.0007 + 3.46e-8*pAtm_)*611.21;
+
+    return p1*exp((18.678 - Tref_/234.5)*Tref_/(Tref_ + 257.14));
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::gamma() const
+{
+    // (APR:Eq. 8)
+    return Cp_*pAtm_/(epsilon_*lambda());
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::lambda() const
+{
+    // (BSG:Eq. 12)
+    const scalar T1 = Tref_ + 273.15;
+    const scalar T2 = Tref_ + 239.24;
+
+    return 1.91846e6*sqr(T1/T2);
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::ra() const
+{
+    // (APR:Eq. 4), (BSG:Eq. 14)
+    const scalar log1 = log((zRefU_ - d_*h_)/(zom_*h_));
+    const scalar log2 = log((zRefH_ - d_*h_)/(zoh_*h_));
+
+    return log1*log2/(sqr(kappa_)*uRef_);
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::rs() const
+{
+    // (BSG:Eq. 13), (APR:Eq. 5; reduced form in p. 4-5)
+    return ri_/(scalar(12)*h_);
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::S
+(
+    const labelList& cells
+) const
+{
+    // Mark fvOption cells within mesh
+    bitSet isZoneCell(mesh().nCells());
+    isZoneCell.set(cells);
+
+    scalar S = 0;
+
+    // Select cells next to specified patches
+    // Sum patch area that is covered by fvOption cells
+    for (const label patchi : patchSet_)
+    {
+        const scalarField& s = mesh().magSf().boundaryField()[patchi];
+
+        const polyPatch& pp = mesh().boundaryMesh()[patchi];
+        const labelList& faceCells = pp.faceCells();
+
+        forAll(faceCells, i)
+        {
+            const bool isCovered = isZoneCell[faceCells[i]];
+
+            if (isCovered)
+            {
+                S += s[i];
+            }
+        }
+    }
+    reduce(S, sumOp<scalar>());
+
+    if (mag(S) < SMALL)
+    {
+        FatalErrorInFunction
+            << "Area coverage of grass cannot be zero. "
+            << "Check whether cellZone has any boundary faces."
+            << exit(FatalError);
+    }
+
+    return S;
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::G
+(
+    const labelList& cells
+) const
+{
+    switch (soilHeatFluxMethod_)
+    {
+        case soilHeatFluxType::PROPORTIONAL_TO_SOLAR_RADIATION:
+        {
+            return Csoil_*q();
+        }
+
+        case soilHeatFluxType::BOUNDARY:
+        {
+            // Retrieve heat flux through patches
+            tmp<FieldField<Field, scalar>> tqBf = this->qBf();
+            const FieldField<Field, scalar>& qBf = tqBf.cref();
+
+            // Mark fvOption cells within mesh
+            bitSet isZoneCell(mesh().nCells());
+            isZoneCell.set(cells);
+
+            scalar G = 0;
+
+            for (const label patchi : patchSet_)
+            {
+                const scalarField& s = mesh().magSf().boundaryField()[patchi];
+
+                const scalarField& qfld = qBf[patchi];
+
+                const polyPatch& pp = mesh().boundaryMesh()[patchi];
+                const labelList& faceCells = pp.faceCells();
+
+                forAll(faceCells, i)
+                {
+                    const bool isCovered = isZoneCell[faceCells[i]];
+
+                    if (isCovered)
+                    {
+                        G += qfld[i]*s[i];
+                    }
+                }
+            }
+            reduce(G, sumOp<scalar>());
+
+            return G/area_;
+        }
+    }
+
+    return -1;
+}
+
+
+Foam::tmp<Foam::FieldField<Foam::Field, Foam::scalar>>
+Foam::evapotranspirationHeatTransferModels::grass::qBf() const
+{
+    const auto& T = mesh().lookupObject<volScalarField>(TName_);
+    const volScalarField::Boundary& Tbf = T.boundaryField();
+
+    auto tq = tmp<FieldField<Field, scalar>>::New(Tbf.size());
+    auto& q = tq.ref();
+
+    forAll(q, patchi)
+    {
+        q.set(patchi, new Field<scalar>(Tbf[patchi].size(), Zero));
+    }
+
+    typedef compressible::turbulenceModel cmpTurbModel;
+
+    if (mesh().foundObject<cmpTurbModel>(cmpTurbModel::propertiesName))
+    {
+        const auto& turb =
+            mesh().lookupObject<cmpTurbModel>(cmpTurbModel::propertiesName);
+
+        // Note: calling he(p,T) instead of he()
+        const volScalarField he(turb.transport().he(turb.transport().p(), T));
+        const volScalarField::Boundary& hebf = he.boundaryField();
+
+        const volScalarField alphaEff(turb.alphaEff());
+        const volScalarField::Boundary& alphaEffbf = alphaEff.boundaryField();
+
+        for (const label patchi : patchSet_)
+        {
+            q[patchi] = alphaEffbf[patchi]*hebf[patchi].snGrad();
+        }
+    }
+    else if (mesh().foundObject<fluidThermo>(fluidThermo::dictName))
+    {
+        const auto& thermo =
+            mesh().lookupObject<fluidThermo>(fluidThermo::dictName);
+
+        // Note: calling he(p,T) instead of he()
+        const volScalarField he(thermo.he(thermo.p(), T));
+        const volScalarField::Boundary& hebf = he.boundaryField();
+
+        const volScalarField& alpha(thermo.alpha());
+        const volScalarField::Boundary& alphabf = alpha.boundaryField();
+
+        for (const label patchi : patchSet_)
+        {
+            q[patchi] = alphabf[patchi]*hebf[patchi].snGrad();
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unable to find a valid thermo model to evaluate q. " << nl
+            << "Database contents are: " << mesh().objectRegistry::sortedToc()
+            << exit(FatalError);
+    }
+
+    // No radiative heat flux contribution is present
+
+    return tq;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::evapotranspirationHeatTransferModels::grass::grass
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    evapotranspirationHeatTransferModel(dict, mesh),
+    soilHeatFluxMethod_
+    (
+        soilHeatFluxTypeNames.getOrDefault
+        (
+            "soilHeatFluxMethod",
+            dict,
+            soilHeatFluxType::PROPORTIONAL_TO_SOLAR_RADIATION
+        )
+    ),
+    Tptr_(nullptr),
+    RHptr_(nullptr),
+    area_(-1),
+    Csoil_(),
+    rho_(),
+    Cp_(),
+    epsilon_(),
+    h_(),
+    kappa_(),
+    uRef_(),
+    zRefU_(),
+    zRefH_(),
+    zom_(),
+    zoh_(),
+    d_(),
+    ri_(),
+    pAtm_(),
+    Tref_(0),
+    TName_(),
+    patchSet_()
+{
+    Info<< "    Activating evapotranspiration heat transfer model: "
+        << typeName << endl;
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::evapotranspirationHeatTransferModels::grass::Q
+(
+    const labelList& cells
+) const
+{
+    if (area_ == -1)
+    {
+        area_ = S(cells);
+    }
+
+    Tref_ = Tptr_->value(mesh().time().timeOutputValue());
+
+    return E(cells)*area_;
+}
+
+
+bool Foam::evapotranspirationHeatTransferModels::grass::read
+(
+    const dictionary& dict
+)
+{
+    if (!evapotranspirationHeatTransferModel::read(dict))
+    {
+        return false;
+    }
+
+    Tptr_.reset(Function1<scalar>::New("Tref", dict, &mesh()));
+    RHptr_.reset(Function1<scalar>::New("RHref", dict, &mesh()));
+
+    area_ = -1;
+
+    const auto range = scalarMinMax::ge(SMALL);
+
+    Csoil_ = dict.getCheckOrDefault<scalar>("Csoil", 0.1, range);
+    rho_ = dict.getCheckOrDefault<scalar>("rho", 1.225, range);
+    Cp_ = dict.getCheckOrDefault<scalar>("Cp", 1013.0, range);
+    epsilon_ = dict.getCheckOrDefault<scalar>("epsilon", 0.622, range);
+    h_ = dict.getCheckOrDefault<scalar>("h", 0.1, range);
+    kappa_ = dict.getCheckOrDefault<scalar>("kappa", 0.41, range);
+    uRef_ = dict.getCheckOrDefault<scalar>("uRef", 2, range);
+    zRefU_ = dict.getCheckOrDefault<scalar>("zRefU", 10, range);
+    zRefH_ = dict.getCheckOrDefault<scalar>("zRefH", 10, range);
+    zom_ = dict.getCheckOrDefault<scalar>("zom", 0.123, range);
+    zoh_ = dict.getCheckOrDefault<scalar>("zoh", 0.0123, range);
+    d_ = dict.getCheckOrDefault<scalar>("d", 2.0/3.0, range);
+    ri_ = dict.getCheckOrDefault<scalar>("ri", 100, range);
+    pAtm_ = dict.getCheckOrDefault<scalar>("pAtm", 101.325, range);
+    TName_ = dict.getOrDefault<word>("T", "T");
+
+    if (soilHeatFluxMethod_ == soilHeatFluxType::BOUNDARY)
+    {
+        patchSet_ =
+            mesh().boundaryMesh().patchSet(dict.get<wordRes>("patches"));
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/grass/grass.H

@@ -0,0 +1,342 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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/>.
+
+Class
+    Foam::evapotranspirationHeatTransferModels::grass
+
+Description
+    Applies sources on temperature (\c T - incompressible) or energy
+    (\c h/e - compressible) equation to incorporate evapotranspiration
+    heat-transfer effects from the specified grass canopy.
+    The model is based on Pemnan-Monteith Equation model.
+
+    Sources applied to either of the below, if exist:
+    \verbatim
+      T         | Temperature                                [K]
+      e         | Internal energy                            [m^2/s^2]
+      h         | Enthalphy                                  [m^2/s^2]
+    \endverbatim
+
+    Required fields:
+    \verbatim
+      T         | Temperature                                [K]
+      e         | Internal energy                            [m^2/s^2]
+      h         | Enthalphy                                  [m^2/s^2]
+      LAD       | Leaf area density                          [m^2/m^3]
+    \endverbatim
+
+    References:
+    \verbatim
+        Governing equations (tag:BSG):
+            Brozovsky, J., Simonsen, A., & Gaitani, N. (2021).
+            Validation of a CFD model for the evaluation of urban
+            microclimate at high latitudes: A case study in Trondheim, Norway.
+            Building and Environment, 205, 108175.
+            DOI:10.1016/j.buildenv.2021.108175
+
+        Governing equations (tag:APR):
+            Allen, R. G., Pereira, L. S., Raes, D., & Smith, M. (1998).
+            Crop evapotranspiration-Guidelines for computing crop
+            water requirements-FAO Irrigation and drainage paper 56.
+            Fao, Rome, 300(9), D05109.
+    \endverbatim
+
+Usage
+    Example by using \c constant/fvOptions:
+    \verbatim
+    evapotranspirationHeatTransfer1
+    {
+        // Inherited entries
+        ...
+
+        // Mandatory entries
+        model               grass;
+        Tref                <Function1<scalar>>;
+        RHref               <Function1<scalar>>;
+
+        // Conditional entries
+
+            // when soilHeatFluxMethod == boundary
+            patches             <wordRes>;
+
+        // Optional entries
+        soilHeatFluxMethod  <word>;
+        Csoil               <scalar>;
+        rho                 <scalar>;
+        Cp                  <scalar>;
+        epsilon             <scalar>;
+        h                   <scalar>;
+        kappa               <scalar>;
+        uRef                <scalar>;
+        zRefU               <scalar>;
+        zRefH               <scalar>;
+        zom                 <scalar>;
+        zoh                 <scalar>;
+        d                   <scalar>;
+        ri                  <scalar>;
+        pAtm                <scalar>;
+        T                   <word>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property | Description                         | Type | Reqd | Deflt
+      model    | Model name: grass                   | word | yes  | -
+      Tref     | Reference weather station air temperature [Celsius] <!--
+               -->                     | Function1\<scalar\> | yes | -
+      RHref    | Reference weather station relative humidity [%]     <!--
+               -->                     | Function1\<scalar\> | yes | -
+      patches  | Names of ground patches           | wordRes | yes | -
+      soilHeatFluxMethod | Method to calculate soil heat flux - see below <!--
+               -->                   | word | no | -
+      Csoil    | Proportionality constant of soil heat-flux wrt <!--
+               --> solar radiation                   | scalar | no | 0.1
+      rho      | Mean air density at constant pressure [kg/m^3]      <!--
+               -->                                   | scalar | no | 1.225
+      Cp       | Specific heat at constant pressure [J/kg/C]         <!--
+               -->                                   | scalar | no | 1013.0
+      epsilon  | Molecular-weight ratio of water vapour/dry air [-]  <!--
+               -->                                  | scalar | no  | 0.622
+      h        | Height of grass layer [m]          | scalar | no  | 0.1
+      kappa    | Von Karman constant [-]            | scalar | no  | 0.41
+      uRef     | Reference velocity magnitude [m/s] | scalar | no  | 2.0
+      zRefU    | Height of wind speed measurements [m] | scalar | no | 10.0
+      zRefH    | Height of humidity measurements [m]   | scalar | no | 10.0
+      zom      | Roughness length governing momentum transfer coefficient <!--
+               -->                                  [-] | scalar | no | 0.123
+      zoh      | Roughness length governing transfer of heat and vapour <!--
+               --> coefficient [-]                      | scalar | no | 0.0123
+      d        | Zero plane displacement height coefficient [-]         <!--
+               -->                                      | scalar | no | 0.666
+      ri       | Bulk stomata resistance of leaf [s/m] | scalar | no | 100.0
+      pAtm     | Atmospheric pressure [kPa]            | scalar | no | 101.325
+      T        | Name of temperature field             | word   | no | T
+    \endtable
+
+    Options for the \c soilHeatFluxMethod entry:
+    \verbatim
+      proportionalToSolarRadiation | Estimate from solar load
+      boundary                     | Obtain soil heat flux from boundary
+    \endverbatim
+
+    The inherited entries are elaborated in:
+      - \link evapotranspirationHeatTransfer.H \endlink
+      - \link evapotranspirationHeatTransferModel.H \endlink
+
+SourceFiles
+    grass.C
+    grassTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_evapotranspirationHeatTransferModels_grass_H
+#define Foam_evapotranspirationHeatTransferModels_grass_H
+
+#include "evapotranspirationHeatTransferModel.H"
+#include "Function1.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace evapotranspirationHeatTransferModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class grass Declaration
+\*---------------------------------------------------------------------------*/
+
+class grass
+:
+    public evapotranspirationHeatTransferModel
+{
+    // Private Enumerations
+
+        //- Options for the soil heat flux
+        enum soilHeatFluxType : char
+        {
+            PROPORTIONAL_TO_SOLAR_RADIATION = 0, //!< "Estimate from solar load"
+            BOUNDARY,            //!< "Obtain soil heat flux from boundary"
+        };
+
+        //- Names for soilHeatFluxType
+        static const Enum<soilHeatFluxType> soilHeatFluxTypeNames;
+
+
+    // Private Data
+
+        //- Soil heat-flux calculation method
+        enum soilHeatFluxType soilHeatFluxMethod_;
+
+        //- Reference weather station air temperature [Celsius]
+        autoPtr<Function1<scalar>> Tptr_;
+
+        //- Reference weather station relative humidity [%]
+        autoPtr<Function1<scalar>> RHptr_;
+
+        //- Area coverage of grass [m^2]
+        mutable scalar area_;
+
+        //- Proportionality constant of soil heat-flux wrt solar radiation [-]
+        scalar Csoil_;
+
+        //- Mean air density at constant pressure [kg/m^3]
+        scalar rho_;
+
+        //- Specific heat at constant pressure [J/kg/C]
+        scalar Cp_;
+
+        //- Molecular-weight ratio of water vapour/dry air [-]
+        scalar epsilon_;
+
+        //- Height of grass layer [m]
+        scalar h_;
+
+        //- Von Karman constant [-]
+        scalar kappa_;
+
+        //- Reference velocity magnitude [m/s]
+        scalar uRef_;
+
+        //- Height of wind speed measurements [m]
+        scalar zRefU_;
+
+        //- Height of humidity measurements [m]
+        scalar zRefH_;
+
+        //- Roughness length governing momentum transfer coefficient [-]
+        scalar zom_;
+
+        //- Roughness length governing transfer of heat and vapour coeff [-]
+        scalar zoh_;
+
+        //- Zero plane displacement height coefficient [-]
+        scalar d_;
+
+        //- Bulk stomata resistance of leaf [s/m]
+        scalar ri_;
+
+        //- Atmospheric pressure [kPa]
+        scalar pAtm_;
+
+        //- Cached reference weather station air temperature [Celsius]
+        mutable scalar Tref_;
+
+        //- Name of temperature field
+        word TName_;
+
+        //- List of patches to calculate boundary heat flux
+        labelHashSet patchSet_;
+
+
+    // Private Member Functions
+
+        //- Return evapotranspiration heat flux through grass layer [W/m^2]
+        tmp<scalarField> E(const labelList& cells) const;
+
+        //- Return slope of the relation between
+        //- vapour pressure-temperature [Pa/K]
+        scalar Delta() const;
+
+        //- Return atmospheric vapour pressure deficit [Pa]
+        scalar D() const;
+
+        //- Return saturated vapour pressure over water [Pa]
+        scalar pSat() const;
+
+        //- Return psychrometric constant [kPa/K]
+        scalar gamma() const;
+
+        //- Return specific latent heat of vaporisation [MJ/kg]
+        scalar lambda() const;
+
+        //- Return bulk aerodynamic resistance [Pa]
+        scalar ra() const;
+
+        //- Return bulk surface resistance [Pa]
+        scalar rs() const;
+
+        //- Return area coverage of grass [m^2]
+        scalar S(const labelList& cells) const;
+
+        //- Return area-averaged heat flux through ground [W/m^2]
+        scalar G(const labelList& cells) const;
+
+        //- Return heat-flux boundary fields
+        tmp<FieldField<Field, scalar>> qBf() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("grass");
+
+
+    // Constructors
+
+        //- Construct from components
+        grass
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- No copy construct
+        grass(const grass&) = delete;
+
+        //- No copy assignment
+        void operator=(const grass&) = delete;
+
+
+    //- Destructor
+    virtual ~grass() = default;
+
+
+    // Member Functions
+
+    // Evaluation
+
+        //- Return heat-transfer rate for speficied cells [J/s]
+        virtual tmp<scalarField> Q(const labelList& cells) const;
+
+
+    // I-O
+
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace evapotranspirationHeatTransferModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/tree/tree.C

@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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 "tree.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace evapotranspirationHeatTransferModels
+{
+    defineTypeNameAndDebug(tree, 0);
+    addToRunTimeSelectionTable
+    (
+        evapotranspirationHeatTransferModel,
+        tree,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::tree::Et() const
+{
+    return a_*q() + b_;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::evapotranspirationHeatTransferModels::tree::leafArea
+(
+    const labelList& cells
+) const
+{
+    const volScalarField& LAD = getOrReadField(LADname_);
+
+    const scalarField& V = mesh().V();
+
+    auto tleafArea = tmp<scalarField>::New(cells.size(), Zero);
+    auto& leafArea = tleafArea.ref();
+
+    forAll(cells, i)
+    {
+        const label celli = cells[i];
+
+        leafArea[i] = LAD[celli]*V[celli];
+    }
+
+    return tleafArea;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::evapotranspirationHeatTransferModels::tree::tree
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    evapotranspirationHeatTransferModel(dict, mesh),
+    a_(),
+    b_(),
+    lambda_(),
+    LADname_()
+{
+    Info<< "    Activating evapotranspiration heat transfer model: "
+        << typeName << endl;
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::evapotranspirationHeatTransferModels::tree::Q
+(
+    const labelList& cells
+) const
+{
+    // Convert units from [MJ g/hr] to [J kg/s]
+    static const scalar unitConverter = scalar(1000)/scalar(3600);
+
+    return unitConverter*lambda_*Et()*leafArea(cells);
+}
+
+
+bool Foam::evapotranspirationHeatTransferModels::tree::read
+(
+    const dictionary& dict
+)
+{
+    if (!evapotranspirationHeatTransferModel::read(dict))
+    {
+        return false;
+    }
+
+    const auto range = scalarMinMax::ge(SMALL);
+
+    a_ = dict.getOrDefault<scalar>("a", 0.3622);
+    b_ = dict.getOrDefault<scalar>("b", 60.758);
+    lambda_ = dict.getCheckOrDefault<scalar>("lambda", 2.44, range);
+    LADname_ = dict.getOrDefault<word>("LAD", "LAD");
+
+    (void) getOrReadField(LADname_);
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/tree/tree.H

@@ -0,0 +1,219 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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/>.
+
+Class
+    Foam::evapotranspirationHeatTransferModels::tree
+
+Description
+    Applies sources on temperature (\c T - incompressible) or energy
+    (\c h/e - compressible) equation to incorporate evapotranspiration
+    heat-transfer effects from the specified tree canopy.
+
+    The heat-transfer is calculated based on the following relations.
+    The heat transfer is given by the first equation and the
+    evapotranspiration is calculated linearly proportional to absorbed
+    solar radiation (direct and diffusive) with the empirical relation
+    shown in the second equation below.
+
+    \f[
+        Q = \lambda E_t \frac{1000}{3600} \int_{cellZone} LAD \, dV
+    \f]
+
+    where
+    \vartable
+      Q         | Heat transfer from tree crown                 [W]
+      \lambda   | Specific latent heat of vaporisation          [MJ/kg]
+      E_t       | Area-averaged evapotranspiration              [g/m^2/hr]
+      LAD       | Leaf area density                             [m^2/m^3]
+      V         | Volume                                        [m^3]
+    \endvartable
+
+    with
+
+    \f[
+        E_t = a R_n + b
+    \f]
+
+    where
+    \vartable
+      R_n    | Heat flux through tree crown                    [W/m^2]
+      a      | Linear-regression coefficient - slope parameter [g/W/hr]
+      b      | Linear-regression coefficient - intercept parameter [g/m^2/hr]
+    \endvartable
+
+    Sources applied to either of the below, if exist:
+    \verbatim
+      T         | Temperature                                [K]
+      e         | Internal energy                            [m^2/s^2]
+      h         | Enthalphy                                  [m^2/s^2]
+    \endverbatim
+
+    Required fields:
+    \verbatim
+      T         | Temperature                                [K]
+      e         | Internal energy                            [m^2/s^2]
+      h         | Enthalphy                                  [m^2/s^2]
+      LAD       | Leaf area density                          [m^2/m^3]
+    \endverbatim
+
+Usage
+    Example by using \c constant/fvOptions:
+    \verbatim
+    evapotranspirationHeatTransfer1
+    {
+        // Inherited entries
+        ...
+
+        // Mandatory entries
+        model               tree;
+
+        // Optional entries
+        a                   <scalar>;
+        b                   <scalar>;
+        lambda              <scalar>;
+        LAD                 <word>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property | Description                         | Type | Reqd | Deflt
+      model    | Model name: tree                    | word | yes  | -
+      a | Linear-regression coefficient - slope parameter [g/W/hr] | scalar <!--
+        -->                                                   | no | 0.3622
+      b | Linear-regression coefficient - intercept parameter [g/m^2/hr] <!--
+        -->                                         | scalar  | no | 60.758
+      lambda   | Specific latent heat of vaporisation [MJ/kg] | scalar <!--
+        -->                                                   | no | 2.44
+      LAD      | Name of leaf-area density field       | word | no | LAD
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link evapotranspirationHeatTransfer.H \endlink
+      - \link evapotranspirationHeatTransferModel.H \endlink
+
+SourceFiles
+    tree.C
+    treeTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_evapotranspirationHeatTransferModels_tree_H
+#define Foam_evapotranspirationHeatTransferModels_tree_H
+
+#include "evapotranspirationHeatTransferModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace evapotranspirationHeatTransferModels
+{
+
+/*---------------------------------------------------------------------------*\
+                            Class tree Declaration
+\*---------------------------------------------------------------------------*/
+
+class tree
+:
+    public evapotranspirationHeatTransferModel
+{
+    // Private Data
+
+        //- Linear-regression coefficient - slope parameter [g/W/hr]
+        scalar a_;
+
+        //- Linear-regression coefficient - intercept parameter [g/m^2/hr]
+        scalar b_;
+
+        //- Specific latent heat of vaporisation [MJ/kg]
+        scalar lambda_;
+
+        //- Name of leaf-area density field
+        word LADname_;
+
+
+    // Private Member Functions
+
+        //- Return area-averaged transpiration [g/m^2/hr]
+        //  Calculated from empirical regressions with heat flux
+        scalar Et() const;
+
+        //- Return volume weighted leaf area density (LAD) [m^2]
+        //  LAD [m^2/m^3]: one-sided area of leaves per unit volume
+        tmp<scalarField> leafArea(const labelList& cells) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("tree");
+
+
+    // Constructors
+
+        //- Construct from components
+        tree
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- No copy construct
+        tree(const tree&) = delete;
+
+        //- No copy assignment
+        void operator=(const tree&) = delete;
+
+
+    //- Destructor
+    virtual ~tree() = default;
+
+
+    // Member Functions
+
+    // Evaluation
+
+        //- Return heat-transfer rate for speficied cells [J/s]
+        virtual tmp<scalarField> Q(const labelList& cells) const;
+
+
+    // I-O
+
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace evapotranspirationHeatTransferModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulence.C

@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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 "treeTurbulence.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(treeTurbulence, 0);
+    addToRunTimeSelectionTable(option, treeTurbulence, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+Foam::volScalarField& Foam::fv::treeTurbulence::getOrReadField
+(
+    const word& fieldName
+) const
+{
+    auto* ptr = mesh_.getObjectPtr<volScalarField>(fieldName);
+
+    if (!ptr)
+    {
+        ptr = new volScalarField
+        (
+            IOobject
+            (
+                fieldName,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh_
+        );
+        mesh_.objectRegistry::store(ptr);
+    }
+
+    return *ptr;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::treeTurbulence::treeTurbulence
+(
+    const word& sourceName,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    fv::cellSetOption(sourceName, modelType, dict, mesh),
+    isEpsilon_(true),
+    betaP_(),
+    betaD_(),
+    Ceps1_(),
+    Ceps2_(),
+    betaStar_(),
+    CdName_(),
+    LADname_()
+{
+    read(dict);
+
+    const auto* turbPtr = mesh_.findObject<turbulenceModel>
+    (
+        turbulenceModel::propertiesName
+    );
+
+    if (!turbPtr)
+    {
+        FatalErrorInFunction
+            << "Unable to find a turbulence model."
+            << abort(FatalError);
+    }
+
+    fieldNames_.resize(2);
+
+    tmp<volScalarField> tepsilon = turbPtr->epsilon();
+    tmp<volScalarField> tomega = turbPtr->omega();
+
+    if (!tepsilon.isTmp())
+    {
+        fieldNames_[0] = tepsilon().name();
+    }
+    else if (!tomega.isTmp())
+    {
+        isEpsilon_ = false;
+        fieldNames_[0] = tomega().name();
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unable to find epsilon or omega field." << nl
+            << abort(FatalError);
+    }
+
+    fieldNames_[1] = turbPtr->k()().name();
+
+    fv::option::resetApplied();
+
+    Log << "    Applying treeTurbulence to: "
+        << fieldNames_[0] << " and " << fieldNames_[1]
+        << endl;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fv::treeTurbulence::addSup
+(
+    fvScalarMatrix& eqn,
+    const label fieldi
+)
+{
+    if (fieldi == 1)
+    {
+        kSource(geometricOneField(), geometricOneField(), eqn, fieldi);
+        return;
+    }
+
+    if (isEpsilon_)
+    {
+        epsilonSource(geometricOneField(), geometricOneField(), eqn, fieldi);
+    }
+    else
+    {
+        omegaSource(geometricOneField(), geometricOneField(), eqn, fieldi);
+    }
+}
+
+
+void Foam::fv::treeTurbulence::addSup
+(
+    const volScalarField& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+)
+{
+    if (fieldi == 1)
+    {
+        kSource(geometricOneField(), rho, eqn, fieldi);
+        return;
+    }
+
+    if (isEpsilon_)
+    {
+        epsilonSource(geometricOneField(), rho, eqn, fieldi);
+    }
+    else
+    {
+        omegaSource(geometricOneField(), rho, eqn, fieldi);
+    }
+}
+
+
+void Foam::fv::treeTurbulence::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+)
+{
+    if (fieldi == 1)
+    {
+        kSource(alpha, rho, eqn, fieldi);
+        return;
+    }
+
+    if (isEpsilon_)
+    {
+        epsilonSource(alpha, rho, eqn, fieldi);
+    }
+    else
+    {
+        omegaSource(alpha, rho, eqn, fieldi);
+    }
+}
+
+
+bool Foam::fv::treeTurbulence::read(const dictionary& dict)
+{
+    if (!fv::cellSetOption::read(dict))
+    {
+        return false;
+    }
+
+    betaP_ = dict.getOrDefault<scalar>("betaP", 1.0);
+    betaD_ = dict.getOrDefault<scalar>("betaD", 4.0);
+    Ceps1_ = dict.getOrDefault<scalar>("Ceps1", 0.9);
+    Ceps2_ = dict.getOrDefault<scalar>("Ceps2", 0.9);
+    betaStar_ = dict.getOrDefault<scalar>("betaStar", 0.09);
+    CdName_ = dict.getOrDefault<word>("Cd", "Cd");
+    LADname_ = dict.getOrDefault<word>("LAD", "LAD");
+
+    (void) getOrReadField(CdName_);
+    (void) getOrReadField(LADname_);
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulence.H

@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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/>.
+
+Class
+    Foam::fv::treeTurbulence
+
+Group
+    grpFvOptionsSources
+
+Description
+    Applies sources on \c k and either \c epsilon or \c omega
+    to incorporate effects of trees on atmospheric boundary layer.
+
+    Corrections applied to:
+    \verbatim
+      k         | Turbulent kinetic energy                  [m^2/s^2]
+    \endverbatim
+
+    Corrections applied to either of the below, if exist:
+    \verbatim
+      epsilon   | Turbulent kinetic energy dissipation rate [m^2/s^3]
+      omega     | Specific dissipation rate                 [1/s]
+    \endverbatim
+
+    Required fields:
+    \verbatim
+      k         | Turbulent kinetic energy                  [m^2/s^2]
+      Cd        | Canopy drag coefficient                   [-]
+      LAD       | Leaf area density                         [m^2/m^3]
+      epsilon   | Turbulent kinetic energy dissipation rate [m^2/s^3]
+      omega     | Specific dissipation rate                 [1/s]
+    \endverbatim
+
+    References:
+    \verbatim
+        Governing equations (tag:BSG):
+            Brozovsky, J., Simonsen, A., & Gaitani, N. (2021).
+            Validation of a CFD model for the evaluation of urban microclimate
+            at high latitudes: A case study in Trondheim, Norway.
+            Building and Environment, 205, 108175.
+            DOI:10.1016/j.buildenv.2021.108175
+    \endverbatim
+
+Usage
+    Example by using \c constant/fvOptions:
+    \verbatim
+    treeTurbulence1
+    {
+        // Mandatory entries
+        type         treeTurbulence;
+
+        // Optional entries
+        Cd           <word>;
+        LAD          <word>;
+        betaP        <scalar>;
+        betaD        <scalar>;
+        Ceps1        <scalar>;
+        Ceps2        <scalar>;
+        betaStar     <scalar>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property | Description                         | Type | Reqd | Deflt
+      type     | Type name: treeTurbulence           | word | yes  | -
+      Cd       | Name of operand canopy drag coefficient field | word | no | Cd
+      LAD      | Name of operand leaf area density field | word | no | LAD
+      betaP    | Model coefficient                 | scalar | no   | 1.0
+      betaD    | Model coefficient                 | scalar | no   | 4.0
+      Ceps1    | Model coefficient                 | scalar | no   | 0.9
+      Ceps2    | Model coefficient                 | scalar | no   | 0.9
+      betaStar | Model coefficient                 | scalar | no   | 0.09
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link fvOption.H \endlink
+      - \link cellSetOption.H \endlink
+
+SourceFiles
+    treeTurbulence.C
+    treeTurbulenceTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fv_treeTurbulence_H
+#define fv_treeTurbulence_H
+
+#include "cellSetOption.H"
+#include "turbulentTransportModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class treeTurbulence Declaration
+\*---------------------------------------------------------------------------*/
+
+class treeTurbulence
+:
+    public fv::cellSetOption
+{
+    // Private Data
+
+        //- Internal flag to determine the working field is epsilon or omega
+        bool isEpsilon_;
+
+        // Model coefficients
+
+            scalar betaP_;
+            scalar betaD_;
+            scalar Ceps1_;
+            scalar Ceps2_;
+            scalar betaStar_;
+
+        //- Name of operand canopy drag coefficient field
+        word CdName_;
+
+        //- Name of operand leaf area density field
+        word LADname_;
+
+
+     // Private Member Functions
+
+        //- Return requested field from the object registry
+        //- or read+register the field to the object registry
+        volScalarField& getOrReadField(const word& fieldName) const;
+
+        //- Apply sources to turbulent kinetic energy
+        template<class AlphaFieldType, class RhoFieldType>
+        void kSource
+        (
+            const AlphaFieldType& alpha,
+            const RhoFieldType& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        ) const;
+
+        //- Apply sources to turbulent kinetic energy dissipation rate
+        template<class AlphaFieldType, class RhoFieldType>
+        void epsilonSource
+        (
+            const AlphaFieldType& alpha,
+            const RhoFieldType& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        ) const;
+
+        //- Apply sources to specific dissipation rate
+        template<class AlphaFieldType, class RhoFieldType>
+        void omegaSource
+        (
+            const AlphaFieldType& alpha,
+            const RhoFieldType& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("treeTurbulence");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        treeTurbulence
+        (
+            const word& sourceName,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- No copy construct
+        treeTurbulence(const treeTurbulence&) = delete;
+
+        //- No copy assignment
+        void operator=(const treeTurbulence&) = delete;
+
+
+    // Member Functions
+
+        //- Add explicit contribution to epsilon or omega equation
+        //- for incompressible flow computations
+         virtual void addSup
+        (
+            fvScalarMatrix& eqn,
+            const label fieldi
+        );
+
+        //- Add explicit contribution to epsilon or omega equation
+        //- for compressible flow computations
+        virtual void addSup
+        (
+            const volScalarField& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        );
+
+        //- Add explicit contribution to epsilon or omega equation
+        //- for multiphase flow computations
+        virtual void addSup
+        (
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        );
+
+        //- Read source dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "treeTurbulenceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulenceTemplates.C

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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 "treeTurbulence.H"
+#include "volFields.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * *  Member Functions * * * * * * * * * * * * * //
+
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::fv::treeTurbulence::kSource
+(
+    const AlphaFieldType& alpha,
+    const RhoFieldType& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+) const
+{
+    const auto* turbPtr = mesh_.findObject<turbulenceModel>
+    (
+        turbulenceModel::propertiesName
+    );
+
+    const volScalarField& k = turbPtr->k();
+    const volVectorField& U = turbPtr->U();
+    const volScalarField magU(mag(U));
+    const volScalarField& Cd = getOrReadField(CdName_);
+    const volScalarField& LAD = getOrReadField(LADname_);
+
+    // (BSG:Eq. 8)
+    eqn +=
+        alpha*rho*LAD*Cd*betaP_*pow3(magU)
+      - fvm::Sp(alpha*rho*LAD*Cd*betaD_*magU, k);
+}
+
+
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::fv::treeTurbulence::epsilonSource
+(
+    const AlphaFieldType& alpha,
+    const RhoFieldType& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+) const
+{
+    const auto* turbPtr = mesh_.findObject<turbulenceModel>
+    (
+        turbulenceModel::propertiesName
+    );
+
+    const volScalarField& epsilon = turbPtr->epsilon();
+    const volScalarField& k = turbPtr->k();
+    const volVectorField& U = turbPtr->U();
+    const volScalarField magU(mag(U));
+    const volScalarField& Cd = getOrReadField(CdName_);
+    const volScalarField& LAD = getOrReadField(LADname_);
+
+    // (BSG:Eq. 9)
+    eqn +=
+        fvm::Sp
+        (
+            alpha*rho*LAD*Cd*(Ceps1_*betaP_/k*pow3(magU) - Ceps2_*betaD_*magU),
+            epsilon
+        );
+}
+
+
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::fv::treeTurbulence::omegaSource
+(
+    const AlphaFieldType& alpha,
+    const RhoFieldType& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+) const
+{
+    const auto* turbPtr = mesh_.findObject<turbulenceModel>
+    (
+        turbulenceModel::propertiesName
+    );
+
+    const volScalarField& omega = turbPtr->omega();
+    const volScalarField& k = turbPtr->k();
+    const volVectorField& U = turbPtr->U();
+    const volScalarField magU(mag(U));
+    const volScalarField& Cd = getOrReadField(CdName_);
+    const volScalarField& LAD = getOrReadField(LADname_);
+
+    // (derived from BSG:Eq. 9 by assuming epsilon = betaStar*omega*k)
+    eqn +=
+        fvm::Sp
+        (
+            alpha*rho*LAD*Cd*betaStar_
+           *(Ceps1_*betaP_*pow3(magU) - Ceps2_*betaD_*k*magU),
+            omega
+        );
+}
+
+
+// ************************************************************************* //