Skip to content

Commit

Permalink
Merge branch 'development' into JAMESinputs
Browse files Browse the repository at this point in the history
  • Loading branch information
AMLattanzi authored Dec 9, 2024
2 parents e2455c1 + 7849906 commit 31b30b1
Show file tree
Hide file tree
Showing 13 changed files with 64 additions and 33 deletions.
8 changes: 6 additions & 2 deletions Source/BoundaryConditions/ERF_ABLMost.H
Original file line number Diff line number Diff line change
Expand Up @@ -368,15 +368,19 @@ public:
update_pblh (const int& lev,
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars,
amrex::MultiFab* z_phys_cc,
const int RhoQv_comp, const int RhoQr_comp);
const int RhoQv_comp,
const int RhoQc_comp,
const int RhoQr_comp);

template <typename PBLHeightEstimator>
void
compute_pblh (const int& lev,
amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars,
amrex::MultiFab* z_phys_cc,
const PBLHeightEstimator& est,
const int RhoQv_comp, const int RhoQr_comp);
const int RhoQv_comp,
const int RhoQc_comp,
const int RhoQr_comp);

void
read_custom_roughness (const int& lev,
Expand Down
12 changes: 8 additions & 4 deletions Source/BoundaryConditions/ERF_ABLMost.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -561,11 +561,13 @@ void
ABLMost::update_pblh (const int& lev,
Vector<Vector<MultiFab>>& vars,
MultiFab* z_phys_cc,
int RhoQv_comp, int RhoQr_comp)
int RhoQv_comp,
int RhoQc_comp,
int RhoQr_comp)
{
if (pblh_type == PBLHeightCalcType::MYNN25) {
MYNNPBLH estimator;
compute_pblh(lev, vars, z_phys_cc, estimator, RhoQv_comp, RhoQr_comp);
compute_pblh(lev, vars, z_phys_cc, estimator, RhoQv_comp, RhoQc_comp, RhoQr_comp);
} else if (pblh_type == PBLHeightCalcType::YSU) {
amrex::Error("YSU PBLH calc not implemented yet");
}
Expand All @@ -577,11 +579,13 @@ ABLMost::compute_pblh (const int& lev,
Vector<Vector<MultiFab>>& vars,
MultiFab* z_phys_cc,
const PBLHeightEstimator& est,
int RhoQv_comp, int RhoQr_comp)
int RhoQv_comp,
int RhoQc_comp,
int RhoQr_comp)
{
est.compute_pblh(m_geom[lev],z_phys_cc, pblh[lev].get(),
vars[lev][Vars::cons],m_lmask_lev[lev][0],
RhoQv_comp, RhoQr_comp);
RhoQv_comp, RhoQc_comp, RhoQr_comp);
}

