Skip to content

Commit

Permalink
fix SAM for moist bubble case, add qsat as a derived quantity and pri…
Browse files Browse the repository at this point in the history 
…nt buoyancy_type at the start (#1629)
  • Loading branch information
@asalmgren
asalmgren authored May 25, 2024
1 parent 120f71f commit fb65da6
Show file tree
Hide file tree
Showing 5 changed files with 105 additions and 49 deletions.
2 changes: 2 additions & 0 deletions Source/DataStructs/DataStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -350,6 +350,8 @@ struct SolverChoice {
amrex::Print() << ")" << std::endl;
}

amrex::Print() << "Buoyancy_type : " << buoyancy_type << std::endl;

advChoice.display();
diffChoice.display();
spongeChoice.display();
Expand Down
7 changes: 5 additions & 2 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -799,7 +799,8 @@ private:
const amrex::Vector<std::string> derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar",
"vorticity_x","vorticity_y","vorticity_z",
"magvel", "divU",
"pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "eq_pot_temp", "num_turb",
"pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens",
"eq_pot_temp", "num_turb",
"dpdx", "dpdy", "pres_hse_x", "pres_hse_y",
"z_phys", "detJ" , "mapfac", "lat_m", "lon_m",
// Time averaged velocity
Expand All @@ -811,8 +812,10 @@ private:
// mynn pbl lengthscale
"Lpbl",
// moisture vars
"qt", "qv", "qc", "qi", "qp", "qrain", "qsnow", "qgraup", "rain_accum", "snow_accum", "graup_accum"
"qt", "qv", "qc", "qi", "qp", "qrain", "qsnow", "qgraup", "qsat",
"rain_accum", "snow_accum", "graup_accum"
#ifdef ERF_COMPUTE_ERROR
// error vars
,"xvel_err", "yvel_err", "zvel_err", "pp_err"
#endif
};
Expand Down
48 changes: 44 additions & 4 deletions Source/IO/Plotfile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,9 +51,15 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>

for (int i = 0; i < ncomp_cons; ++i) {
if ( containerHasElement(plot_var_names, cons_names[i]) ) {
tmp_plot_names.push_back(cons_names[i]);
if ( (solverChoice.moisture_type == MoistureType::SAM) || (derived_names[i] != "rhoQ4" &&
derived_names[i] != "rhoQ5" &&
derived_names[i] != "rhoQ6") )
{
tmp_plot_names.push_back(cons_names[i]);
} // moisture_type
}
}

// check for velocity since it's not in cons_names
// if we are asked for any velocity component, we will need them all
if (containerHasElement(plot_var_names, "x_velocity") ||
Expand All @@ -63,13 +69,26 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
tmp_plot_names.push_back("y_velocity");
tmp_plot_names.push_back("z_velocity");
}

//
// If the model we are running doesn't have the variable listed in the inputs file,
// just ignore it rather than aborting
//
for (int i = 0; i < derived_names.size(); ++i) {
if ( containerHasElement(plot_var_names, derived_names[i]) ) {
if (solverChoice.use_terrain || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) {
tmp_plot_names.push_back(derived_names[i]);
}
}
if ( (solverChoice.moisture_type == MoistureType::SAM) || (derived_names[i] != "qi" &&
derived_names[i] != "qsnow" &&
derived_names[i] != "qgraup" &&
derived_names[i] != "snow_accum" &&
derived_names[i] != "graup_accum") )
{
tmp_plot_names.push_back(derived_names[i]);
} // moisture_type
} // use_terrain?
} // hasElement
}

#ifdef ERF_USE_PARTICLES
const auto& particles_namelist( particleData.getNamesUnalloc() );
for (auto it = particles_namelist.cbegin(); it != particles_namelist.cend(); ++it) {
Expand Down Expand Up @@ -967,6 +986,27 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
mf_comp += 1;
}

if (containerHasElement(plot_var_names, "qsat"))
{
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const Array4<Real>& derdat = mf[lev].array(mfi);
const Array4<Real const>& S_arr = vars_new[lev][Vars::cons].const_array(mfi);
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Real qv = S_arr(i,j,k,RhoQ1_comp) / S_arr(i,j,k,Rho_comp);
Real T = getTgivenRandRTh(S_arr(i,j,k,Rho_comp), S_arr(i,j,k,RhoTheta_comp), qv);
Real pressure = getPgivenRTh(S_arr(i,j,k,RhoTheta_comp), qv);
erf_qsatw(T, pressure, derdat(i,j,k,mf_comp));
});
}
mf_comp ++;
}

