Merge branch 'development' into refactor_prob_init

ewquon · Sep 8, 2023 · 04ae88c · 04ae88c
2 parents 5bb0c60 + 6f05140
commit 04ae88c
Show file tree

Hide file tree

Showing 6 changed files with 486 additions and 13 deletions.
diff --git a/Docs/sphinx_doc/BoundaryConditions.rst b/Docs/sphinx_doc/BoundaryConditions.rst
@@ -311,6 +311,14 @@ In ERF, when the MOST boundary condition is applied, velocity and temperature in
 
 MOST Inputs
 ~~~~~~~~~~~~~~~~~~~
+To evaluate the fluxes with MOST, the surface rougness parameter :math:`z_{0}` must be specified. This quantity may be considered a constant or may be parameterized through the friction velocity :math:`u_{\star}`. ERF supports three methods for parameterizing the surface roughness: ``constant``, ``charnock``, and ``modified_charnock``. The latter two methods parameterize :math:`z_{0} = f(u_{\star})` and are described in `Jimenez & Dudhia, American Meteorological Society, 2018 <https://doi.org/10.1175/JAMC-D-17-0137.1>`_. The rougness calculation method may be specified with
+
+::
+
+   erf.most.roughness_type    = STRING    #Z_0 type (constant, charnock, modified_charnock)
+
+If the ``charnock`` method is employed, the :math:`a` constant may be specified with ``erf.most.charnock_constant`` (defaults to 0.0185). If the ``modified_charnock`` method is employed, the depth :math:`d` may be specified with ``erf.most.modified_charnock_depth`` (defaults to 30 m).
+
 When computing an average :math:`\overline{\phi}` for the MOST boundary, where :math:`\phi` denotes a generic variable, ERF supports a variety of approaches. Specifically, ``planar averages`` and ``local region averages`` may be computed with or without ``time averaging``. With each averaging methodology, the query point :math:`z` may be determined from the following procedures: specified vertical distance :math:`z_{ref}` from the bottom surface, specified :math:`k_{index}`, or (when employing terrain-fit coordinates) specified normal vector length :math:`z_{ref}`. The available inputs to the MOST boundary and their associated data types are
 
 ::

diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H
@@ -62,6 +62,21 @@ public:
             }
         }
 
+        // Specify how to compute the flux
+        std::string rough_string{"constant"};
+        pp.query("most.roughness_type", rough_string);
+        if (rough_string == "constant") {
+            rough_type = RoughCalcType::CONSTANT;
+        } else if (rough_string == "charnock") {
+            rough_type = RoughCalcType::CHARNOCK;
+            pp.query("most.charnock_constant",cnk_a);
+        } else if (rough_string == "modified_charnock") {
+            rough_type = RoughCalcType::MODIFIED_CHARNOCK;
+            pp.query("most.modified_charnock_depth",depth);
+        } else {
+            amrex::Abort("Undefined MOST roughness type!");
+        }
+
         int nlevs = m_geom.size();
         z_0.resize(nlevs);
         u_star.resize(nlevs);
@@ -164,13 +179,22 @@ public:
         SURFACE_TEMPERATURE ///< Surface temperature specified
     };
 
-    ThetaCalcType theta_type;
+    enum struct RoughCalcType {
+        CONSTANT = 0,      ///< Constant z0
+        CHARNOCK,
+        MODIFIED_CHARNOCK
+    };
+
     FluxCalcType  flux_type;
+    ThetaCalcType theta_type;
+    RoughCalcType rough_type;
 
 private:
     amrex::Real z0_const;
     amrex::Real surf_temp;
     amrex::Real surf_temp_flux{0};
+    amrex::Real cnk_a{0.0185};
+    amrex::Real depth{30.0};
     amrex::Vector<amrex::Geometry>  m_geom;
     amrex::Vector<amrex::FArrayBox> z_0;
 

diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp
@@ -19,16 +19,40 @@ ABLMost::update_fluxes (int lev,
     // Iterate the fluxes if moeng type
     if (flux_type == FluxCalcType::MOENG) {
         if (theta_type == ThetaCalcType::HEAT_FLUX) {
-            surface_flux most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux);
-            compute_fluxes(lev, max_iters, most_flux);
+            if (rough_type == RoughCalcType::CONSTANT) {
+                surface_flux most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux);
+                compute_fluxes(lev, max_iters, most_flux);
+            } else if (rough_type == RoughCalcType::CHARNOCK) {
+                surface_flux_charnock most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux, cnk_a);
+                compute_fluxes(lev, max_iters, most_flux);
+            } else {
+                surface_flux_mod_charnock most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux, depth);
+                compute_fluxes(lev, max_iters, most_flux);
+            }
         } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) {
-            surface_temp most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux);
-            compute_fluxes(lev, max_iters, most_flux);
+            if (rough_type == RoughCalcType::CONSTANT) {
+                surface_temp most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux);
+                compute_fluxes(lev, max_iters, most_flux);
+            } else if (rough_type == RoughCalcType::CHARNOCK) {
+                surface_temp_charnock most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux, cnk_a);
+                compute_fluxes(lev, max_iters, most_flux);
+            } else {
+                surface_temp_mod_charnock most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux, depth);
+                compute_fluxes(lev, max_iters, most_flux);
+            }
         } else {
-            adiabatic most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux);
-            compute_fluxes(lev, max_iters, most_flux);
-        }
-    }
+            if (rough_type == RoughCalcType::CONSTANT) {
+                adiabatic most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux);
+                compute_fluxes(lev, max_iters, most_flux);
+            } else if (rough_type == RoughCalcType::CHARNOCK) {
+                adiabatic_charnock most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux, cnk_a);
+                compute_fluxes(lev, max_iters, most_flux);
+            } else {
+                adiabatic_mod_charnock most_flux(m_ma.get_zref(),surf_temp, surf_temp_flux, depth);
+                compute_fluxes(lev, max_iters, most_flux);
+            }
+        } // theta flux
+    } // Moeng flux
 }
 
 

diff --git a/Source/BoundaryConditions/MOSTAverage.H b/Source/BoundaryConditions/MOSTAverage.H
@@ -164,7 +164,7 @@ protected:
     int m_navg{4};                                                   // 4 averages for U/V/T/Umag
     int m_maxlev{0};                                                 // Total number of levels
     int m_policy{0};                                                 // Policy for type of averaging
-    amrex::Real m_zref{-1.0};                                        // Height above surface for MOST BC
+    amrex::Real m_zref{10.0};                                        // Height above surface for MOST BC
     std::string m_pp_prefix {"erf"};                                 // ParmParse prefix
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> m_x_pos;         // Ptr to 2D mf to hold x position (maxlev)
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> m_y_pos;         // Ptr to 2D mf to hold y position (maxlev)

diff --git a/Source/BoundaryConditions/MOSTAverage.cpp b/Source/BoundaryConditions/MOSTAverage.cpp
@@ -246,11 +246,15 @@ MOSTAverage::set_k_indices_N()
     if (read_z) {
         for (int lev(0); lev < m_maxlev; lev++) {
             Real m_zlo = m_geom[lev].ProbLo(2);
+            Real m_zhi = m_geom[lev].ProbHi(2);
             Real m_dz  = m_geom[lev].CellSize(2);
 
             AMREX_ASSERT_WITH_MESSAGE(m_zref >= m_zlo + 0.5 * m_dz,
                                       "Query point must be past first z-cell!");
 
+            AMREX_ASSERT_WITH_MESSAGE(m_zref <= m_zhi - 0.5 * m_dz,
+                                      "Query point must be below the last z-cell!");
+
             int lk = static_cast<int>(floor((m_zref - m_zlo) / m_dz - 0.5));
 
             m_zref = (lk + 0.5) * m_dz + m_zlo;