void
Expand Down
6 changes: 6 additions & 0 deletions Source/DataStructs/ERF_DataStruct.H
Original file line number Diff line number Diff line change
Expand Up @@ -107,17 +107,22 @@ struct SolverChoice {
pp.query_enum_case_insensitive("moisture_model",moisture_type);
if (moisture_type == MoistureType::SAM) {
RhoQv_comp = RhoQ1_comp;
RhoQc_comp = RhoQ2_comp;
RhoQr_comp = RhoQ4_comp;
} else if (moisture_type == MoistureType::SAM_NoIce) {
RhoQv_comp = RhoQ1_comp;
RhoQc_comp = RhoQ2_comp;
RhoQr_comp = RhoQ4_comp;
} else if (moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
RhoQv_comp = RhoQ1_comp;
RhoQc_comp = RhoQ2_comp;
} else if (moisture_type == MoistureType::Kessler) {
RhoQv_comp = RhoQ1_comp;
RhoQc_comp = RhoQ2_comp;
RhoQr_comp = RhoQ3_comp;
} else if (moisture_type == MoistureType::Kessler_NoRain) {
RhoQv_comp = RhoQ1_comp;
RhoQc_comp = RhoQ2_comp;
}

// TODO: should we set default for dry??
Expand Down Expand Up @@ -637,6 +642,7 @@ struct SolverChoice {
bool do_precip {true};
bool use_moist_background {false};
int RhoQv_comp {-1};
int RhoQc_comp {-1};

// This component will be model-dependent:
// if a model with no rain, this will stay -1
Expand Down
3 changes: 2 additions & 1 deletion Source/Diffusion/ERF_ComputeTurbulentViscosity.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -489,7 +489,8 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel ,
ComputeDiffusivityMYNN25(xvel, yvel, cons_in, eddyViscosity,
geom, turbChoice, most, use_moisture,
level, bc_ptr, vert_only, z_phys_nd,
solverChoice.RhoQv_comp, solverChoice.RhoQr_comp);
solverChoice.RhoQv_comp, solverChoice.RhoQc_comp,
solverChoice.RhoQr_comp);
} else if (turbChoice.pbl_type == PBLType::YSU) {
ComputeDiffusivityYSU(xvel, yvel, cons_in, eddyViscosity,
geom, turbChoice, most, use_moisture,
Expand Down
3 changes: 2 additions & 1 deletion Source/Diffusion/ERF_DiffusionSrcForState_N.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -713,6 +713,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1;

const int rhoqv_comp = solverChoice.RhoQv_comp;
const int rhoqc_comp = solverChoice.RhoQc_comp;
const int rhoqr_comp = solverChoice.RhoQr_comp;

ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
Expand All @@ -728,7 +729,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
mu_turb,cellSizeInv,domain,
pbl_mynn_B1_l,tm_arr(i,j,0),
rhoqv_comp, rhoqr_comp,
rhoqv_comp, rhoqc_comp, rhoqr_comp,
c_ext_dir_on_zlo, c_ext_dir_on_zhi,
u_ext_dir_on_zlo, u_ext_dir_on_zhi,
v_ext_dir_on_zlo, v_ext_dir_on_zhi,
Expand Down
3 changes: 2 additions & 1 deletion Source/Diffusion/ERF_DiffusionSrcForState_T.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -824,6 +824,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
auto pbl_mynn_B1_l = turbChoice.pbl_mynn.B1;

const int rhoqv_comp = solverChoice.RhoQv_comp;
const int rhoqc_comp = solverChoice.RhoQc_comp;
const int rhoqr_comp = solverChoice.RhoQr_comp;

ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
Expand All @@ -839,7 +840,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain,
cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
mu_turb,cellSizeInv,domain,
pbl_mynn_B1_l,tm_arr(i,j,0),
rhoqv_comp, rhoqr_comp,
rhoqv_comp, rhoqc_comp, rhoqr_comp,
c_ext_dir_on_zlo, c_ext_dir_on_zhi,
u_ext_dir_on_zlo, u_ext_dir_on_zhi,
v_ext_dir_on_zlo, v_ext_dir_on_zhi,
Expand Down
4 changes: 3 additions & 1 deletion Source/ERF.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1061,7 +1061,9 @@ ERF::InitData_post ()
// we don't want to call update_fluxes multiple times because
// it will change u* and theta* from their previous values
m_most->update_pblh(lev, vars_new, z_phys_cc[lev].get(),
solverChoice.RhoQv_comp, solverChoice.RhoQr_comp);
solverChoice.RhoQv_comp,
solverChoice.RhoQc_comp,
solverChoice.RhoQr_comp);
m_most->update_fluxes(lev, time);
}
}
Expand Down
3 changes: 2 additions & 1 deletion Source/PBL/ERF_ComputeDiffusivityMYNN25.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel,
bool /*vert_only*/,
const std::unique_ptr<MultiFab>& z_phys_nd,
const int RhoQv_comp,
const int RhoQc_comp,
const int RhoQr_comp)
{
const bool use_terrain = (z_phys_nd != nullptr);
Expand Down Expand Up @@ -140,7 +141,7 @@ ComputeDiffusivityMYNN25 (const MultiFab& xvel,
u_ext_dir_on_zlo, u_ext_dir_on_zhi,
v_ext_dir_on_zlo, v_ext_dir_on_zhi,
dthetadz, dudz, dvdz,
RhoQv_comp, RhoQr_comp, use_most);
RhoQv_comp, RhoQc_comp, RhoQr_comp, use_most);

