Forest canopy

CMake/BuildERFExe.cmake

@@ -169,6 +169,7 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/SourceTerms/ERF_make_sources.cpp
        ${SRC_DIR}/SourceTerms/ERF_moist_set_rhs.cpp
        ${SRC_DIR}/SourceTerms/ERF_NumericalDiffusion.cpp
+       ${SRC_DIR}/SourceTerms/ERF_ForestDrag.cpp
        ${SRC_DIR}/TimeIntegration/ERF_ComputeTimestep.cpp
        ${SRC_DIR}/TimeIntegration/ERF_Advance.cpp
        ${SRC_DIR}/TimeIntegration/ERF_TimeStep.cpp

Docs/sphinx_doc/Inputs.rst

@@ -1096,6 +1096,9 @@ List of Parameters
 | **erf.input_sounding_file**         | Name(s) of the         | String(s)         | input_sounding_file |
 |                                     | input sounding file(s) |                   |                     |
 +-------------------------------------+------------------------+-------------------+---------------------+
+| **erf.forest_file**                 | Name(s) of the         | String            | None                |
+|                                     | canopy forest file     |                   |                     |
++-------------------------------------+------------------------+-------------------+---------------------+
 | **erf.input_sounding_time**         | Time(s) of the         | Real(s)           | false               |
 |                                     | input sounding file(s) |                   |                     |
 +-------------------------------------+------------------------+-------------------+---------------------+

Docs/sphinx_doc/index.rst

@@ -57,6 +57,7 @@ In addition to this documentation, there is API documentation for ERF generated
    theory/Forcings.rst
    Particles.rst
    theory/WindFarmModels.rst
+   theory/Forest.rst
    theory/UnitsAndConstants.rst
 
 .. toctree::

Docs/sphinx_doc/theory/Forest.rst

@@ -0,0 +1,31 @@
+Forest Model
+--------------
+
+The forest model provides an option to include the drag from forested regions to be included in the momentum equation. The
+drag force is calculated as follows:
+
+.. math::
+
+   F_i= - C_d L(x,y,z) U_i | U_i |
+
+
+Here :math:`C_d` is the coefficient of drag for the forested region and :math:`L(x,y,z)` is the leaf area density (LAD) for the
+forested region. A three-dimensional model for the LAD is usually unavailable and is also cumbersome to use if there are thousands
+of trees. Two different models are available as an alternative:
+
+.. math::
+   L=\frac{LAI}{h}
+
+.. math::
+   L(z)=L_m \left(\frac{h - z_m}{h - z}\right)^n  exp\left[n \left(1 -\frac{h - z_m}{h - z}\right )\right]
+
+Here :math:`LAI` is the leaf area index and is available from measurements, :math:`h` is the height of the tree, :math:`z_m` is the location
+of the maximum LAD, :math:`L_m` is the maximum value of LAD at :math:`z_m` and :math:`n` is a model constant with values  6 (below :math:`z_m`) and 0.5
+(above :math:`z_m`), respectively. :math:`L_m` is computed by integrating the following equation (see `Lalic and Mihailovic (2004)
+<https://doi.org/10.1175/1520-0450(2004)043<0641:AERDLD>2.0.CO;2>`_):
+
+.. math::
+   LAI = \int_{0}^{h} L(z) dz
+
+The simplified model with uniform LAD is recommended for forested regions with no knowledge of the individual trees. LAI values can be used from
+climate model look-up tables for different regions around the world if no local remote sensing data is available.

Exec/ABL/erf_forest_def

@@ -0,0 +1,4 @@
+1  1024 512 45 200  0.2 6 0.8
+1  1024 1024 35 200  0.2 6 0.8
+1  1024 1224 75 200  0.2 6 0.8 
+2  1024 1524 120 200 0.2 10 0.8

Exec/ABL/inputs_canopy