if(solverChoice.moisture_type == MoistureType::Kessler){
if (containerHasElement(plot_var_names, "rain_accum"))
{
Expand Down
43 changes: 20 additions & 23 deletions Source/Microphysics/Kessler/Kessler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,6 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice)
});
}


for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi);
auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi);
Expand Down Expand Up @@ -118,30 +117,28 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice)
//Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));

// If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature
if(qv_array(i,j,k) > qsat){
if (qv_array(i,j,k) > qsat) {
dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac));
}


// If water vapor is less than the saturated value, then the cloud water can evaporate, leading to evaporative cooling and
// reducing temperature
if(qv_array(i,j,k) < qsat && qc_array(i,j,k) > 0.0){
// If water vapor is less than the saturated value, then the cloud water can evaporate,
// leading to evaporative cooling and reducing temperature
if (qv_array(i,j,k) < qsat && qc_array(i,j,k) > 0.0) {
dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
}
}

if(qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) {
if (qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) {
Real C = 1.6 + 124.9*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.2046);
dq_rain_to_vapor = 1.0/(0.001*rho_array(i,j,k))*(1.0 - qv_array(i,j,k)/qsat)*C*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.525)/
(5.4e5 + 2.55e6/(pressure*qsat))*dtn;
// The negative sign is to make this variable (vapor formed from evaporation)
// a poistive quantity (as qv/qs < 1)
// a positive quantity (as qv/qs < 1)
dq_rain_to_vapor = std::min({qp_array(i,j,k), dq_rain_to_vapor});

// Removing latent heat due to evaporation from rain water to water vapor, reduces the (potential) temperature
}

// If there is cloud water present then do accretion and autoconversion to rain

if (qc_array(i,j,k) > 0.0) {
qcc = qc_array(i,j,k);

Expand All @@ -163,12 +160,12 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice)
Real dq_sed = dtn * dJinv * (1.0/rho_array(i,j,k)) * (fz_array(i,j,k+1) - fz_array(i,j,k))/dz;
if(std::fabs(dq_sed) < 1e-14) dq_sed = 0.0;

qv_array(i,j,k) = qv_array(i,j,k) - dq_vapor_to_clwater + dq_clwater_to_vapor + dq_rain_to_vapor;
qc_array(i,j,k) = qc_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain;
qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor;
qv_array(i,j,k) += -dq_vapor_to_clwater + dq_clwater_to_vapor + dq_rain_to_vapor;
qc_array(i,j,k) += dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain;
qp_array(i,j,k) += dq_sed + dq_clwater_to_rain - dq_rain_to_vapor;

theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond *
(dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor);
Real theta_over_T = theta_array(i,j,k)/tabs_array(i,j,k);
theta_array(i,j,k) += theta_over_T * d_fac_cond * (dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor);

qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
Expand All @@ -179,8 +176,7 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice)
}
}


if(solverChoice.moisture_type == MoistureType::Kessler_NoRain){
if (solverChoice.moisture_type == MoistureType::Kessler_NoRain){

// get the temperature, dentisy, theta, qt and qc from input
for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
Expand Down Expand Up @@ -217,20 +213,21 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice)
Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k));