// Spatially varying MOST
Real theta0 = tm_arr(i,j,0);
Expand Down
3 changes: 2 additions & 1 deletion Source/PBL/ERF_ComputeDiffusivityYSU.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -204,13 +204,14 @@ ComputeDiffusivityYSU (const MultiFab& xvel,
constexpr Real prandtl_max = 4.0;
Real dthetadz, dudz, dvdz;
const int RhoQv_comp = -1;
const int RhoQc_comp = -1;
const int RhoQr_comp = -1;
ComputeVerticalDerivativesPBL(i, j, k,
uvel, vvel, cell_data, izmin, izmax, 1.0/dz_terrain,
c_ext_dir_on_zlo, c_ext_dir_on_zhi,
u_ext_dir_on_zlo, u_ext_dir_on_zhi,
v_ext_dir_on_zlo, v_ext_dir_on_zhi,
dthetadz, dudz, dvdz, RhoQv_comp, RhoQr_comp);
dthetadz, dudz, dvdz, RhoQv_comp, RhoQc_comp, RhoQr_comp);
const Real shear_squared = dudz*dudz + dvdz*dvdz + 1.0e-9; // 1.0e-9 from WRF to avoid divide by zero
const Real theta = cell_data(i,j,k,RhoTheta_comp) / cell_data(i,j,k,Rho_comp);
Real richardson = CONST_GRAV / theta * dthetadz / shear_squared;
Expand Down
15 changes: 8 additions & 7 deletions Source/PBL/ERF_PBLHeight.H
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ struct MYNNPBLH {
const amrex::MultiFab& cons,
const amrex::iMultiFab* lmask,
const int RhoQv_comp,
const int RhoQc_comp,
const int RhoQr_comp) const
{
#if 0
Expand All @@ -35,7 +36,7 @@ struct MYNNPBLH {
auto thetav_min = amrex::ReduceToPlane<amrex::ReduceOpMin,amrex::Real>(dir, bxlow, cons,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) -> amrex::Real
{
return Thetav(i,j,k,cons_arrs[box_no],RhoQv_comp,RhoQr_comp);
return Thetav(i,j,k,cons_arrs[box_no],RhoQv_comp,RhoQc_comp,RhoQr_comp);
});
#endif

Expand Down Expand Up @@ -88,7 +89,7 @@ struct MYNNPBLH {
int jj = amrex::max(amrex::min(j,jmax),jmin);

if (zphys_arr(ii,jj,k) < thetamin_height) {
amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQr_comp);
amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
if (min_thv_arr(i,j,0) > thv) min_thv_arr(i,j,0) = thv;
}
});
Expand All @@ -104,8 +105,8 @@ struct MYNNPBLH {
//
// Find PBL height based on thetav increase (best for CBLs)
//
amrex::Real thv = Thetav(i,j,k ,cons_arr,RhoQv_comp,RhoQr_comp);
amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQr_comp);
amrex::Real thv = Thetav(i,j,k ,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);

int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
if (is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_land)
Expand Down Expand Up @@ -164,7 +165,7 @@ struct MYNNPBLH {

ParallelFor(gtbxlow, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQr_comp);
amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
if (min_thv_arr(i,j,0) > thv) min_thv_arr(i,j,0) = thv;
});

Expand All @@ -176,8 +177,8 @@ struct MYNNPBLH {
//
// Find PBL height based on thetav increase (best for CBLs)
//
amrex::Real thv = Thetav(i,j,k ,cons_arr,RhoQv_comp,RhoQr_comp);
amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQr_comp);
amrex::Real thv = Thetav(i,j,k ,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);

int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
if (is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_land)
Expand Down
28 changes: 16 additions & 12 deletions Source/PBL/ERF_PBLModels.H
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,9 @@ ComputeDiffusivityMYNN25 (const amrex::MultiFab& xvel,
const amrex::BCRec* bc_ptr,
bool /*vert_only*/,
const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
const int RhoQv_comp, const int RhoQr_comp);
const int RhoQv_comp,
const int RhoQc_comp,
const int RhoQr_comp);