@@ -0,0 +1,64 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 50
+
+amrex.fpe_trap_invalid = 1
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_extent =  2048     2048    1024
+amr.n_cell           =    64       64      32
+
+geometry.is_periodic = 1 1 0
+
+zlo.type = "SlipWall"
+zhi.type = "SlipWall"
+
+# TIME STEP CONTROL
+erf.fixed_dt = 2.0  # fixed time step depending on grid resolution
+erf.cfl      = 0.5
+
+# DIAGNOSTICS & VERBOSITY
+erf.sum_interval   = 1       # timesteps between computing mass
+erf.v              = 1       # verbosity in ERF.cpp
+amr.v              = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+
+# CHECKPOINT FILES
+erf.check_file      = chk        # root name of checkpoint file
+erf.check_int       = 50         # number of timesteps between checkpoints
+
+# PLOTFILES
+erf.plot_file_1     = plt        # prefix of plotfile name
+erf.plot_int_1      = 50        # number of timesteps between plotfiles
+erf.plot_vars_1     = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta
+
+# SOLVER CHOICE
+erf.alpha_T = 0.0
+erf.alpha_C = 1.0
+erf.use_gravity = false
+
+erf.molec_diff_type = "None"
+erf.les_type        = "Smagorinsky"
+erf.Cs              = 0.1
+
+erf.init_type = "uniform"
+
+#erf.forest_file = erf_forest_def
+
+# PROBLEM PARAMETERS
+prob.rho_0 = 1.0
+prob.A_0 = 1.0
+
+prob.U_0 = 10.0
+prob.V_0 = 0.0
+prob.W_0 = 0.0
+prob.T_0 = 300.0
+
+# Higher values of perturbations lead to instability
+# Instability seems to be coming from BC
+prob.U_0_Pert_Mag = 0.08
+prob.V_0_Pert_Mag = 0.08 #
+prob.W_0_Pert_Mag = 0.0

Source/BoundaryConditions/ERF_ABLMost.H

@@ -431,15 +431,15 @@ public:
     get_t_surf (const int& lev) { return t_surf[lev].get(); }
 
     amrex::Real
-    get_zref () {return m_ma.get_zref();}
+    get_zref () { return m_ma.get_zref(); }
 
     const amrex::FArrayBox*
-    get_z0 (const int& lev) {return &z_0[lev];}
+    get_z0 (const int& lev) { return &z_0[lev]; }
 
     bool have_variable_sea_roughness() { return m_var_z0; }
 
     const amrex::iMultiFab*
-    get_lmask(const int& lev) {return m_lmask_lev[lev][0];}
+    get_lmask(const int& lev) { return m_lmask_lev[lev][0]; }
 
     int
     lmask_min_reduce (amrex::iMultiFab& lmask, const int& nghost)

Source/DataStructs/ERF_DataStruct.H

@@ -644,5 +644,8 @@ struct SolverChoice {
     amrex::Real turb_disk_angle = -1.0;
     amrex::Real windfarm_x_shift = -1.0;
     amrex::Real windfarm_y_shift = -1.0;
+
+    // Flag for valid canopy model
+    bool do_forest {false};
 };
 #endif

Source/ERF.H

@@ -45,6 +45,7 @@
 #include <ERF_PhysBCFunct.H>
 #include <ERF_FillPatcher.H>
 #include <ERF_SampleData.H>
+#include <ERF_ForestDrag.H>
 
 #ifdef ERF_USE_PARTICLES
 #include "ERF_ParticleData.H"
@@ -1123,6 +1124,7 @@ private:
     std::unique_ptr<WriteBndryPlanes> m_w2d  = nullptr;
     std::unique_ptr<ReadBndryPlanes>  m_r2d  = nullptr;
     std::unique_ptr<ABLMost>          m_most = nullptr;
+    amrex::Vector<std::unique_ptr<ForestDrag>> m_forest;
 
     //
     // Holds info for dynamically generated tagging criteria

Source/ERF.cpp

@@ -30,7 +30,7 @@ Real ERF::change_max    =  1.1;
 int  ERF::fixed_mri_dt_ratio = 0;
 
 // Dictate verbosity in screen output
