Skip to content

Commit

Permalink
Implement ABLMost::get_q_star and buoyancy flux calculation in MYNN
Browse files Browse the repository at this point in the history 
Verified dry solution does not change
  • Loading branch information
@ewquon
ewquon committed Aug 14, 2024
1 parent 5c8c978 commit 1eb6435
Show file tree
Hide file tree
Showing 5 changed files with 31 additions and 9 deletions.
3 changes: 3 additions & 0 deletions Source/BoundaryConditions/ABLMost.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(); }

Expand Down
5 changes: 4 additions & 1 deletion Source/Diffusion/ComputeTurbulentViscosity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel,
const Geometry& geom,
const TurbChoice& turbChoice,
std::unique_ptr<ABLMost>& most,
bool use_moisture,
int level,
const BCRec* bc_ptr,
bool /*vert_only*/,
Expand Down Expand Up @@ -482,6 +483,7 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel ,
const TurbChoice& turbChoice, const Real const_grav,
std::unique_ptr<ABLMost>& most,
const bool& exp_most,
const bool& use_moisture,
int level,
const BCRec* bc_ptr,
bool vert_only)
Expand Down Expand Up @@ -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);
}
}
1 change: 1 addition & 0 deletions Source/Diffusion/EddyViscosity.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab&
const TurbChoice& turbChoice, const amrex::Real const_grav,
std::unique_ptr<ABLMost>& most,
const bool& exp_most,
const bool& use_moisture,
int level,
const amrex::BCRec* bc_ptr,
bool vert_only = false);
Expand Down
28 changes: 21 additions & 7 deletions Source/Diffusion/PBLModels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel,
const Geometry& geom,
const TurbChoice& turbChoice,
std::unique_ptr<ABLMost>& most,
bool use_moisture,
int level,
const BCRec* bc_ptr,
bool /*vert_only*/,
Expand Down Expand Up @@ -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<Real>{};

const Array4<Real const> z_nd_arr = use_terrain ? z_phys_nd->array(mfi) : Array4<Real>{};

Expand All @@ -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<Real>::max();
}
Expand Down
3 changes: 2 additions & 1 deletion Source/TimeIntegration/ERF_advance_dycore.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
}

// ***********************************************************************************************
Expand Down

0 comments on commit 1eb6435

Please sign in to comment.