/**
* Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
Expand Down Expand Up @@ -94,23 +96,24 @@ ComputeVerticalDerivativesPBL (int i, int j, int k,
amrex::Real& dudz,
amrex::Real& dvdz,
const int RhoQv_comp,
const int RhoQc_comp,
const int RhoQr_comp,
const bool use_most=true)
{
if ( k==izmax && c_ext_dir_on_zhi ) {
dthetadz = (1.0/3.0)*( -Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQr_comp)
- 3.0 * Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQr_comp)
+ 4.0 * Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv;
dthetadz = (1.0/3.0)*( -Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
- 3.0 * Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
+ 4.0 * Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv;
} else if ( k==izmin && c_ext_dir_on_zlo ) {
dthetadz = (1.0/3.0)*( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp)
+ 3.0 * Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQr_comp)
- 4.0 * Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv;
dthetadz = (1.0/3.0)*( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
+ 3.0 * Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
- 4.0 * Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv;
} else if ( k==izmin && use_most ) {
dthetadz = ( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp)
- Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv;
dthetadz = ( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
- Thetav(i,j,k ,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv;
} else {
dthetadz = 0.5*( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQr_comp)
- Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQr_comp) )*dz_inv;
dthetadz = 0.5*( Thetav(i,j,k+1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp)
- Thetav(i,j,k-1,cell_data,RhoQv_comp,RhoQc_comp,RhoQr_comp) )*dz_inv;
}

if ( k==izmax && u_ext_dir_on_zhi ) {
Expand Down Expand Up @@ -165,6 +168,7 @@ ComputeQKESourceTerms (int i, int j, int k,
amrex::Real pbl_mynn_B1_l,
const amrex::Real theta_mean,
const int RhoQv_comp,
const int RhoQc_comp,
const int RhoQr_comp,
bool c_ext_dir_on_zlo,
bool c_ext_dir_on_zhi,
Expand Down Expand Up @@ -195,7 +199,7 @@ ComputeQKESourceTerms (int i, int j, int k,
u_ext_dir_on_zlo, u_ext_dir_on_zhi,
v_ext_dir_on_zlo, v_ext_dir_on_zhi,
dthetadz, dudz, dvdz,
RhoQv_comp, RhoQr_comp, use_most);
RhoQv_comp, RhoQc_comp, RhoQr_comp, use_most);

// Note: Transport terms due to turbulence and pressure are included when
// DiffusionSrcForState_* is called from ERF_slow_rhs_post.
Expand Down
4 changes: 3 additions & 1 deletion Source/TimeIntegration/ERF_Advance.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,9 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/)
// Reassign the field ptrs for MAC avg computation.
m_most->update_mac_ptrs(lev, vars_old, Theta_prim, Qv_prim, Qr_prim);
m_most->update_pblh(lev, vars_old, z_phys_cc[lev].get(),
solverChoice.RhoQv_comp, solverChoice.RhoQr_comp);
solverChoice.RhoQv_comp,
solverChoice.RhoQc_comp,
solverChoice.RhoQr_comp);
m_most->update_fluxes(lev, time);
}
}
Expand Down
5 changes: 4 additions & 1 deletion Source/Utils/ERF_Thetav.H
Original file line number Diff line number Diff line change
Expand Up @@ -10,15 +10,18 @@ amrex::Real
Thetav (int i, int j, int k,
const amrex::Array4<const amrex::Real>& cell_data,
const int RhoQv_comp,
const int RhoQc_comp,
const int RhoQr_comp)
{
amrex::Real thetav = cell_data(i,j,k,RhoTheta_comp) / cell_data(i,j,k,Rho_comp);

if (RhoQr_comp > 0) {
thetav *= (1.0 + 0.61 * cell_data(i,j,k,RhoQv_comp) / cell_data(i,j,k,Rho_comp)
- cell_data(i,j,k,RhoQc_comp) / cell_data(i,j,k,Rho_comp)
- cell_data(i,j,k,RhoQr_comp) / cell_data(i,j,k,Rho_comp));
} else if (RhoQv_comp > 0) {
thetav *= (1.0 + 0.61 * cell_data(i,j,k,RhoQ1_comp) / cell_data(i,j,k,Rho_comp));
thetav *= (1.0 + 0.61 * cell_data(i,j,k,RhoQ1_comp) / cell_data(i,j,k,Rho_comp)
- cell_data(i,j,k,RhoQc_comp) / cell_data(i,j,k,Rho_comp));
}

return thetav;
Expand Down

0 comments on commit 31b30b1

Please sign in to comment.