-int ERF::verbose        = 0;
+int  ERF::verbose       = 0;
 int  ERF::mg_verbose    = 0;
 bool ERF::use_fft       = false;
 
@@ -129,6 +129,10 @@ ERF::ERF_shared ()
     lsm_data.resize(nlevs_max);
     lsm_flux.resize(nlevs_max);
 
+    // NOTE: size canopy model before readparams (if file exists, we construct)
+    m_forest.resize(nlevs_max);
+    for (int lev = 0; lev < max_level; ++lev) { m_forest[lev] = nullptr; }
+
     ReadParameters();
     initializeMicrophysics(nlevs_max);
 
@@ -300,6 +304,7 @@ ERF::ERF_shared ()
         Hwave_onegrid[lev] = nullptr;
         Lwave_onegrid[lev] = nullptr;
     }
+
     // Theta prim for MOST
     Theta_prim.resize(nlevs_max);
 
@@ -1588,6 +1593,15 @@ ERF::ReadParameters ()
         pp.query("cf_width", cf_width);
         pp.query("cf_set_width", cf_set_width);
 
+        // Query the canopy model file name
+        std::string forestfile;
+        solverChoice.do_forest = pp.query("forest_file", forestfile);
+        if (solverChoice.do_forest) {
+            for (int lev = 0; lev <= max_level; ++lev) {
+                m_forest[lev] = std::make_unique<ForestDrag>(forestfile);
+            }
+        }
+
         // AmrMesh iterate on grids?
         bool iterate(true);
         pp_amr.query("iterate_grids",iterate);

Source/ERF_make_new_level.cpp

@@ -131,6 +131,11 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in,
         }
     }
 
+    // ********************************************************************************************
+    // Build the data structures for canopy model (depends upon z_phys)
+    // ********************************************************************************************
+    if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
+
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
@@ -205,6 +210,11 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
         }
     }
 
+    // ********************************************************************************************
+    // Build the data structures for canopy model (depends upon z_phys)
+    // ********************************************************************************************
+    if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
+
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
@@ -329,6 +339,11 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
         }
     }
 
+    // ********************************************************************************************
+    // Build the data structures for canopy model (depends upon z_phys)
+    // ********************************************************************************************
+    if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
+
     // *****************************************************************************************************
     // Create the physbcs objects (after initializing the terrain but before calling FillCoarsePatch
     // *****************************************************************************************************

Source/SourceTerms/ERF_ForestDrag.H

@@ -0,0 +1,40 @@
+#ifndef ERF_FORESTDRAG_H_
+#define ERF_FORESTDRAG_H_
+
+#include <memory>
+#include <AMReX_MultiFab.H>
+
+/*
+  ForestDrag flow physics adapted from:
+  Lalic & Mihailovic (2004)
+  https://doi.org/10.1175/1520-0450(2004)043<0641:AERDLD>2.0.CO;2
+ */
+class ForestDrag
+{
+public:
+
+    explicit ForestDrag (std::string forestfile);
+
+    ~ForestDrag () = default;
+
+    void
+    define_drag_field (const amrex::BoxArray& ba,
+                       const amrex::DistributionMapping& dm,
+                       amrex::Geometry& geom,
+                       amrex::MultiFab* z_phys_nd);
+
+    amrex::MultiFab*
+    get_drag_field () { return m_forest_drag.get(); }
+
+private:
+    amrex::Vector<amrex::Real> m_type_forest;
+    amrex::Vector<amrex::Real> m_x_forest;
+    amrex::Vector<amrex::Real> m_y_forest;
+    amrex::Vector<amrex::Real> m_height_forest;
+    amrex::Vector<amrex::Real> m_diameter_forest;
+    amrex::Vector<amrex::Real> m_cd_forest;
+    amrex::Vector<amrex::Real> m_lai_forest;
+    amrex::Vector<amrex::Real> m_laimax_forest;
+    std::unique_ptr<amrex::MultiFab> m_forest_drag;
+};
+#endif