From 1eb6435e63784ff75bc9f06f81c6bf907c480ea7 Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Wed, 14 Aug 2024 15:34:17 -0600 Subject: [PATCH] Implement ABLMost::get_q_star and buoyancy flux calculation in MYNN Verified dry solution does not change --- Source/BoundaryConditions/ABLMost.H | 3 ++ .../Diffusion/ComputeTurbulentViscosity.cpp | 5 +++- Source/Diffusion/EddyViscosity.H | 1 + Source/Diffusion/PBLModels.cpp | 28 ++++++++++++++----- Source/TimeIntegration/ERF_advance_dycore.cpp | 3 +- 5 files changed, 31 insertions(+), 9 deletions(-) diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index 3a555bd17..7f99dd44c 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -303,6 +303,9 @@ public: const amrex::MultiFab* get_t_star (const int& lev) { return t_star[lev].get(); } + const amrex::MultiFab* + get_q_star (const int& lev) { return q_star[lev].get(); } + const amrex::MultiFab* get_olen (const int& lev) { return olen[lev].get(); } diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index 1d7f0266c..46e815582 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -16,6 +16,7 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, const Geometry& geom, const TurbChoice& turbChoice, std::unique_ptr& most, + bool use_moisture, int level, const BCRec* bc_ptr, bool /*vert_only*/, @@ -482,6 +483,7 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel , const TurbChoice& turbChoice, const Real const_grav, std::unique_ptr& most, const bool& exp_most, + const bool& use_moisture, int level, const BCRec* bc_ptr, bool vert_only) @@ -521,6 +523,7 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel , if (turbChoice.pbl_type != PBLType::None) { // NOTE: state_new is passed in for Cons_old (due to ptr swap in advance) ComputeTurbulentViscosityPBL(xvel, yvel, cons_in, eddyViscosity, - geom, turbChoice, most, level, bc_ptr, vert_only, z_phys_nd); + geom, turbChoice, most, use_moisture, + level, bc_ptr, vert_only, z_phys_nd); } } diff --git a/Source/Diffusion/EddyViscosity.H b/Source/Diffusion/EddyViscosity.H index 94607eeee..55e7b783d 100644 --- a/Source/Diffusion/EddyViscosity.H +++ b/Source/Diffusion/EddyViscosity.H @@ -20,6 +20,7 @@ ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab& const TurbChoice& turbChoice, const amrex::Real const_grav, std::unique_ptr& most, const bool& exp_most, + const bool& use_moisture, int level, const amrex::BCRec* bc_ptr, bool vert_only = false); diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp index 9dc0676bf..410339e12 100644 --- a/Source/Diffusion/PBLModels.cpp +++ b/Source/Diffusion/PBLModels.cpp @@ -26,6 +26,7 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, const Geometry& geom, const TurbChoice& turbChoice, std::unique_ptr& most, + bool use_moisture, int level, const BCRec* bc_ptr, bool /*vert_only*/, @@ -139,13 +140,17 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, Real d_kappa = KAPPA; Real d_gravity = CONST_GRAV; - const auto& t_mean_mf = most->get_mac_avg(level,4); // This is theta_v - const auto& u_star_mf = most->get_u_star(level); // Use desired level - const auto& t_star_mf = most->get_t_star(level); // Use desired level + const auto& t_mean_mf = most->get_mac_avg(level,4); // theta_v + const auto& q_mean_mf = most->get_mac_avg(level,3); // q_v + const auto& u_star_mf = most->get_u_star(level); + const auto& t_star_mf = most->get_t_star(level); + const auto& q_star_mf = most->get_q_star(level); const auto& tm_arr = t_mean_mf->array(mfi); + const auto& qm_arr = t_mean_mf->array(mfi); const auto& u_star_arr = u_star_mf->array(mfi); const auto& t_star_arr = t_star_mf->array(mfi); + const auto& q_star_arr = (use_moisture) ? q_star_mf->array(mfi) : Array4{}; const Array4 z_nd_arr = use_terrain ? z_phys_nd->array(mfi) : Array4{}; @@ -163,12 +168,21 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel, dthetadz, dudz, dvdz); // Spatially varying MOST - Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0); // TODO: need to account for moisture - Real theta0 = tm_arr(i,j,0); + Real theta0 = tm_arr(i,j,0); + Real qv0 = qm_arr(i,j,0); + Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0); + Real surface_latent_heat{0}; + if (use_moisture) { + // Compute buoyancy flux (Stull Eqn. 4.4.5d) + surface_latent_heat = -u_star_arr(i,j,0) * q_star_arr(i,j,0); + surface_heat_flux *= (1.0 + 0.61*qv0); + surface_heat_flux += 0.61 * theta0 * surface_latent_heat; + } + Real l_obukhov; if (std::abs(surface_heat_flux) > eps) { - l_obukhov = ( theta0 * u_star_arr(i,j,0) * u_star_arr(i,j,0) ) / - ( d_kappa * d_gravity * t_star_arr(i,j,0) ); + l_obukhov = -( theta0 * u_star_arr(i,j,0)*u_star_arr(i,j,0)*u_star_arr(i,j,0) ) + / ( d_kappa * d_gravity * surface_heat_flux ); } else { l_obukhov = std::numeric_limits::max(); } diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 702d9081e..209f9be4c 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -93,6 +93,7 @@ void ERF::advance_dycore(int level, (tc.pbl_type != PBLType::None) ); bool l_use_kturb = ( (tc.les_type != LESType::None) || (tc.pbl_type != PBLType::None) ); + bool l_use_moisture = ( solverChoice.moisture_type != MoistureType::None ); const bool use_most = (m_most != nullptr); const bool exp_most = (solverChoice.use_explicit_most); @@ -219,7 +220,7 @@ void ERF::advance_dycore(int level, *eddyDiffs, *Hfx1, *Hfx2, *Hfx3, *Diss, // to be updated fine_geom, *mapfac_u[level], *mapfac_v[level], z_phys_nd[level], tc, solverChoice.gravity, - m_most, exp_most, level, bc_ptr_d); + m_most, exp_most, l_use_moisture, level, bc_ptr_d); } // ***********************************************************************************************