// If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature
if(qv_array(i,j,k) > qsat){
if (qv_array(i,j,k) > qsat){
dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac));
}
// If water vapor is less than the saturated value, then the cloud water can evaporate, leading to evaporative cooling and
// reducing temperature
if(qv_array(i,j,k) < qsat && qc_array(i,j,k) > 0.0){
if (qv_array(i,j,k) < qsat && qc_array(i,j,k) > 0.0){
dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
}

qv_array(i,j,k) = qv_array(i,j,k) - dq_vapor_to_clwater + dq_clwater_to_vapor;
qc_array(i,j,k) = qc_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor;
qv_array(i,j,k) += -dq_vapor_to_clwater + dq_clwater_to_vapor;
qc_array(i,j,k) += dq_vapor_to_clwater - dq_clwater_to_vapor;

Real theta_over_T = theta_array(i,j,k)/tabs_array(i,j,k);

theta_array(i,j,k) = theta_array(i,j,k) + theta_array(i,j,k)/tabs_array(i,j,k)*d_fac_cond *
(dq_vapor_to_clwater - dq_clwater_to_vapor);
theta_array(i,j,k) += theta_over_T * d_fac_cond * (dq_vapor_to_clwater - dq_clwater_to_vapor);

qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
Expand Down
54 changes: 34 additions & 20 deletions Source/Microphysics/SAM/Cloud_SAM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,12 +22,12 @@ void SAM::Cloud (const SolverChoice& solverChoice) {

SolverChoice sc = solverChoice;

int moisture_type = 1;
int SAM_moisture_type = 1;

if(solverChoice.moisture_type == MoistureType::SAM) {
moisture_type = 1;
SAM_moisture_type = 1;
} else if(solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
moisture_type = 2;
SAM_moisture_type = 2;
}

for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]), TilingIfNotGPU()); mfi.isValid(); ++mfi) {
Expand Down Expand Up @@ -65,9 +65,9 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
// before the Newton iteration, which assumes it is.

omn = 1.0;
if(moisture_type == 1){
if (SAM_moisture_type == 1){
// Cloud ice not permitted (melt to form water)
if(tabs_array(i,j,k) >= tbgmax) {
if (tabs_array(i,j,k) >= tbgmax) {
omn = 1.0;
delta_qi = qci_array(i,j,k);
qci_array(i,j,k) = 0.0;
Expand All @@ -79,7 +79,7 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
pres_array(i,j,k) *= 0.01;
}
// Cloud water not permitted (freeze to form ice)
else if(tabs_array(i,j,k) <= tbgmin) {
else if (tabs_array(i,j,k) <= tbgmin) {
omn = 0.0;
delta_qc = qcl_array(i,j,k);
qcl_array(i,j,k) = 0.0;
Expand All @@ -104,7 +104,8 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
pres_array(i,j,k) *= 0.01;
}
}
else if(moisture_type == 2){
else if (SAM_moisture_type == 2)
{
// No ice. ie omn = 1.0
delta_qc = qcl_array(i,j,k) - qn_array(i,j,k);
delta_qi = 0.0;
Expand Down Expand Up @@ -145,7 +146,7 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
erf_dtqsatw(tabs, pres, dqsatw);
erf_dtqsati(tabs, pres, dqsati);

if(moisture_type == 1) {
if (SAM_moisture_type == 1) {
// Cloud ice not permitted (condensation & fusion)
if(tabs >= tbgmax) {
omn = 1.0;
Expand All @@ -160,7 +161,7 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
omn = an*tabs-bn;
domn = an;
}
} else if(moisture_type == 2) {
} else if (SAM_moisture_type == 2) {
omn = 1.0;
domn = 0.0;
}
Expand All @@ -183,16 +184,18 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
dtabs = -fff/dfff;
tabs = tabs+dtabs;

// Update the pressure
pres = rho_array(i,j,k) * R_d * tabs
* (1.0 + R_v/R_d * qsat) * 0.01;
// For now at least we perform this iteration at constant pressure
// For the moist bubble case, the results are indistinguisable
// between running with this used vs commented out
// pres = rho_array(i,j,k) * R_d * tabs
// * (1.0 + R_v/R_d * qsat) * 0.01;

// Update iteration
niter = niter+1;
} while(std::abs(dtabs) > tol && niter < 20);

// Update qsat from last iteration (dq = dq/dt * dt)
qsat = qsat + dqsat*dtabs;
qsat += dqsat*dtabs;

// Changes in each component
delta_qv = qv_array(i,j,k) - qsat;
Expand All @@ -209,19 +212,30 @@ void SAM::Cloud (const SolverChoice& solverChoice) {
// Update temperature
tabs_array(i,j,k) = tabs;

// Update pressure
pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
* (1.0 + R_v/R_d * qv_array(i,j,k));
// Update theta from temperature (it is essential to do this BEFORE the pressure is updated)
// This would be inconsistent with updating the pressure as part of the iteration above.
// Empirically based on the moist bubble rise case, getting the correct theta here
// depends on using the old (unchanged by the phase changes) pressure.

// Update theta from temperature
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);

// Update pressure to be consistent with updated theta_array
pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k));

// This was used in the earlier implmentation when we updated theta using this new pressure
// pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
// * (1.0 + R_v/R_d * qv_array(i,j,k));

// Pressure unit conversion
pres_array(i,j,k) *= 0.01;

}
//
// We cannot blindly relax to qsat, but we can convert qc/qi -> qv
else {
// The concept here is that if we assume the temperature change due to this conversion
// doesn't change qsat, we can safely put all the moisture into qv without reaching qv > qsat
// because we are in the "else" part of the test on (qt_array(i,j,k) > qsat)
//
} else {
// Changes in each component
delta_qv = qcl_array(i,j,k) + qci_array(i,j,k);
delta_qc = qcl_array(i,j,k);
Expand Down

0 comments on commit fb65da6

Please sign in to comment.