diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp
index dc32173d9..8ac3006f1 100644
--- a/Source/BoundaryConditions/ERF_FillPatch.cpp
+++ b/Source/BoundaryConditions/ERF_FillPatch.cpp
@@ -15,7 +15,6 @@ PhysBCFunctNoOp null_bc;
* @param[in] time time at which the data should be filled
* @param[out] mfs Vector of MultiFabs to be filled containing, in order: cons, xvel, yvel, and zvel data
*/
-
void
ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs, bool fillset)
{
@@ -100,6 +99,10 @@ ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs, bool fillset)
// We call this even if init_type == real because this routine will fill the vertical bcs
(*physbcs[lev])(mfs,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,init_type,cons_only,BCVars::cons_bc,time);
+
+ // Update vars and pointers if we have a micro model
+ if (solverChoice.moisture_type != MoistureType::None)
+ micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
}
/*
@@ -109,7 +112,6 @@ ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs, bool fillset)
* @param[in] time time at which the data should be filled
* @param[out] mf MultiFab to be filled (qmoist[lev])
*/
-
void
ERF::FillPatchMoistVars (int lev, MultiFab& mf)
{
@@ -151,7 +153,6 @@ ERF::FillPatchMoistVars (int lev, MultiFab& mf)
* @param[in] eddyDiffs diffusion coefficients for LES turbulence models
* @param[in] allow_most_bcs if true then use MOST bcs at the low boundary
*/
-
void
ERF::FillIntermediatePatch (int lev, Real time,
const Vector<MultiFab*>& mfs,
@@ -260,9 +261,11 @@ ERF::FillIntermediatePatch (int lev, Real time,
// MOST boundary conditions
if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs)
m_most->impose_most_bcs(lev,mfs,eddyDiffs_lev[lev].get(),z_phys_nd[lev].get());
-}
-//
+ // Update vars and pointers if we have a micro model
+ if (solverChoice.moisture_type != MoistureType::None)
+ micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
+}
/*
* Fill valid and ghost data.
@@ -274,7 +277,6 @@ ERF::FillIntermediatePatch (int lev, Real time,
* @param[in] time time at which the data should be filled
* @param[out] mfs Vector of MultiFabs to be filled containing, in order: cons, xvel, yvel, and zvel data
*/
-
void
ERF::FillCoarsePatch (int lev, Real time, const Vector<MultiFab*>& mfs)
{
diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H
index 01276a7b8..a3b31199a 100644
--- a/Source/DataStructs/DataStruct.H
+++ b/Source/DataStructs/DataStruct.H
@@ -87,7 +87,11 @@ struct SolverChoice {
// Set a different default for moist vs dry
if (moisture_type != MoistureType::None) {
- buoyancy_type = 2; // uses Tprime
+ if (moisture_type == MoistureType::FastEddy) {
+ buoyancy_type = 1; // asserted in make buoyancy
+ } else {
+ buoyancy_type = 2; // uses Tprime
+ }
}
// Is the terrain static or moving?
diff --git a/Source/ERF.H b/Source/ERF.H
index 122cb0116..c6ff0ebe7 100644
--- a/Source/ERF.H
+++ b/Source/ERF.H
@@ -525,7 +525,7 @@ private:
amrex::Vector<amrex::MultiFab> rW_new;
Microphysics micro;
- amrex::Vector<std::unique_ptr<amrex::MultiFab>> qmoist; // This has up to 6 components: qv, qc, qi, qr, qs, qg
+ amrex::Vector<amrex::Vector<amrex::MultiFab*>> qmoist; // (lev,ncomp) This has up to 6 components: qv, qc, qi, qr, qs, qg
#if defined(ERF_USE_RRTMGP)
Radiation rad;
diff --git a/Source/ERF.cpp b/Source/ERF.cpp
index 000066d67..53a1f65c1 100644
--- a/Source/ERF.cpp
+++ b/Source/ERF.cpp
@@ -109,6 +109,12 @@ ERF::ERF ()
amrex::Print() << "\n";
}
+ int nlevs_max = max_level + 1;
+
+ // NOTE: size micro before readparams (chooses the model at all levels)
+ micro.ReSize(nlevs_max);
+ qmoist.resize(nlevs_max);
+
ReadParameters();
const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1);
const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2);
@@ -128,8 +134,6 @@ ERF::ERF ()
// No valid BoxArray and DistributionMapping have been defined.
// But the arrays for them have been resized.
- int nlevs_max = max_level + 1;
-
istep.resize(nlevs_max, 0);
nsubsteps.resize(nlevs_max, 1);
for (int lev = 1; lev <= max_level; ++lev) {
@@ -157,9 +161,6 @@ ERF::ERF ()
vars_old[lev].resize(Vars::NumTypes);
}
- micro.ReSize(nlevs_max);
- qmoist.resize(nlevs_max);
-
mri_integrator_mem.resize(nlevs_max);
physbcs.resize(nlevs_max);
@@ -218,8 +219,7 @@ ERF::ERF ()
}
}
-ERF::~ERF ()
-= default;
+ERF::~ERF () = default;
// advance solution to final time
void
@@ -530,7 +530,6 @@ ERF::InitData ()
}
}
-
if (solverChoice.coupling_type == CouplingType::TwoWay) {
int src_comp_reflux = 0;
int num_comp_reflux = vars_new[0][Vars::cons].nComp();
@@ -549,7 +548,7 @@ ERF::InitData ()
if (solverChoice.moisture_type != MoistureType::None)
{
for (int lev = 0; lev <= finest_level; lev++) {
- FillPatchMoistVars(lev, *(qmoist[lev]));
+ FillPatchMoistVars(lev, *(qmoist[lev][0])); // qv component
}
}
}
@@ -579,20 +578,6 @@ ERF::InitData ()
}
}
-
- for (int lev = 0; lev <= finest_level; ++lev)
- {
- // Initialize microphysics here
- micro.define(lev, solverChoice);
-
- // If not restarting we need to fill qmoist given qt and qp.
- if (solverChoice.moisture_type != MoistureType::None)
- {
- micro.Init(lev, vars_new[lev][Vars::cons], *(qmoist[lev]), grids[lev], Geom(lev), 0.0); // dummy value, not needed just to diagnose
- micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
- }
- }
-
// Configure ABLMost params if used MostWall boundary condition
// NOTE: we must set up the MOST routine before calling WritePlotFile because
// WritePlotFile calls FillPatch in order to compute gradients
@@ -815,6 +800,7 @@ ERF::InitData ()
void
ERF::restart ()
{
+ // TODO: This could be deleted since ba/dm are not created yet?
for (int lev = 0; lev <= finest_level; ++lev)
{
auto& lev_new = vars_new[lev];
@@ -857,8 +843,6 @@ ERF::init_only (int lev, Real time)
lev_new[Vars::yvel].setVal(0.0); lev_old[Vars::yvel].setVal(0.0);
lev_new[Vars::zvel].setVal(0.0); lev_old[Vars::zvel].setVal(0.0);
- qmoist[lev]->setVal(0.);
-
// Initialize background flow (optional)
if (init_type == "input_sounding") {
// The base state is initialized by integrating vertically through the
@@ -1175,13 +1159,11 @@ ERF::MakeHorizontalAverages ()
if (use_moisture)
{
- MultiFab qv(*(qmoist[lev]), make_alias, 0, 1);
-
for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
const Box& bx = mfi.validbox();
auto fab_arr = mf.array(mfi);
auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi);
- auto const qv_arr = qv.const_array(mfi);
+ auto const qv_arr = qmoist[lev][0]->const_array(mfi);
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
Real dens = cons_arr(i, j, k, Rho_comp);
@@ -1462,6 +1444,12 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
amrex::Print() << "\n";
}
+ int nlevs_max = max_level + 1;
+
+ // NOTE: size micro before readparams (chooses the model at all levels)
+ micro.ReSize(nlevs_max);
+ qmoist.resize(nlevs_max);
+
ReadParameters();
const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1);
const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2);
@@ -1473,8 +1461,6 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
// No valid BoxArray and DistributionMapping have been defined.
// But the arrays for them have been resized.
- int nlevs_max = max_level + 1;
-
istep.resize(nlevs_max, 0);
nsubsteps.resize(nlevs_max, 1);
for (int lev = 1; lev <= max_level; ++lev) {
@@ -1502,9 +1488,6 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in,
rV_old.resize(nlevs_max);
rW_old.resize(nlevs_max);
- micro.ReSize(nlevs_max);
- qmoist.resize(nlevs_max);
-
mri_integrator_mem.resize(nlevs_max);
physbcs.resize(nlevs_max);
diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp
index cbd27688e..7d9c77e7a 100644
--- a/Source/ERF_make_new_level.cpp
+++ b/Source/ERF_make_new_level.cpp
@@ -62,7 +62,15 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
// Microphysics
// *******************************************************************************************
int q_size = micro.Get_Qmoist_Size();
- qmoist[lev] = std::make_unique<MultiFab>(ba, dm, q_size, ngrow_state); // qv, qc, qi, qr, qs, qg
+ qmoist[lev].resize(q_size);
+ micro.Define(lev, solverChoice);
+ if (solverChoice.moisture_type != MoistureType::None)
+ {
+ micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
+ }
+ for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
+ qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
+ }
// ********************************************************************************************
// Build the data structures for calculating diffusive/turbulent terms
@@ -198,9 +206,16 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
//********************************************************************************************
// Microphysics
// *******************************************************************************************
- int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1;
int q_size = micro.Get_Qmoist_Size();
- qmoist[lev] = std::make_unique<MultiFab>(ba, dm, q_size, ngrow_state);
+ qmoist[lev].resize(q_size);
+ micro.Define(lev, solverChoice);
+ if (solverChoice.moisture_type != MoistureType::None)
+ {
+ micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
+ }
+ for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
+ qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
+ }
init_stuff(lev, ba, dm);
@@ -305,7 +320,15 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
// Microphysics
// *******************************************************************************************
int q_size = micro.Get_Qmoist_Size();
- qmoist[lev]->define(ba, dm, q_size, ngrow_state); // qv, qc, qi, qr, qs, qg
+ qmoist[lev].resize(q_size);
+ micro.Define(lev, solverChoice);
+ if (solverChoice.moisture_type != MoistureType::None)
+ {
+ micro.Init(lev, vars_new[lev][Vars::cons], grids[lev], Geom(lev), 0.0); // dummy dt value
+ }
+ for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
+ qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
+ }
init_stuff(lev,ba,dm);
@@ -556,9 +579,6 @@ ERF::ClearLevel (int lev)
rW_new[lev].clear();
rW_old[lev].clear();
- // Clears the qmoist memory
- qmoist[lev].reset();
-
// Clears the integrator memory
mri_integrator_mem[lev].reset();
physbcs[lev].reset();
diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp
index d44153128..dce2ee8a1 100644
--- a/Source/IO/Checkpoint.cpp
+++ b/Source/IO/Checkpoint.cpp
@@ -130,11 +130,14 @@ ERF::WriteCheckpointFile () const
IntVect ng;
if (solverChoice.moisture_type != MoistureType::None) {
- // We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars
- ng = qmoist[lev]->nGrowVect();
- MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev]->nComp(),ng);
- MultiFab::Copy(moist_vars,*(qmoist[lev]),0,0,qmoist[lev]->nComp(),ng);
- VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars"));
+ for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
+ // We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars
+ ng = qmoist[lev][mvar]->nGrowVect();
+ int nvar = qmoist[lev][mvar]->nComp();
+ MultiFab moist_vars(grids[lev],dmap[lev],nvar,ng);
+ MultiFab::Copy(moist_vars,*(qmoist[lev][mvar]),0,0,nvar,ng);
+ VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars"));
+ }
}
// Note that we write the ghost cells of the base state (unlike above)
@@ -300,7 +303,6 @@ ERF::ReadCheckpointFile ()
}
for (int lev = 0; lev <= finest_level; ++lev) {
-
// read in level 'lev' BoxArray from Header
BoxArray ba;
ba.readFrom(is);
@@ -338,13 +340,14 @@ ERF::ReadCheckpointFile ()
MultiFab::Copy(vars_new[lev][Vars::zvel],zvel,0,0,1,0);
IntVect ng;
- if (solverChoice.moisture_type == MoistureType::None) {
- qmoist[lev]->setVal(0.0);
- } else {
- ng = qmoist[lev]->nGrowVect();
- MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev]->nComp(),ng);
- VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars"));
- MultiFab::Copy(*(qmoist[lev]),moist_vars,0,0,qmoist[lev]->nComp(),ng);
+ if (solverChoice.moisture_type != MoistureType::None) {
+ for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
+ ng = qmoist[lev][mvar]->nGrowVect();
+ int nvar = qmoist[lev][mvar]->nComp();
+ MultiFab moist_vars(grids[lev],dmap[lev],nvar,ng);
+ VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars"));
+ MultiFab::Copy(*(qmoist[lev][mvar]),moist_vars,0,0,nvar,ng);
+ }
}
// Note that we read the ghost cells of the base state (unlike above)
diff --git a/Source/IO/ERF_Write1DProfiles.cpp b/Source/IO/ERF_Write1DProfiles.cpp
index 3b5c6ca24..bd06f78f2 100644
--- a/Source/IO/ERF_Write1DProfiles.cpp
+++ b/Source/IO/ERF_Write1DProfiles.cpp
@@ -246,18 +246,17 @@ ERF::derive_diag_profiles(Gpu::HostVector<Real>& h_avg_u , Gpu::HostVector<Rea
if (use_moisture)
{
- MultiFab qv(*(qmoist[lev]), make_alias, 0, 1);
for ( MFIter mfi(mf_cons,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
- const Array4<Real>& fab_arr = mf_out.array(mfi);
+ const Array4<Real>& fab_arr = mf_out.array(mfi);
const Array4<Real>& cons_arr = mf_cons.array(mfi);
const Array4<Real>& u_cc_arr = u_cc.array(mfi);
const Array4<Real>& v_cc_arr = v_cc.array(mfi);
const Array4<Real>& w_cc_arr = w_cc.array(mfi);
const Array4<Real>& p0_arr = p_hse.array(mfi);
- const Array4<Real>& qv_arr = qv.array(mfi);
+ const Array4<Real>& qv_arr = qmoist[0][0]->array(mfi); // TODO: Is this written only on lev 0?
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp
index 3b1347b8d..b237edd0e 100644
--- a/Source/IO/Plotfile.cpp
+++ b/Source/IO/Plotfile.cpp
@@ -119,6 +119,14 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
fillset);
}
+ // Get qmoist pointers if using moisture
+ bool l_use_moisture = (solverChoice.moisture_type != MoistureType::None);
+ for (int lev = 0; lev <= finest_level; ++lev) {
+ for (int mvar(0); mvar<qmoist[lev].size(); ++mvar) {
+ qmoist[lev][mvar] = micro.Get_Qmoist_Ptr(lev,mvar);
+ }
+ }
+
Vector<MultiFab> mf(finest_level+1);
for (int lev = 0; lev <= finest_level; ++lev) {
mf[lev].define(grids[lev], dmap[lev], ncomp_mf, 0);
@@ -211,11 +219,12 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
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);
- const Array4<Real const>& qv_arr = qmoist[0]->const_array(mfi);
+ const Array4<Real const>& qv_arr = (qmoist[lev][0]) ? qmoist[lev][0]->const_array(mfi) :
+ Array4<Real const> {};
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
- Real qv_for_p = qv_arr(i,j,k);
+ Real qv_for_p = (qv_arr) ? qv_arr(i,j,k) : 0;
const Real rhotheta = S_arr(i,j,k,RhoTheta_comp);
derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta,qv_for_p);
});
@@ -233,11 +242,12 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
const Array4<Real>& derdat = mf[lev].array(mfi);
const Array4<Real const>& p0_arr = p_hse.const_array(mfi);
const Array4<Real const>& S_arr = vars_new[lev][Vars::cons].const_array(mfi);
- const Array4<Real const> & qv_arr = qmoist[0]->const_array(mfi);
+ const Array4<Real const>& qv_arr = (qmoist[lev][0]) ? qmoist[lev][0]->const_array(mfi) :
+ Array4<Real const> {};
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
- Real qv_for_p = qv_arr(i,j,k);
+ Real qv_for_p = (qv_arr) ? qv_arr(i,j,k) : 0;
const Real rhotheta = S_arr(i,j,k,RhoTheta_comp);
derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta,qv_for_p) - p0_arr(i,j,k);
});
@@ -600,69 +610,71 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
}
// NOTE: Protect against accessing non-existent data
- int q_size = qmoist[lev]->nComp();
- MultiFab qv_mf(*(qmoist[lev]), make_alias, 0, 1);
+ if (l_use_moisture) {
+ int q_size = qmoist[lev].size();
+ MultiFab qv_mf(*(qmoist[lev][0]), make_alias, 0, 1);
- if (containerHasElement(plot_var_names, "qt") && (q_size >= 3))
- {
- MultiFab qc_mf(*(qmoist[lev]), make_alias, 1, 1);
- MultiFab qi_mf(*(qmoist[lev]), make_alias, 2, 1);
- MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0);
- MultiFab::Add (mf[lev],qc_mf,0,mf_comp,1,0);
- MultiFab::Add (mf[lev],qi_mf,0,mf_comp,1,0);
- mf_comp += 1;
- }
+ if (containerHasElement(plot_var_names, "qt") && (q_size >= 3))
+ {
+ MultiFab qc_mf(*(qmoist[lev][1]), make_alias, 0, 1);
+ MultiFab qi_mf(*(qmoist[lev][2]), make_alias, 0, 1);
+ MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0);
+ MultiFab::Add (mf[lev],qc_mf,0,mf_comp,1,0);
+ MultiFab::Add (mf[lev],qi_mf,0,mf_comp,1,0);
+ mf_comp += 1;
+ }
- if (containerHasElement(plot_var_names, "qp") && (q_size >= 6))
- {
- MultiFab qr_mf(*(qmoist[lev]), make_alias, 3, 1);
- MultiFab qs_mf(*(qmoist[lev]), make_alias, 4, 1);
- MultiFab qg_mf(*(qmoist[lev]), make_alias, 5, 1);
- MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0);
- MultiFab::Add (mf[lev],qs_mf,0,mf_comp,1,0);
- MultiFab::Add (mf[lev],qg_mf,0,mf_comp,1,0);
- mf_comp += 1;
- }
+ if (containerHasElement(plot_var_names, "qp") && (q_size >= 6))
+ {
+ MultiFab qr_mf(*(qmoist[lev][3]), make_alias, 0, 1);
+ MultiFab qs_mf(*(qmoist[lev][4]), make_alias, 0, 1);
+ MultiFab qg_mf(*(qmoist[lev][5]), make_alias, 0, 1);
+ MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0);
+ MultiFab::Add (mf[lev],qs_mf,0,mf_comp,1,0);
+ MultiFab::Add (mf[lev],qg_mf,0,mf_comp,1,0);
+ mf_comp += 1;
+ }
- if (containerHasElement(plot_var_names, "qv"))
- {
- MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0);
- mf_comp += 1;
- }
+ if (containerHasElement(plot_var_names, "qv"))
+ {
+ MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0);
+ mf_comp += 1;
+ }
- if (containerHasElement(plot_var_names, "qc") && (q_size >= 2))
- {
- MultiFab qc_mf(*(qmoist[lev]), make_alias, 1, 1);
- MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0);
- mf_comp += 1;
- }
+ if (containerHasElement(plot_var_names, "qc") && (q_size >= 2))
+ {
+ MultiFab qc_mf(*(qmoist[lev][1]), make_alias, 0, 1);
+ MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0);
+ mf_comp += 1;
+ }
- if (containerHasElement(plot_var_names, "qi") && (q_size >= 3))
- {
- MultiFab qi_mf(*(qmoist[lev]), make_alias, 2, 1);
- MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0);
- mf_comp += 1;
- }
+ if (containerHasElement(plot_var_names, "qi") && (q_size >= 3))
+ {
+ MultiFab qi_mf(*(qmoist[lev][2]), make_alias, 0, 1);
+ MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0);
+ mf_comp += 1;
+ }
- if (containerHasElement(plot_var_names, "qrain") && (q_size >= 4))
- {
- MultiFab qr_mf(*(qmoist[lev]), make_alias, 3, 1);
- MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0);
- mf_comp += 1;
- }
+ if (containerHasElement(plot_var_names, "qrain") && (q_size >= 4))
+ {
+ MultiFab qr_mf(*(qmoist[lev][3]), make_alias, 0, 1);
+ MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0);
+ mf_comp += 1;
+ }
- if (containerHasElement(plot_var_names, "qsnow") && (q_size >= 5))
- {
- MultiFab qs_mf(*(qmoist[lev]), make_alias, 4, 1);
- MultiFab::Copy(mf[lev],qs_mf,0,mf_comp,1,0);
- mf_comp += 1;
- }
+ if (containerHasElement(plot_var_names, "qsnow") && (q_size >= 5))
+ {
+ MultiFab qs_mf(*(qmoist[lev][4]), make_alias, 0, 1);
+ MultiFab::Copy(mf[lev],qs_mf,0,mf_comp,1,0);
+ mf_comp += 1;
+ }
- if (containerHasElement(plot_var_names, "qgraup") && (q_size >= 6))
- {
- MultiFab qg_mf(*(qmoist[lev]), make_alias, 5, 1);
- MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0);
- mf_comp += 1;
+ if (containerHasElement(plot_var_names, "qgraup") && (q_size >= 6))
+ {
+ MultiFab qg_mf(*(qmoist[lev][5]), make_alias, 0, 1);
+ MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0);
+ mf_comp += 1;
+ }
}
#ifdef ERF_USE_PARTICLES
diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp
index 2e2fb0586..94a832c7f 100644
--- a/Source/Initialization/ERF_init_custom.cpp
+++ b/Source/Initialization/ERF_init_custom.cpp
@@ -24,8 +24,7 @@ using namespace amrex;
void
ERF::init_custom (int lev)
{
- auto& lev_new = vars_new[lev];
- auto& qmoist_new = *(qmoist[lev]);
+ auto& lev_new = vars_new[lev];
MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component
MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component
@@ -42,12 +41,12 @@ ERF::init_custom (int lev)
yvel_pert.setVal(0.);
zvel_pert.setVal(0.);
- MultiFab qmoist_pert(qmoist[lev]->boxArray(), qmoist[lev]->DistributionMap(), 3, qmoist[lev]->nGrow());
- qmoist_pert.setVal(0.);
-
- MultiFab qv_pert(qmoist_pert, amrex::make_alias, 0, 1);
- MultiFab qc_pert(qmoist_pert, amrex::make_alias, 1, 1);
- MultiFab qi_pert(qmoist_pert, amrex::make_alias, 2, 1);
+ MultiFab qv_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), 1, lev_new[Vars::cons].nGrow());
+ MultiFab qc_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), 1, lev_new[Vars::cons].nGrow());
+ MultiFab qi_pert(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), 1, lev_new[Vars::cons].nGrow());
+ qv_pert.setVal(0.);
+ qc_pert.setVal(0.);
+ qi_pert.setVal(0.);
int fix_random_seed = 0;
ParmParse pp("erf"); pp.query("fix_random_seed", fix_random_seed);
@@ -108,11 +107,11 @@ ERF::init_custom (int lev)
if (solverChoice.moisture_type != MoistureType::None) {
MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ1_comp, RhoQ1_comp, 1, cons_pert.nGrow());
MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ2_comp, RhoQ2_comp, 1, cons_pert.nGrow());
- MultiFab::Add( qmoist_new, qmoist_pert, 0, 0, 3, qmoist_pert.nGrow()); // qv, qc, qi
-#if defined(ERF_USE_WARM_NO_PRECIP)
- MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQv_comp, RhoQv_comp, 1, cons_pert.nGrow());
- MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQc_comp, RhoQc_comp, 1, cons_pert.nGrow());
-#endif
+
+ MultiFab::Add(*(qmoist[lev][0]) , qv_pert , 0, 0, 1, qv_pert.nGrow());
+ MultiFab::Add(*(qmoist[lev][1]) , qc_pert , 0, 0, 1, qc_pert.nGrow());
+ if (qmoist[lev].size() > 2)
+ MultiFab::Add(*(qmoist[lev][2]) , qi_pert , 0, 0, 1, qi_pert.nGrow());
}
MultiFab::Add(lev_new[Vars::xvel], xvel_pert, 0, 0, 1, xvel_pert.nGrowVect());
diff --git a/Source/Microphysics/FastEddy/FastEddy.H b/Source/Microphysics/FastEddy/FastEddy.H
index dc54a8225..11a59a1d7 100644
--- a/Source/Microphysics/FastEddy/FastEddy.H
+++ b/Source/Microphysics/FastEddy/FastEddy.H
@@ -51,7 +51,7 @@ public:
// Set up for first time
void
- define (SolverChoice& sc) override
+ Define (SolverChoice& sc) override
{
docloud = sc.do_cloud;
doprecip = sc.do_precip;
@@ -65,7 +65,6 @@ public:
// init
void
Init (const amrex::MultiFab& cons_in,
- amrex::MultiFab& qmoist,
const amrex::BoxArray& grids,
const amrex::Geometry& geom,
const amrex::Real& dt_advance) override;
@@ -98,8 +97,10 @@ public:
// wrapper to do all the updating
void
- Advance ( ) override
+ Advance (const amrex::Real& dt_advance) override
{
+ dt = dt_advance;
+
this->AdvanceFE();
this->Diagnose();
}
@@ -107,17 +108,20 @@ public:
amrex::MultiFab*
Qmoist_Ptr (const int& varIdx) override
{
- AMREX_ALWAYS_ASSERT(varIdx <= m_qmoist_size);
- return mic_fab_vars[varIdx].get();
+ AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
+ return mic_fab_vars[MicVarMap[varIdx]].get();
}
int
- Qmoist_Size ( ) override { return FastEddy::m_qmoist_size; }
+ Qmoist_Size () override { return FastEddy::m_qmoist_size; }
private:
// Number of qmoist variables
int m_qmoist_size = 2;
+ // MicVar map (Qmoist indices -> MicVar enum)
+ amrex::Vector<int> MicVarMap;
+
// geometry
amrex::Geometry m_geom;
diff --git a/Source/Microphysics/FastEddy/Init_FE.cpp b/Source/Microphysics/FastEddy/Init_FE.cpp
index c31d75390..1e27d4e49 100644
--- a/Source/Microphysics/FastEddy/Init_FE.cpp
+++ b/Source/Microphysics/FastEddy/Init_FE.cpp
@@ -19,7 +19,7 @@ using namespace amrex;
* @param[in] geom Geometry associated with these MultiFabs and grids
* @param[in] dt_advance Timestep for the advance
*/
-void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist,
+void FastEddy::Init (const MultiFab& cons_in,
const BoxArray& grids,
const Geometry& geom,
const Real& dt_advance)
@@ -28,9 +28,13 @@ void FastEddy::Init (const MultiFab& cons_in, MultiFab& qmoist,
m_geom = geom;
m_gtoe = grids;
+ MicVarMap.resize(m_qmoist_size);
+ MicVarMap = {MicVar_FE::qv, MicVar_FE::qc};
+
// initialize microphysics variables
for (auto ivar = 0; ivar < MicVar_FE::NumVars; ++ivar) {
- mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect());
+ mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
+ 1, cons_in.nGrowVect());
mic_fab_vars[ivar]->setVal(0.);
}
diff --git a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp
index 1faf625de..55b613ca3 100644
--- a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp
+++ b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp
@@ -5,8 +5,8 @@
* Computes diagnostic quantities like cloud ice/liquid and precipitation ice/liquid
* from the existing Microphysics variables.
*/
-void Kessler::Diagnose () {
-
+void Kessler::Diagnose ()
+{
auto qt = mic_fab_vars[MicVar_Kess::qt];
auto qp = mic_fab_vars[MicVar_Kess::qp];
auto qv = mic_fab_vars[MicVar_Kess::qv];
@@ -47,8 +47,8 @@ void Kessler::Diagnose () {
* Computes diagnostic quantities like temperature and pressure
* from the existing Microphysics variables.
*/
-void Kessler::Compute_Tabs_and_Pres () {
-
+void Kessler::Compute_Tabs_and_Pres ()
+{
auto rho = mic_fab_vars[MicVar_Kess::rho];
auto theta = mic_fab_vars[MicVar_Kess::theta];
auto qv = mic_fab_vars[MicVar_Kess::qv];
diff --git a/Source/Microphysics/Kessler/Init_Kessler.cpp b/Source/Microphysics/Kessler/Init_Kessler.cpp
index 5e68d3e76..56f8c6be4 100644
--- a/Source/Microphysics/Kessler/Init_Kessler.cpp
+++ b/Source/Microphysics/Kessler/Init_Kessler.cpp
@@ -21,7 +21,6 @@ using namespace amrex;
* @param[in] dt_advance Timestep for the advance
*/
void Kessler::Init (const MultiFab& cons_in,
- MultiFab& qmoist,
const BoxArray& grids,
const Geometry& geom,
const Real& dt_advance)
@@ -30,15 +29,19 @@ void Kessler::Init (const MultiFab& cons_in,
m_geom = geom;
m_gtoe = grids;
+ MicVarMap.resize(m_qmoist_size);
+ MicVarMap = {MicVar_Kess::qv, MicVar_Kess::qcl, MicVar_Kess::qci,
+ MicVar_Kess::qr, MicVar_Kess::qpi, MicVar_Kess::qg};
+
// initialize microphysics variables
- for (auto ivar = 0; ivar < MicVar_FE::NumVars; ++ivar) {
- mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(), 1, cons_in.nGrowVect());
+ for (auto ivar = 0; ivar < MicVar_Kess::NumVars; ++ivar) {
+ mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
+ 1, cons_in.nGrowVect());
mic_fab_vars[ivar]->setVal(0.);
}
// Set class data members
for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) {
-
const auto& box3d = mfi.tilebox();
const auto& lo = amrex::lbound(box3d);
diff --git a/Source/Microphysics/Kessler/Kessler.H b/Source/Microphysics/Kessler/Kessler.H
index 22bbd4629..00ade80ef 100644
--- a/Source/Microphysics/Kessler/Kessler.H
+++ b/Source/Microphysics/Kessler/Kessler.H
@@ -40,6 +40,7 @@ namespace MicVar_Kess {
qcl, // cloud water
qpl, // precip rain
qpi, // precip ice
+ qg, // graupel
// temporary variable
omega,
NumVars
@@ -68,7 +69,7 @@ public:
// Set up for first time
void
- define (SolverChoice& sc) override
+ Define (SolverChoice& sc) override
{
docloud = sc.do_cloud;
doprecip = sc.do_precip;
@@ -82,7 +83,6 @@ public:
// init
void
Init (const amrex::MultiFab& cons_in,
- amrex::MultiFab& qmoist,
const amrex::BoxArray& grids,
const amrex::Geometry& geom,
const amrex::Real& dt_advance) override;
@@ -118,8 +118,10 @@ public:
// wrapper to advance micro vars
void
- Advance ( ) override
+ Advance (const amrex::Real& dt_advance) override
{
+ dt = dt_advance;
+
this->AdvanceKessler();
this->Diagnose();
this->Compute_Tabs_and_Pres();
@@ -128,20 +130,23 @@ public:
amrex::MultiFab*
Qmoist_Ptr (const int& varIdx) override
{
- AMREX_ALWAYS_ASSERT(varIdx <= m_qmoist_size);
- return mic_fab_vars[varIdx].get();
+ AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
+ return mic_fab_vars[MicVarMap[varIdx]].get();
}
void
Compute_Tabs_and_Pres ();
int
- Qmoist_Size ( ) override { return Kessler::m_qmoist_size; }
+ Qmoist_Size () override { return Kessler::m_qmoist_size; }
- private:
+private:
// Number of qmoist variables
int m_qmoist_size = 6;
+ // MicVar map (Qmoist indices -> MicVar enum)
+ amrex::Vector<int> MicVarMap;
+
// geometry
amrex::Geometry m_geom;
diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp
index 9696e0981..cb5ca7dfe 100644
--- a/Source/Microphysics/Kessler/Kessler.cpp
+++ b/Source/Microphysics/Kessler/Kessler.cpp
@@ -11,8 +11,8 @@ using namespace amrex;
/**
* Compute Precipitation-related Microphysics quantities.
*/
-void Kessler::AdvanceKessler () {
-
+void Kessler::AdvanceKessler ()
+{
auto qt = mic_fab_vars[MicVar_Kess::qt];
auto qp = mic_fab_vars[MicVar_Kess::qp];
auto qn = mic_fab_vars[MicVar_Kess::qn];
@@ -42,28 +42,28 @@ void Kessler::AdvanceKessler () {
ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
- Real rho_avg, qp_avg;
-
- if (k==k_lo) {
- rho_avg = rho_array(i,j,k);
- qp_avg = qp_array(i,j,k);
- } else if (k==k_hi+1) {
- rho_avg = rho_array(i,j,k-1);
- qp_avg = qp_array(i,j,k-1);
- } else {
- rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3
- qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k));
- }
+ Real rho_avg, qp_avg;
+
+ if (k==k_lo) {
+ rho_avg = rho_array(i,j,k);
+ qp_avg = qp_array(i,j,k);
+ } else if (k==k_hi+1) {
+ rho_avg = rho_array(i,j,k-1);
+ qp_avg = qp_array(i,j,k-1);
+ } else {
+ rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3
+ qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k));
+ }
- qp_avg = std::max(0.0, qp_avg);
+ qp_avg = std::max(0.0, qp_avg);
- Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s
+ Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s
- fz_array(i,j,k) = rho_avg*V_terminal*qp_avg;
+ fz_array(i,j,k) = rho_avg*V_terminal*qp_avg;
- /*if(k==0){
- fz_array(i,j,k) = 0;
- }*/
+ /*if(k==0){
+ fz_array(i,j,k) = 0;
+ }*/
});
}
diff --git a/Source/Microphysics/Kessler/Update_Kessler.cpp b/Source/Microphysics/Kessler/Update_Kessler.cpp
index 8f6194262..cc085f9e0 100644
--- a/Source/Microphysics/Kessler/Update_Kessler.cpp
+++ b/Source/Microphysics/Kessler/Update_Kessler.cpp
@@ -80,7 +80,7 @@ void Kessler::Copy_Micro_to_State (amrex::MultiFab& cons)
// get potential total density, temperature, qt, qp
amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
- states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k);
+ //states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k);
states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k);
states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k);
states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qp_arr(i,j,k);
diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H
index 12efedf80..4113fe8db 100644
--- a/Source/Microphysics/Microphysics.H
+++ b/Source/Microphysics/Microphysics.H
@@ -12,7 +12,7 @@ public:
Microphysics () { }
- ~Microphysics () = default;
+ ~Microphysics () = default;
void
ReSize (const int& nlev) { m_moist_model.resize(nlev); }
@@ -27,27 +27,26 @@ public:
}
void
- define (const int& lev,
+ Define (const int& lev,
SolverChoice& sc)
{
- m_moist_model[lev]->define(sc);
+ m_moist_model[lev]->Define(sc);
}
void
Init (const int& lev,
const amrex::MultiFab& cons_in,
- amrex::MultiFab& qmoist,
const amrex::BoxArray& grids,
const amrex::Geometry& geom,
const amrex::Real& dt_advance)
{
- m_moist_model[lev]->Init(cons_in, qmoist, grids, geom, dt_advance);
+ m_moist_model[lev]->Init(cons_in, grids, geom, dt_advance);
}
void
- Advance (const int& lev)
+ Advance (const int& lev, const amrex::Real& dt_advance)
{
- m_moist_model[lev]->Advance();
+ m_moist_model[lev]->Advance(dt_advance);
}
void
@@ -80,7 +79,7 @@ public:
Get_Qmoist_Ptr (const int& lev, const int& varIdx) { return m_moist_model[lev]->Qmoist_Ptr(varIdx); }
int
- Get_Qmoist_Size ( ) { return m_moist_model[0]->Qmoist_Size(); }
+ Get_Qmoist_Size () { return m_moist_model[0]->Qmoist_Size(); }
private:
amrex::Vector<std::unique_ptr<NullMoist>> m_moist_model;
diff --git a/Source/Microphysics/Null/NullMoist.H b/Source/Microphysics/Null/NullMoist.H
index 1c47de49c..5dfb0836e 100644
--- a/Source/Microphysics/Null/NullMoist.H
+++ b/Source/Microphysics/Null/NullMoist.H
@@ -14,18 +14,17 @@ class NullMoist {
virtual
void
- define (SolverChoice& /*sc*/) { }
+ Define (SolverChoice& /*sc*/) { }
virtual
void Init (const amrex::MultiFab& /*cons_in*/,
- amrex::MultiFab& /*qmoist*/,
const amrex::BoxArray& /*grids*/,
const amrex::Geometry& /*geom*/,
const amrex::Real& /*dt_advance*/) { }
virtual
void
- Advance ( ) { }
+ Advance (const amrex::Real& /*dt_advance*/) { }
virtual
void
@@ -42,7 +41,7 @@ class NullMoist {
virtual
void
- Diagnose ( ) { }
+ Diagnose () { }
virtual
void
@@ -58,7 +57,7 @@ class NullMoist {
virtual
int
- Qmoist_Size ( ) { return NullMoist::m_qmoist_size; }
+ Qmoist_Size () { return NullMoist::m_qmoist_size; }
private:
int m_qmoist_size = 1;
diff --git a/Source/Microphysics/SAM/Diagnose_SAM.cpp b/Source/Microphysics/SAM/Diagnose_SAM.cpp
index 784f21a8e..544a202f6 100644
--- a/Source/Microphysics/SAM/Diagnose_SAM.cpp
+++ b/Source/Microphysics/SAM/Diagnose_SAM.cpp
@@ -1,4 +1,5 @@
#include "Microphysics.H"
+#include "EOS.H"
/**
* Computes diagnostic quantities like cloud ice/liquid and precipitation ice/liquid
@@ -14,6 +15,7 @@ void SAM::Diagnose () {
auto qci = mic_fab_vars[MicVar::qci];
auto qpl = mic_fab_vars[MicVar::qpl];
auto qpi = mic_fab_vars[MicVar::qpi];
+ auto qg = mic_fab_vars[MicVar::qg];
auto tabs = mic_fab_vars[MicVar::tabs];
for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
@@ -25,35 +27,55 @@ void SAM::Diagnose () {
auto qci_array = qci->array(mfi);
auto qpl_array = qpl->array(mfi);
auto qpi_array = qpi->array(mfi);
+ auto qg_array = qg->array(mfi);
auto tabs_array = tabs->array(mfi);
const auto& box3d = mfi.tilebox();
- amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
- qv_array(i,j,k) = qt_array(i,j,k) - qn_array(i,j,k);
- amrex::Real omn = std::max(0.0, std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg));
- qcl_array(i,j,k) = qn_array(i,j,k)*omn;
- qci_array(i,j,k) = qn_array(i,j,k)*(1.0-omn);
- amrex::Real omp = std::max(0.0, std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr));
- qpl_array(i,j,k) = qp_array(i,j,k)*omp;
- qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp);
+ amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+ {
+ qv_array(i,j,k) = qt_array(i,j,k) - qn_array(i,j,k);
+ amrex::Real omn = std::max(0.0, std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg));
+ qcl_array(i,j,k) = qn_array(i,j,k)*omn;
+ qci_array(i,j,k) = qn_array(i,j,k)*(1.0-omn);
+ amrex::Real omp = std::max(0.0, std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr));
+ qpl_array(i,j,k) = qp_array(i,j,k)*omp;
+ qpi_array(i,j,k) = qp_array(i,j,k)*(1.0-omp);
+
+ // TODO: If omp above is bound by [0, 1], shouldn't this always be 0?
+ // Graupel == precip total - rain - snow (but must be >= 0)
+ qg_array(i,j,k) = std::max(0.0, qp_array(i,j,k)-qpl_array(i,j,k)-qpi_array(i,j,k));
});
}
}
+
/**
- * Wrapper for the PrecipFall, Cloud, Precipitation, and Diagnostics routines.
+ * Computes diagnostic quantities like temperature and pressure
+ * from the existing Microphysics variables.
*/
-void SAM::Proc () {
+void SAM::Compute_Tabs_and_Pres ()
+{
+ auto rho = mic_fab_vars[MicVar::rho];
+ auto theta = mic_fab_vars[MicVar::theta];
+ auto qv = mic_fab_vars[MicVar::qv];
+ auto tabs = mic_fab_vars[MicVar::tabs];
+ auto pres = mic_fab_vars[MicVar::pres];
- MicroPrecipFall();
+ for ( amrex::MFIter mfi(*tabs, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+ const auto& box3d = mfi.tilebox();
- if (docloud) {
- Cloud();
- if (doprecip) {
- Precip();
- }
- }
+ auto rho_array = rho->array(mfi);
+ auto theta_array = theta->array(mfi);
+ auto qv_array = qv->array(mfi);
+ auto tabs_array = tabs->array(mfi);
+ auto pres_array = pres->array(mfi);
- Diagnose();
+ amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+ {
+ amrex::Real RhoTheta = rho_array(i,j,k)*theta_array(i,j,k);
+ tabs_array(i,j,k) = getTgivenRandRTh(rho_array(i,j,k), RhoTheta, qv_array(i,j,k));
+ pres_array(i,j,k) = getPgivenRTh(RhoTheta, qv_array(i,j,k))/100.;
+ });
+ }
}
diff --git a/Source/Microphysics/SAM/Init_SAM.cpp b/Source/Microphysics/SAM/Init_SAM.cpp
index ad5f0c484..19970f4f3 100644
--- a/Source/Microphysics/SAM/Init_SAM.cpp
+++ b/Source/Microphysics/SAM/Init_SAM.cpp
@@ -8,6 +8,240 @@
using namespace amrex;
+/**
+ * Initializes the Microphysics module.
+ *
+ * @param[in] cons_in Conserved variables input
+ * @param[in] qc_in Cloud variables input
+ * @param[in,out] qv_in Vapor variables input
+ * @param[in] qi_in Ice variables input
+ * @param[in] grids The boxes on which we will evolve the solution
+ * @param[in] geom Geometry associated with these MultiFabs and grids
+ * @param[in] dt_advance Timestep for the advance
+ */
+void SAM::Init (const MultiFab& cons_in,
+ const BoxArray& grids,
+ const Geometry& geom,
+ const Real& dt_advance)
+{
+ dt = dt_advance;
+ m_geom = geom;
+ m_gtoe = grids;
+
+ MicVarMap.resize(m_qmoist_size);
+ MicVarMap = {MicVar::qv, MicVar::qcl, MicVar::qci,
+ MicVar::qr, MicVar::qpi, MicVar::qg};
+
+ // initialize microphysics variables
+ for (auto ivar = 0; ivar < MicVar::NumVars; ++ivar) {
+ mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
+ 1, cons_in.nGrowVect());
+ mic_fab_vars[ivar]->setVal(0.);
+ }
+
+ // Set class data members
+ for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) {
+ const auto& box3d = mfi.tilebox();
+
+ const auto& lo = amrex::lbound(box3d);
+ const auto& hi = amrex::ubound(box3d);
+
+ nlev = box3d.length(2);
+ zlo = lo.z;
+ zhi = hi.z;
+
+ // parameters
+ accrrc.resize({zlo}, {zhi});
+ accrsi.resize({zlo}, {zhi});
+ accrsc.resize({zlo}, {zhi});
+ coefice.resize({zlo}, {zhi});
+ evaps1.resize({zlo}, {zhi});
+ evaps2.resize({zlo}, {zhi});
+ accrgi.resize({zlo}, {zhi});
+ accrgc.resize({zlo}, {zhi});
+ evapg1.resize({zlo}, {zhi});
+ evapg2.resize({zlo}, {zhi});
+ evapr1.resize({zlo}, {zhi});
+ evapr2.resize({zlo}, {zhi});
+
+ // data (input)
+ rho1d.resize({zlo}, {zhi});
+ pres1d.resize({zlo}, {zhi});
+ tabs1d.resize({zlo}, {zhi});
+ gamaz.resize({zlo}, {zhi});
+ zmid.resize({zlo}, {zhi});
+
+ // output
+ qifall.resize({zlo}, {zhi});
+ tlatqi.resize({zlo}, {zhi});
+ }
+}
+
+
+/**
+ * Initializes the Microphysics module.
+ *
+ * @param[in] cons_in Conserved variables input
+ */
+void SAM::Copy_State_to_Micro (const MultiFab& cons_in)
+{
+ // Get the temperature, density, theta, qt and qp from input
+ for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
+ const auto& box3d = mfi.tilebox();
+
+ auto states_array = cons_in.array(mfi);
+
+ auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi);
+ auto qp_array = mic_fab_vars[MicVar::qp]->array(mfi);
+ auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi);
+ auto qc_array = mic_fab_vars[MicVar::qcl]->array(mfi);
+ auto qi_array = mic_fab_vars[MicVar::qci]->array(mfi);
+ auto rho_array = mic_fab_vars[MicVar::rho]->array(mfi);
+ auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi);
+
+ // Get pressure, theta, temperature, density, and qt, qp
+ amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+ {
+ rho_array(i,j,k) = states_array(i,j,k,Rho_comp);
+ theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp);
+ qt_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp);
+ qp_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp);
+ qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k);
+ //qv_array(i,j,k) = qv_array_from_moist(i,j,k);
+ //temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k));
+ //pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/100.;
+ });
+ }
+}
+
+void SAM::Compute_Coefficients ()
+{
+ auto dz = m_geom.CellSize(2);
+ auto lowz = m_geom.ProbLo(2);
+
+ auto accrrc_t = accrrc.table();
+ auto accrsi_t = accrsi.table();
+ auto accrsc_t = accrsc.table();
+ auto coefice_t = coefice.table();
+ auto evaps1_t = evaps1.table();
+ auto evaps2_t = evaps2.table();
+ auto accrgi_t = accrgi.table();
+ auto accrgc_t = accrgc.table();
+ auto evapg1_t = evapg1.table();
+ auto evapg2_t = evapg2.table();
+ auto evapr1_t = evapr1.table();
+ auto evapr2_t = evapr2.table();
+
+ auto rho1d_t = rho1d.table();
+ auto pres1d_t = pres1d.table();
+ auto tabs1d_t = tabs1d.table();
+
+ auto gamaz_t = gamaz.table();
+ auto zmid_t = zmid.table();
+
+ Real gam3 = erf_gammafff(3.0 );
+ Real gamr1 = erf_gammafff(3.0+b_rain );
+ Real gamr2 = erf_gammafff((5.0+b_rain)/2.0);
+ // Real gamr3 = erf_gammafff(4.0+b_rain );
+ Real gams1 = erf_gammafff(3.0+b_snow );
+ Real gams2 = erf_gammafff((5.0+b_snow)/2.0);
+ // Real gams3 = erf_gammafff(4.0+b_snow );
+ Real gamg1 = erf_gammafff(3.0+b_grau );
+ Real gamg2 = erf_gammafff((5.0+b_grau)/2.0);
+ // Real gamg3 = erf_gammafff(4.0+b_grau );
+
+ // calculate the plane average variables
+ PlaneAverage rho_ave(mic_fab_vars[MicVar::rho].get(), m_geom, m_axis);
+ PlaneAverage theta_ave(mic_fab_vars[MicVar::theta].get(), m_geom, m_axis);
+ rho_ave.compute_averages(ZDir(), rho_ave.field());
+ theta_ave.compute_averages(ZDir(), theta_ave.field());
+
+ // get host variable rho, and rhotheta
+ int ncell = rho_ave.ncell_line();
+
+ Gpu::HostVector<Real> rho_h(ncell), theta_h(ncell);
+ rho_ave.line_average(0, rho_h);
+ theta_ave.line_average(0, theta_h);
+
+ // copy data to device
+ Gpu::DeviceVector<Real> rho_d(ncell), theta_d(ncell);
+ Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin());
+ Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin());
+ Gpu::streamSynchronize();
+
+ Real* rho_dptr = rho_d.data();
+ Real* theta_dptr = theta_d.data();
+
+ Real gOcp = m_gOcp;
+
+ // TODO: None of this is qv aware...?
+ amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept
+ {
+ Real RhoTheta = rho_dptr[k]*theta_dptr[k];
+ Real pressure = getPgivenRTh(RhoTheta);
+ rho1d_t(k) = rho_dptr[k];
+ pres1d_t(k) = pressure/100.;
+ tabs1d_t(k) = getTgivenRandRTh(rho_dptr[k],RhoTheta);
+ zmid_t(k) = lowz + (k+0.5)*dz;
+ gamaz_t(k) = gOcp*zmid_t(k);
+ });
+
+ if(round(gam3) != 2) {
+ std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl;
+ std::exit(-1);
+ }
+
+ // Populate all the coefficients
+ amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept
+ {
+ Real pratio = sqrt(1.29 / rho1d_t(k));
+ Real rrr1=393.0/(tabs1d_t(k)+120.0)*std::pow((tabs1d_t(k)/273.0),1.5);
+ Real rrr2=std::pow((tabs1d_t(k)/273.0),1.94)*(1000.0/pres1d_t(k));
+ Real estw = 100.0*erf_esatw(tabs1d_t(k));
+ Real esti = 100.0*erf_esati(tabs1d_t(k));
+
+ // accretion by snow:
+ Real coef1 = 0.25 * PI * nzeros * a_snow * gams1 * pratio/pow((PI * rhos * nzeros/rho1d_t(k) ) , ((3.0+b_snow)/4.0));
+ Real coef2 = exp(0.025*(tabs1d_t(k) - 273.15));
+ accrsi_t(k) = coef1 * coef2 * esicoef;
+ accrsc_t(k) = coef1 * esccoef;
+ coefice_t(k) = coef2;
+
+ // evaporation of snow:
+ coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k));
+ coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti);
+ evaps1_t(k) = 0.65*4.0*nzeros/sqrt(PI*rhos*nzeros)/(coef1+coef2)/sqrt(rho1d_t(k));
+ evaps2_t(k) = 0.49*4.0*nzeros*gams2*sqrt(a_snow/(muelq*rrr1))/pow((PI*rhos*nzeros) , ((5.0+b_snow)/8.0)) /
+ (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_snow)/8.0))*sqrt(pratio);
+
+ // accretion by graupel:
+ coef1 = 0.25*PI*nzerog*a_grau*gamg1*pratio/pow((PI*rhog*nzerog/rho1d_t(k)) , ((3.0+b_grau)/4.0));
+ coef2 = exp(0.025*(tabs1d_t(k) - 273.15));
+ accrgi_t(k) = coef1 * coef2 * egicoef;
+ accrgc_t(k) = coef1 * egccoef;
+
+ // evaporation of graupel:
+ coef1 =(lsub/(tabs1d_t(k)*rv)-1.0)*lsub/(therco*rrr1*tabs1d_t(k));
+ coef2 = rv*tabs1d_t(k)/(diffelq*rrr2*esti);
+ evapg1_t(k) = 0.65*4.0*nzerog/sqrt(PI*rhog*nzerog)/(coef1+coef2)/sqrt(rho1d_t(k));
+ evapg2_t(k) = 0.49*4.0*nzerog*gamg2*sqrt(a_grau/(muelq*rrr1))/pow((PI * rhog * nzerog) , ((5.0+b_grau)/8.0)) /
+ (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_grau)/8.0))*sqrt(pratio);
+
+ // accretion by rain:
+ accrrc_t(k)= 0.25 * PI * nzeror * a_rain * gamr1 * pratio/pow((PI * rhor * nzeror / rho1d_t(k)) , ((3+b_rain)/4.))* erccoef;
+
+ // evaporation of rain:
+ coef1 =(lcond/(tabs1d_t(k)*rv)-1.0)*lcond/(therco*rrr1*tabs1d_t(k));
+ coef2 = rv*tabs1d_t(k)/(diffelq * rrr2 * estw);
+ evapr1_t(k) = 0.78 * 2.0 * PI * nzeror /
+ sqrt(PI * rhor * nzeror) / (coef1+coef2) / sqrt(rho1d_t(k));
+ evapr2_t(k) = 0.31 * 2.0 * PI * nzeror * gamr2 * 0.89 * sqrt(a_rain/(muelq*rrr1))/
+ pow((PI * rhor * nzeror) , ((5.0+b_rain)/8.0)) /
+ (coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio);
+ });
+}
+
+#if 0
/**
* Initializes the Microphysics module.
*
@@ -243,3 +477,4 @@ void SAM::Init (const MultiFab& cons_in, MultiFab& qmoist,
(coef1+coef2) * pow(rho1d_t(k) , ((1.0+b_rain)/8.0))*sqrt(pratio);
});
}
+#endif
diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H
index 31f63ee1c..4cccf0995 100644
--- a/Source/Microphysics/SAM/SAM.H
+++ b/Source/Microphysics/SAM/SAM.H
@@ -40,6 +40,7 @@ namespace MicVar {
qcl, // cloud water
qpl, // precip rain
qpi, // precip ice
+ qg, // graupel
// temporary variable
omega,
NumVars
@@ -75,15 +76,12 @@ public:
// micro interface for precip fall
void MicroPrecipFall ();
- // process microphysics
- void Proc ();
-
// diagnose
void Diagnose () override;
// Set up for first time
void
- define (SolverChoice& sc) override
+ Define (SolverChoice& sc) override
{
docloud = sc.do_cloud;
doprecip = sc.do_precip;
@@ -97,18 +95,17 @@ public:
// init
void
Init (const amrex::MultiFab& cons_in,
- amrex::MultiFab& qmoist,
const amrex::BoxArray& grids,
const amrex::Geometry& geom,
const amrex::Real& dt_advance) override;
// Copy state into micro vars
void
- Copy_State_to_Micro (const amrex::MultiFab& /*cons_in*/) override { }
+ Copy_State_to_Micro (const amrex::MultiFab& cons_in) override;
// Copy state into micro vars
void
- Copy_Micro_to_State (amrex::MultiFab& /*cons_in*/) override { }
+ Copy_Micro_to_State (amrex::MultiFab& cons_in) override;
// update ERF variables
void
@@ -120,7 +117,8 @@ public:
{
this->Copy_State_to_Micro(cons_in);
this->Diagnose();
- //this->Compute_Tabs_and_Pres();
+ this->Compute_Tabs_and_Pres();
+ this->Compute_Coefficients();
}
void
@@ -131,8 +129,10 @@ public:
// wrapper to do all the updating
void
- Advance ( ) override
+ Advance (const amrex::Real& dt_advance) override
{
+ dt = dt_advance;
+
this->Cloud();
this->Diagnose();
this->IceFall();
@@ -143,19 +143,32 @@ public:
amrex::MultiFab*
Qmoist_Ptr (const int& varIdx) override
{
- AMREX_ALWAYS_ASSERT(varIdx <= m_qmoist_size);
- return mic_fab_vars[varIdx].get();
+ AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
+ return mic_fab_vars[MicVarMap[varIdx]].get();
}
+ void
+ Compute_Tabs_and_Pres ( );
+
+ void
+ Compute_Coefficients ( );
+
int
Qmoist_Size ( ) override { return SAM::m_qmoist_size; }
- private:
+ // independent variables
+ amrex::Array<FabPtr, MicVar::NumVars> mic_fab_vars;
+
+private:
// Number of qmoist variables
int m_qmoist_size = 6;
+ // MicVar map (Qmoist indices -> MicVar enum)
+ amrex::Vector<int> MicVarMap;
+
// geometry
amrex::Geometry m_geom;
+
// valid boxes on which to evolve the solution
amrex::BoxArray m_gtoe;
@@ -199,9 +212,6 @@ public:
amrex::TableData<amrex::Real, 1> qv1d;
amrex::TableData<amrex::Real, 1> qn1d;
- // independent variables
- amrex::Array<FabPtr, MicVar::NumVars> mic_fab_vars;
-
amrex::TableData<amrex::Real, 1> gamaz;
amrex::TableData<amrex::Real, 1> zmid; // mid value of vertical coordinate in physical domain
diff --git a/Source/Microphysics/SAM/Update_SAM.cpp b/Source/Microphysics/SAM/Update_SAM.cpp
index 7e391d1dd..bd7531d7f 100644
--- a/Source/Microphysics/SAM/Update_SAM.cpp
+++ b/Source/Microphysics/SAM/Update_SAM.cpp
@@ -57,4 +57,38 @@ void SAM::Update (amrex::MultiFab& cons,
qmoist.FillBoundary(m_geom.periodicity());
}
+/**
+ * Updates conserved and microphysics variables in the provided MultiFabs from
+ * the internal MultiFabs that store Microphysics module data.
+ *
+ * @param[out] cons Conserved variables
+ * @param[out] qmoist: qv, qc, qi, qr, qs, qg
+ */
+void SAM::Copy_Micro_to_State (amrex::MultiFab& cons)
+{
+ // Get the temperature, density, theta, qt and qp from input
+ for ( amrex::MFIter mfi(cons,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+ const auto& box3d = mfi.tilebox();
+
+ auto states_arr = cons.array(mfi);
+
+ auto rho_arr = mic_fab_vars[MicVar_Kess::rho]->array(mfi);
+ auto theta_arr = mic_fab_vars[MicVar_Kess::theta]->array(mfi);
+ auto qt_arr = mic_fab_vars[MicVar_Kess::qt]->array(mfi);
+ auto qp_arr = mic_fab_vars[MicVar_Kess::qp]->array(mfi);
+
+ // get potential total density, temperature, qt, qp
+ amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+ {
+ //states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k);
+ states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k);
+ states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k);
+ states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qp_arr(i,j,k);
+ });
+ }
+
+ // Fill interior ghost cells and periodic boundaries
+ cons.FillBoundary(m_geom.periodicity());
+}
+
diff --git a/Source/Radiation/Radiation.H b/Source/Radiation/Radiation.H
index 0eb4e1beb..0c072ff5f 100644
--- a/Source/Radiation/Radiation.H
+++ b/Source/Radiation/Radiation.H
@@ -34,56 +34,57 @@
// Radiation code interface class
class Radiation {
public:
- Radiation() {
+ Radiation () {
// First, make sure yakl has been initialized
if (!yakl::isInitialized()) yakl::init();
}
- ~Radiation() = default;
+ ~Radiation () = default;
// init
- void initialize(const amrex::MultiFab& cons_in, amrex::MultiFab& qmoist,
- const amrex::BoxArray& grids,
- const amrex::Geometry& geom,
- const amrex::Real& dt_advance,
- const bool& do_sw_rad,
- const bool& do_lw_rad,
- const bool& do_aero_rad,
- const bool& do_snow_opt,
- const bool& is_cmip6_volcano);
+ void initialize (const amrex::MultiFab& cons_in,
+ amrex::MultiFab& qmoist,
+ const amrex::BoxArray& grids,
+ const amrex::Geometry& geom,
+ const amrex::Real& dt_advance,
+ const bool& do_sw_rad,
+ const bool& do_lw_rad,
+ const bool& do_aero_rad,
+ const bool& do_snow_opt,
+ const bool& is_cmip6_volcano);
// run radiation model
- void run();
+ void run ();
// call back
- void on_complete();
-
- void radiation_driver_lw(int ncol, int nlev,
- const real3d& gas_vmr,
- const real2d& pmid, const real2d& pint, const real2d& tmid, const real2d& tint,
- const real3d& cld_tau_gpt, const real3d& aer_tau_bnd,
- FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
- const real2d& qrl, const real2d& qrlc);
-
- void radiation_driver_sw(int ncol,
- const real3d& gas_vmr, const real2d& pmid, const real2d& pint, const real2d& tmid,
- const real2d& albedo_dir, const real2d& albedo_dif, const real1d& coszrs,
- const real3d& cld_tau_gpt, const real3d& cld_ssa_gpt, const real3d& cld_asm_gpt,
- const real3d& aer_tau_bnd, const real3d& aer_ssa_bnd, const real3d& aer_asm_bnd,
- FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
+ void on_complete ();
+
+ void radiation_driver_lw (int ncol, int nlev,
+ const real3d& gas_vmr,
+ const real2d& pmid, const real2d& pint, const real2d& tmid, const real2d& tint,
+ const real3d& cld_tau_gpt, const real3d& aer_tau_bnd,
+ FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
+ const real2d& qrl, const real2d& qrlc);
+
+ void radiation_driver_sw (int ncol,
+ const real3d& gas_vmr, const real2d& pmid, const real2d& pint, const real2d& tmid,
+ const real2d& albedo_dir, const real2d& albedo_dif, const real1d& coszrs,
+ const real3d& cld_tau_gpt, const real3d& cld_ssa_gpt, const real3d& cld_asm_gpt,
+ const real3d& aer_tau_bnd, const real3d& aer_ssa_bnd, const real3d& aer_asm_bnd,
+ FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
const real2d& qrs, const real2d& qrsc);
- void set_daynight_indices(const real1d& coszrs,
- const int1d& day_indices,
- const int1d& night_indices);
+ void set_daynight_indices (const real1d& coszrs,
+ const int1d& day_indices,
+ const int1d& night_indices);
- void get_gas_vmr(const std::vector<std::string>& gas_names,
- const real3d& gas_vmr);
+ void get_gas_vmr (const std::vector<std::string>& gas_names,
+ const real3d& gas_vmr);
- void calculate_heating_rate(const real2d& flux_up,
- const real2d& flux_dn,
- const real2d& pint,
- const real2d& heating_rate);
+ void calculate_heating_rate (const real2d& flux_up,
+ const real2d& flux_dn,
+ const real2d& pint,
+ const real2d& heating_rate);
private:
// geometry
amrex::Geometry m_geom;
diff --git a/Source/Radiation/Radiation.cpp b/Source/Radiation/Radiation.cpp
index 8767ea9a2..56bee93c3 100644
--- a/Source/Radiation/Radiation.cpp
+++ b/Source/Radiation/Radiation.cpp
@@ -74,7 +74,8 @@ namespace internal {
}
// init
-void Radiation::initialize(const MultiFab& cons_in, MultiFab& qmoist,
+void Radiation::initialize(const MultiFab& cons_in,
+ MultiFab& qmoist,
const BoxArray& grids,
const Geometry& geom,
const Real& dt_advance,
diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp
index 362464d16..b7f987cbb 100644
--- a/Source/TimeIntegration/ERF_Advance.cpp
+++ b/Source/TimeIntegration/ERF_Advance.cpp
@@ -56,7 +56,8 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/
&vars_old[lev][Vars::yvel], &vars_old[lev][Vars::zvel]});
if (solverChoice.moisture_type != MoistureType::None) {
- FillPatchMoistVars(lev, *qmoist[lev]);
+ // TODO: This is only qv
+ FillPatchMoistVars(lev, *(qmoist[lev][0]));
}
MultiFab* S_crse;
diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp
index 3ed147cc6..101918fd9 100644
--- a/Source/TimeIntegration/ERF_advance_dycore.cpp
+++ b/Source/TimeIntegration/ERF_advance_dycore.cpp
@@ -256,7 +256,7 @@ void ERF::advance_dycore(int level,
Gpu::DeviceVector<Real> qv_d(ncell), qc_d(ncell), qi_d(ncell);
if (q_size >= 1) {
- MultiFab qvapor (*(qmoist[level]), make_alias, 0, 1);
+ MultiFab qvapor (*(qmoist[level][0]), make_alias, 0, 1);
PlaneAverage qv_ave(&qvapor, geom[level], solverChoice.ave_plane);
qv_ave.compute_averages(ZDir(), qv_ave.field());
@@ -267,7 +267,7 @@ void ERF::advance_dycore(int level,
}
if (q_size >= 2) {
- MultiFab qcloud (*(qmoist[level]), make_alias, 1, 1);
+ MultiFab qcloud (*(qmoist[level][1]), make_alias, 0, 1);
PlaneAverage qc_ave(&qcloud, geom[level], solverChoice.ave_plane);
qc_ave.compute_averages(ZDir(), qc_ave.field());
@@ -278,7 +278,7 @@ void ERF::advance_dycore(int level,
}
if (q_size >= 3) {
- MultiFab qice (*(qmoist[level]), make_alias, 2, 1);
+ MultiFab qice (*(qmoist[level][2]), make_alias, 0, 1);
PlaneAverage qi_ave(&qice , geom[level], solverChoice.ave_plane);
qi_ave.compute_averages(ZDir(), qi_ave.field());
diff --git a/Source/TimeIntegration/ERF_advance_microphysics.cpp b/Source/TimeIntegration/ERF_advance_microphysics.cpp
index fedc99b46..e726bab88 100644
--- a/Source/TimeIntegration/ERF_advance_microphysics.cpp
+++ b/Source/TimeIntegration/ERF_advance_microphysics.cpp
@@ -8,7 +8,7 @@ void ERF::advance_microphysics (int lev,
{
if (solverChoice.moisture_type != MoistureType::None) {
micro.Update_Micro_Vars_Lev(lev, cons);
- micro.Advance(lev);
+ micro.Advance(lev, dt_advance);
micro.Update_State_Vars_Lev(lev, cons);
}
}
diff --git a/Source/TimeIntegration/ERF_advance_radiation.cpp b/Source/TimeIntegration/ERF_advance_radiation.cpp
index c42222b60..826c0d934 100644
--- a/Source/TimeIntegration/ERF_advance_radiation.cpp
+++ b/Source/TimeIntegration/ERF_advance_radiation.cpp
@@ -14,7 +14,8 @@ void ERF::advance_radiation (int lev,
bool do_snow_opt {true};
bool is_cmip6_volcano {true};
- rad.initialize(cons, qmoist[lev],
+ // TODO: Only passing qv from the qmoist vector!!!
+ rad.initialize(cons, *(qmoist[lev][0]),
grids[lev],
Geom(lev),
dt_advance,
diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp
index 22c540c4c..c6dea21cb 100644
--- a/Source/TimeIntegration/ERF_make_buoyancy.cpp
+++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp
@@ -36,7 +36,7 @@ using namespace amrex;
void make_buoyancy (Vector<MultiFab>& S_data,
const MultiFab& S_prim,
MultiFab& buoyancy,
- MultiFab* qmoist,
+ Vector<MultiFab*> qmoist,
Gpu::DeviceVector<Real> qv_d,
Gpu::DeviceVector<Real> qc_d,
Gpu::DeviceVector<Real> qi_d,
@@ -184,9 +184,9 @@ void make_buoyancy (Vector<MultiFab>& S_data,
// Base state density
const Array4<const Real>& r0_arr = r0->const_array(mfi);
- const Array4<const Real> & qv_data = qmoist->array(mfi,0);
- const Array4<const Real> & qc_data = qmoist->array(mfi,1);
- const Array4<const Real> & qi_data = qmoist->array(mfi,2);
+ const Array4<const Real> & qv_data = qmoist[0]->array(mfi);
+ const Array4<const Real> & qc_data = qmoist[1]->array(mfi);
+ const Array4<const Real> & qi_data = qmoist[2]->array(mfi);
amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
@@ -252,9 +252,9 @@ void make_buoyancy (Vector<MultiFab>& S_data,
const Array4<const Real> & cell_data = S_data[IntVar::cons].array(mfi);
const Array4<const Real> & cell_prim = S_prim.array(mfi);
- const Array4<const Real> & qv_data = qmoist->const_array(mfi,0);
- const Array4<const Real> & qc_data = qmoist->const_array(mfi,1);
- const Array4<const Real> & qi_data = qmoist->const_array(mfi,2);
+ const Array4<const Real> & qv_data = qmoist[0]->const_array(mfi);
+ const Array4<const Real> & qc_data = qmoist[1]->const_array(mfi);
+ const Array4<const Real> & qi_data = qmoist[2]->const_array(mfi);
amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
@@ -318,9 +318,9 @@ void make_buoyancy (Vector<MultiFab>& S_data,
const Array4<const Real> & cell_data = S_data[IntVar::cons].array(mfi);
const Array4<const Real> & cell_prim = S_prim.array(mfi);
- const Array4<const Real> & qv_data = qmoist->const_array(mfi,0);
- const Array4<const Real> & qc_data = qmoist->const_array(mfi,1);
- const Array4<const Real> & qi_data = qmoist->const_array(mfi,2);
+ const Array4<const Real> & qv_data = qmoist[0]->const_array(mfi);
+ const Array4<const Real> & qc_data = qmoist[1]->const_array(mfi);
+ const Array4<const Real> & qi_data = qmoist[2]->const_array(mfi);
amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp
index 4b673ac08..0a0e042eb 100644
--- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp
+++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp
@@ -80,7 +80,7 @@ void erf_slow_rhs_pre (int level, int finest_level,
const MultiFab& xvel,
const MultiFab& yvel,
const MultiFab& zvel,
- const MultiFab* qv,
+ Vector<MultiFab*> qmoist,
std::unique_ptr<MultiFab>& z_t_mf,
MultiFab& Omega,
const MultiFab& source,
@@ -543,7 +543,7 @@ void erf_slow_rhs_pre (int level, int finest_level,
FArrayBox pprime; pprime.resize(gbx,1);
Elixir pp_eli = pprime.elixir();
const Array4<Real > & pp_arr = pprime.array();
- const Array4<Real const> & qv_arr = (use_moisture) ? qv->const_array(mfi) : Array4<Real> {};
+ const Array4<Real const> & qv_arr = (use_moisture) ? qmoist[0]->const_array(mfi) : Array4<Real> {};
{
BL_PROFILE("slow_rhs_pre_pprime");
ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
diff --git a/Source/TimeIntegration/TI_headers.H b/Source/TimeIntegration/TI_headers.H
index 01e72c9b4..6e72b11d4 100644
--- a/Source/TimeIntegration/TI_headers.H
+++ b/Source/TimeIntegration/TI_headers.H
@@ -22,7 +22,7 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk,
const amrex::MultiFab& xvel,
const amrex::MultiFab& yvel,
const amrex::MultiFab& zvel,
- const amrex::MultiFab* qmoist,
+ amrex::Vector<amrex::MultiFab*> qmoist,
std::unique_ptr<amrex::MultiFab>& z_t,
amrex::MultiFab& Omega,
const amrex::MultiFab& source,
@@ -218,7 +218,7 @@ void make_fast_coeffs (int level,
void make_buoyancy (amrex::Vector< amrex::MultiFab>& S_data,
const amrex::MultiFab & S_prim,
amrex::MultiFab& buoyancy,
- amrex::MultiFab* qmoist,
+ amrex::Vector<amrex::MultiFab*> qmoist,
const amrex::Gpu::DeviceVector<amrex::Real> qv_d,
const amrex::Gpu::DeviceVector<amrex::Real> qc_d,
const amrex::Gpu::DeviceVector<amrex::Real> qi_d,
diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H
index 5ac193aa0..f783597e5 100644
--- a/Source/TimeIntegration/TI_slow_rhs_fun.H
+++ b/Source/TimeIntegration/TI_slow_rhs_fun.H
@@ -90,12 +90,12 @@
MultiFab* p0_new = &p_hse_new;
make_buoyancy(S_data, S_prim, buoyancy,
- qmoist[level].get(), qv_d, qc_d, qi_d,
+ qmoist[level], qv_d, qc_d, qi_d,
fine_geom, solverChoice, r0_new);
erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch,
xvel_new, yvel_new, zvel_new,
- qmoist[level].get(),
+ qmoist[level],
z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12,
Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs,
Hfx3, Diss,
@@ -183,12 +183,12 @@
// If not moving_terrain
make_buoyancy(S_data, S_prim, buoyancy,
- qmoist[level].get(), qv_d, qc_d, qi_d,
+ qmoist[level], qv_d, qc_d, qi_d,
fine_geom, solverChoice, r0);
erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_data, S_prim, S_scratch,
xvel_new, yvel_new, zvel_new,
- qmoist[level].get(),
+ qmoist[level],
z_t_rk[level], Omega, source, buoyancy, Tau11, Tau22, Tau33, Tau12,
Tau13, Tau21, Tau23, Tau31, Tau32, SmnSmn, eddyDiffs,
Hfx3, Diss,