ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Public Member Functions | Static Public Member Functions | Private Types | Private Attributes | Static Private Attributes | List of all members
SAM Class Reference

#include <ERF_SAM.H>

Inheritance diagram for SAM:
Inheritance graph
[legend]
Collaboration diagram for SAM:
Collaboration graph
[legend]

Public Member Functions

 SAM ()
 
virtual ~SAM ()=default
 
void Cloud (const SolverChoice &sc)
 
void IceFall (const SolverChoice &sc)
 
void Precip (const SolverChoice &sc)
 
void PrecipFall (const SolverChoice &sc)
 
void Define (SolverChoice &sc) override
 
void Init (const amrex::MultiFab &cons_in, const amrex::BoxArray &grids, const amrex::Geometry &geom, const amrex::Real &dt_advance, std::unique_ptr< amrex::MultiFab > &z_phys_nd, std::unique_ptr< amrex::MultiFab > &detJ_cc) override
 
void Set_dzmin (const amrex::Real dz_min) override
 
void Copy_State_to_Micro (const amrex::MultiFab &cons_in) override
 
void Copy_Micro_to_State (amrex::MultiFab &cons_in) override
 
void Update_Micro_Vars (amrex::MultiFab &cons_in) override
 
void Update_State_Vars (amrex::MultiFab &cons_in) override
 
void Advance (const amrex::Real &dt_advance, const SolverChoice &sc) override
 
amrex::MultiFab * Qmoist_Ptr (const int &varIdx) override
 
void Compute_Coefficients ()
 
int Qmoist_Size () override
 
int Qstate_Moist_Size () override
 
void Set_RealWidth (const int real_width) override
 
void Qmoist_Restart_Vars (const SolverChoice &, std::vector< int > &a_idx, std::vector< std::string > &a_names) const override
 
- Public Member Functions inherited from NullMoist
 NullMoist ()
 
virtual ~NullMoist ()=default
 
virtual int Qstate_NonMoist_Size ()
 
virtual void GetPlotVarNames (amrex::Vector< std::string > &a_vec) const
 
virtual void GetPlotVar (const std::string &, amrex::MultiFab &) const
 

Static Public Member Functions

AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real NewtonIterSat (int &i, int &j, int &k, const int &SAM_moisture_type, const amrex::Real &fac_cond, const amrex::Real &fac_fus, const amrex::Real &fac_sub, const amrex::Real &an, const amrex::Real &bn, const amrex::Array4< amrex::Real > &tabs_array, const amrex::Array4< amrex::Real > &pres_array, const amrex::Array4< amrex::Real > &qv_array, const amrex::Array4< amrex::Real > &qc_array, const amrex::Array4< amrex::Real > &qi_array, const amrex::Array4< amrex::Real > &qn_array, const amrex::Array4< amrex::Real > &qt_array)
 

Private Types

using FabPtr = std::shared_ptr< amrex::MultiFab >
 

Private Attributes

int m_qmoist_size = 3
 
int n_qstate_moist_size = 6
 
amrex::Vector< int > MicVarMap
 
amrex::Geometry m_geom
 
int m_real_width {0}
 
amrex::Real dt
 
int nlev
 
int zlo
 
int zhi
 
int m_axis
 
amrex::Real m_fac_cond
 
amrex::Real m_fac_fus
 
amrex::Real m_fac_sub
 
amrex::Real m_rdOcp
 
bool m_do_cond
 
amrex::Real m_dzmin
 
amrex::MultiFab * m_z_phys_nd
 
amrex::MultiFab * m_detJ_cc
 
amrex::Array< FabPtr, MicVar::NumVars > mic_fab_vars
 
amrex::TableData< amrex::Real, 1 > accrrc
 
amrex::TableData< amrex::Real, 1 > accrsi
 
amrex::TableData< amrex::Real, 1 > accrsc
 
amrex::TableData< amrex::Real, 1 > coefice
 
amrex::TableData< amrex::Real, 1 > evaps1
 
amrex::TableData< amrex::Real, 1 > evaps2
 
amrex::TableData< amrex::Real, 1 > accrgi
 
amrex::TableData< amrex::Real, 1 > accrgc
 
amrex::TableData< amrex::Real, 1 > evapg1
 
amrex::TableData< amrex::Real, 1 > evapg2
 
amrex::TableData< amrex::Real, 1 > evapr1
 
amrex::TableData< amrex::Real, 1 > evapr2
 
amrex::TableData< amrex::Real, 1 > rho1d
 
amrex::TableData< amrex::Real, 1 > pres1d
 
amrex::TableData< amrex::Real, 1 > tabs1d
 
amrex::TableData< amrex::Real, 1 > qt1d
 
amrex::TableData< amrex::Real, 1 > qv1d
 
amrex::TableData< amrex::Real, 1 > qn1d
 

Static Private Attributes

static constexpr amrex::Real CFL_MAX = myhalf
 

Member Typedef Documentation

◆ FabPtr

using SAM::FabPtr = std::shared_ptr<amrex::MultiFab>
private

Constructor & Destructor Documentation

◆ SAM()

SAM::SAM ( )
inline
59 {}

◆ ~SAM()

virtual SAM::~SAM ( )
virtualdefault

Member Function Documentation

◆ Advance()

void SAM::Advance ( const amrex::Real &  dt_advance,
const SolverChoice &  sc 
)
inlineoverridevirtual

Reimplemented from NullMoist.

129  {
130  dt = dt_advance;
131 
132  this->Cloud(sc);
133  this->IceFall(sc);
134  this->Precip(sc);
135  this->PrecipFall(sc);
136  }
SAM::Precip
void Precip(const SolverChoice &sc)
Definition: ERF_Precip.cpp:10
SAM::IceFall
void IceFall(const SolverChoice &sc)
Definition: ERF_IceFall.cpp:11
SAM::Cloud
void Cloud(const SolverChoice &sc)
Definition: ERF_CloudSAM.cpp:12
SAM::PrecipFall
void PrecipFall(const SolverChoice &sc)
Definition: ERF_PrecipFall.cpp:16
SAM::dt
amrex::Real dt
Definition: ERF_SAM.H:307
Here is the call graph for this function:

◆ Cloud()

void SAM::Cloud ( const SolverChoice &  sc)

Split cloud components according to saturation pressures; source theta from latent heat.

13 {
14  if (!m_do_cond) { return; }
15 
16  constexpr Real an = one/(tbgmax-tbgmin);
17  constexpr Real bn = tbgmin*an;
18 
19  Real fac_cond = m_fac_cond;
20  Real fac_sub = m_fac_sub;
21  Real fac_fus = m_fac_fus;
22  Real rdOcp = m_rdOcp;
23 
24  int SAM_moisture_type = 1;
25  if (sc.moisture_type == MoistureType::SAM_NoIce ||
26  sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
27  SAM_moisture_type = 2;
28  }
29 
30  auto domain = m_geom.Domain();
31  int i_lo = domain.smallEnd(0);
32  int i_hi = domain.bigEnd(0);
33  int j_lo = domain.smallEnd(1);
34  int j_hi = domain.bigEnd(1);
35 
36  for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]), TilingIfNotGPU()); mfi.isValid(); ++mfi) {
37  auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi);
38  auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi);
39  auto qv_array = mic_fab_vars[MicVar::qv]->array(mfi);
40  auto qcl_array = mic_fab_vars[MicVar::qcl]->array(mfi);
41  auto qci_array = mic_fab_vars[MicVar::qci]->array(mfi);
42 
43  auto tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi);
44  auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi);
45  auto pres_array = mic_fab_vars[MicVar::pres]->array(mfi);
46 
47  auto tbx = mfi.tilebox();
48  if (tbx.smallEnd(0) == i_lo) { tbx.growLo(0,-m_real_width); }
49  if (tbx.bigEnd(0) == i_hi) { tbx.growHi(0,-m_real_width); }
50  if (tbx.smallEnd(1) == j_lo) { tbx.growLo(1,-m_real_width); }
51  if (tbx.bigEnd(1) == j_hi) { tbx.growHi(1,-m_real_width); }
52 
53  ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
54  {
55  // Saturation moisture fractions
56  Real omn;
57  Real qsat;
58  Real qsatw;
59  Real qsati;
60 
61  // Newton iteration vars
62  Real delta_qv, delta_qc, delta_qi;
63 
64  // NOTE: Conversion before iterations is necessary to
65  // convert cloud water to ice or vice versa.
66  // This ensures the omn splitting is enforced
67  // before the Newton iteration, which assumes it is.
68 
69  omn = one;
70  if (SAM_moisture_type == 1){
71  // Cloud ice not permitted (melt to form water)
72  if (tabs_array(i,j,k) >= tbgmax) {
73  omn = one;
74  delta_qi = qci_array(i,j,k);
75  qci_array(i,j,k) = zero;
76  qcl_array(i,j,k) += delta_qi;
77  tabs_array(i,j,k) -= fac_fus * delta_qi;
78  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
79  }
80  // Cloud water not permitted (freeze to form ice)
81  else if (tabs_array(i,j,k) <= tbgmin) {
82  omn = zero;
83  delta_qc = qcl_array(i,j,k);
84  qcl_array(i,j,k) = zero;
85  qci_array(i,j,k) += delta_qc;
86  tabs_array(i,j,k) += fac_fus * delta_qc;
87  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
88  }
89  // Mixed cloud phase (split according to omn)
90  else {
91  omn = an*tabs_array(i,j,k)-bn;
92  delta_qc = qcl_array(i,j,k) - qn_array(i,j,k) * omn;
93  delta_qi = qci_array(i,j,k) - qn_array(i,j,k) * (one - omn);
94  qcl_array(i,j,k) = qn_array(i,j,k) * omn;
95  qci_array(i,j,k) = qn_array(i,j,k) * (one - omn);
96  tabs_array(i,j,k) += fac_fus * delta_qc;
97  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
98  }
99  }
100  else if (SAM_moisture_type == 2)
101  {
102  // No ice. ie omn = one
103  delta_qc = qcl_array(i,j,k) - qn_array(i,j,k);
104  delta_qi = zero;
105  qcl_array(i,j,k) = qn_array(i,j,k);
106  qci_array(i,j,k) = zero;
107  tabs_array(i,j,k) += fac_cond * delta_qc;
108  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
109  }
110 
111  // Saturation moisture fractions
112  erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw);
113  erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati);
114  qsat = omn * qsatw + (one-omn) * qsati;
115 
116  // We have enough total moisture to relax to equilibrium
117  if (qt_array(i,j,k) > qsat) {
118 
119  // Update temperature
120  tabs_array(i,j,k) = NewtonIterSat(i, j, k , SAM_moisture_type ,
121  fac_cond , fac_fus , fac_sub ,
122  an , bn ,
123  tabs_array, pres_array,
124  qv_array , qcl_array , qci_array,
125  qn_array , qt_array);
126 
127  // Update theta
128  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), Real(100.0)*pres_array(i,j,k), rdOcp);
129 
130  //
131  // We cannot blindly relax to qsat, but we can convert qc/qi -> qv.
132  // The concept here is that if we put all the moisture into qv and modify
133  // the temperature, we can then check if qv > qsat occurs (for final T/P/qv).
134  // If the reduction in T/qsat and increase in qv does trigger the
135  // aforementioned condition, we can do Newton iteration to drive qv = qsat.
136  //
137  } else {
138  // Changes in each component
139  delta_qv = qcl_array(i,j,k) + qci_array(i,j,k);
140  delta_qc = qcl_array(i,j,k);
141  delta_qi = qci_array(i,j,k);
142 
143  // Partition the change in non-precipitating q
144  qv_array(i,j,k) += delta_qv;
145  qcl_array(i,j,k) = zero;
146  qci_array(i,j,k) = zero;
147  qn_array(i,j,k) = zero;
148  qt_array(i,j,k) = qv_array(i,j,k);
149 
150  // Update temperature (endothermic since we evap/sublime)
151  tabs_array(i,j,k) -= fac_cond * delta_qc + fac_sub * delta_qi;
152 
153  // Update theta
154  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), Real(100.0)*pres_array(i,j,k), rdOcp);
155 
156  // Verify assumption that qv > qsat does not occur
157  erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw);
158  erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati);
159  qsat = omn * qsatw + (one-omn) * qsati;
160  if (qt_array(i,j,k) > qsat) {
161 
162  // Update temperature
163  tabs_array(i,j,k) = NewtonIterSat(i, j, k , SAM_moisture_type ,
164  fac_cond , fac_fus , fac_sub ,
165  an , bn ,
166  tabs_array, pres_array,
167  qv_array , qcl_array , qci_array,
168  qn_array , qt_array);
169 
170  // Update theta
171  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), Real(100.0)*pres_array(i,j,k), rdOcp);
172 
173  }
174  }
175  });
176  } // mfi
177 }
tbgmax
constexpr amrex::Real tbgmax
Definition: ERF_Constants.H:42
one
constexpr amrex::Real one
Definition: ERF_Constants.H:7
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
tbgmin
constexpr amrex::Real tbgmin
Definition: ERF_Constants.H:41
getThgivenTandP
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getThgivenTandP(const amrex::Real T, const amrex::Real P, const amrex::Real rdOcp)
Definition: ERF_EOS.H:18
rdOcp
const Real rdOcp
Definition: ERF_InitCustomPert_Bomex.H:16
ParallelFor
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const auto prob_lo=geomdata.ProbLo();const auto dx=geomdata.CellSize();const Real x=(prob_lo[0]+(i+myhalf) *dx[0])/mf_m(i, j, 0);const Real z=z_cc(i, j, k);Real L=std::sqrt(std::pow((x - x_c)/x_r, 2)+std::pow((z - z_c)/z_r, 2));if(L<=one) { Real dT=T_pert *(std::cos(PI *L)+one)/two;Real Tbar_hse=p_hse(i, j, k)/(R_d *r_hse(i, j, k));Real theta_perturbed=(Tbar_hse+dT) *std::pow(p_0/p_hse(i, j, k), rdOcp);Real theta_0=(Tbar_hse) *std::pow(p_0/p_hse(i, j, k), rdOcp);if(const_rho) { state_pert(i, j, k, RhoTheta_comp)=r_hse(i, j, k) *(theta_perturbed - theta_0);} else { state_pert(i, j, k, Rho_comp)=getRhoThetagivenP(p_hse(i, j, k))/theta_perturbed - r_hse(i, j, k);} } })
erf_qsatw
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_MicrophysicsUtils.H:171
erf_qsati
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsati(amrex::Real t, amrex::Real p, amrex::Real &qsati)
Definition: ERF_MicrophysicsUtils.H:160
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
SAM::NewtonIterSat
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real NewtonIterSat(int &i, int &j, int &k, const int &SAM_moisture_type, const amrex::Real &fac_cond, const amrex::Real &fac_fus, const amrex::Real &fac_sub, const amrex::Real &an, const amrex::Real &bn, const amrex::Array4< amrex::Real > &tabs_array, const amrex::Array4< amrex::Real > &pres_array, const amrex::Array4< amrex::Real > &qv_array, const amrex::Array4< amrex::Real > &qc_array, const amrex::Array4< amrex::Real > &qi_array, const amrex::Array4< amrex::Real > &qn_array, const amrex::Array4< amrex::Real > &qt_array)
Definition: ERF_SAM.H:160
SAM::mic_fab_vars
amrex::Array< FabPtr, MicVar::NumVars > mic_fab_vars
Definition: ERF_SAM.H:330
SAM::m_rdOcp
amrex::Real m_rdOcp
Definition: ERF_SAM.H:319
SAM::m_fac_fus
amrex::Real m_fac_fus
Definition: ERF_SAM.H:317
SAM::m_do_cond
bool m_do_cond
Definition: ERF_SAM.H:320
SAM::m_fac_cond
amrex::Real m_fac_cond
Definition: ERF_SAM.H:316
SAM::m_geom
amrex::Geometry m_geom
Definition: ERF_SAM.H:301
SAM::m_fac_sub
amrex::Real m_fac_sub
Definition: ERF_SAM.H:318
SAM::m_real_width
int m_real_width
Definition: ERF_SAM.H:304
MicVar::pres
@ pres
Definition: ERF_SAM.H:33
MicVar::qci
@ qci
Definition: ERF_SAM.H:39
MicVar::qv
@ qv
Definition: ERF_SAM.H:37
MicVar::qt
@ qt
Definition: ERF_SAM.H:35
MicVar::qn
@ qn
Definition: ERF_SAM.H:36
MicVar::theta
@ theta
Definition: ERF_SAM.H:31
MicVar::qcl
@ qcl
Definition: ERF_SAM.H:38
MicVar::tabs
@ tabs
Definition: ERF_SAM.H:32
SolverChoice::moisture_type
MoistureType moisture_type
Definition: ERF_DataStruct.H:1206

Referenced by Advance().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Compute_Coefficients()

void SAM::Compute_Coefficients ( )
132 {
133  auto accrrc_t = accrrc.table();
134  auto accrsi_t = accrsi.table();
135  auto accrsc_t = accrsc.table();
136  auto coefice_t = coefice.table();
137  auto evaps1_t = evaps1.table();
138  auto evaps2_t = evaps2.table();
139  auto accrgi_t = accrgi.table();
140  auto accrgc_t = accrgc.table();
141  auto evapg1_t = evapg1.table();
142  auto evapg2_t = evapg2.table();
143  auto evapr1_t = evapr1.table();
144  auto evapr2_t = evapr2.table();
145 
146  auto rho1d_t = rho1d.table();
147  auto pres1d_t = pres1d.table();
148  auto tabs1d_t = tabs1d.table();
149 
150  Real gam3 = erf_gammafff(three );
151  Real gamr1 = erf_gammafff(three+b_rain );
152  Real gamr2 = erf_gammafff((Real(5.0)+b_rain)/two);
153  Real gams1 = erf_gammafff(three+b_snow );
154  Real gams2 = erf_gammafff((Real(5.0)+b_snow)/two);
155  Real gamg1 = erf_gammafff(three+b_grau );
156  Real gamg2 = erf_gammafff((Real(5.0)+b_grau)/two);
157 
158  // calculate the plane average variables
159  PlaneAverage rho_ave(mic_fab_vars[MicVar::rho].get(), m_geom, m_axis);
160  PlaneAverage theta_ave(mic_fab_vars[MicVar::theta].get(), m_geom, m_axis);
161  PlaneAverage qv_ave(mic_fab_vars[MicVar::qv].get(), m_geom, m_axis);
162  rho_ave.compute_averages(ZDir(), rho_ave.field());
163  theta_ave.compute_averages(ZDir(), theta_ave.field());
164  qv_ave.compute_averages(ZDir(), qv_ave.field());
165 
166  // get host variable rho, and rhotheta
167  int ncell = rho_ave.ncell_line();
168 
169  Gpu::HostVector<Real> rho_h(ncell), theta_h(ncell), qv_h(ncell);
170  rho_ave.line_average(0, rho_h);
171  theta_ave.line_average(0, theta_h);
172  qv_ave.line_average(0, qv_h);
173 
174  // copy data to device
175  Gpu::DeviceVector<Real> rho_d(ncell), theta_d(ncell), qv_d(ncell);
176  Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin());
177  Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin());
178  Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin());
179  Gpu::streamSynchronize();
180 
181  Real* rho_dptr = rho_d.data();
182  Real* theta_dptr = theta_d.data();
183  Real* qv_dptr = qv_d.data();
184 
185  ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept
186  {
187  Real RhoTheta = rho_dptr[k]*theta_dptr[k];
188  Real pressure = getPgivenRTh(RhoTheta, qv_dptr[k]);
189  rho1d_t(k) = rho_dptr[k];
190  pres1d_t(k) = pressure*Real(0.01);
191  // NOTE: Limit the temperature to the melting point of ice to avoid a divide by
192  // 0 condition when computing the cold evaporation coefficients. This should
193  // not affect results since evaporation requires snow/graupel to be present
194  // and thus T<Real(273.16)
195  tabs1d_t(k) = std::min(getTgivenRandRTh(rho_dptr[k], RhoTheta, qv_dptr[k]),Real(273.16));
196  });
197 
198  if(round(gam3) != 2) {
199  std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl;
200  std::exit(-1);
201  }
202 
203  // Populate all the coefficients
204  ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept
205  {
206  Real Prefactor;
207  Real pratio = sqrt(Real(1.29) / rho1d_t(k));
208  //Real rrr1 = Real(393.0)/(tabs1d_t(k)+Real(120.0))*std::pow((tabs1d_t(k)/Real(273.0)),Real(1.5));
209  //Real rrr2 = std::pow((tabs1d_t(k)/Real(273.0)),Real(1.94))*(Real(1000.0)/pres1d_t(k));
210  Real estw = Real(100.0)*erf_esatw(tabs1d_t(k));
211  Real esti = Real(100.0)*erf_esati(tabs1d_t(k));
212 
213  // accretion by snow:
214  Real coef1 = fourth * PI * nzeros * a_snow * gams1 * pratio/pow((PI * rhos * nzeros/rho1d_t(k) ) , ((three+b_snow)/Real(4.0)));
215  Real coef2 = exp(Real(0.025)*(tabs1d_t(k) - Real(273.15)));
216  accrsi_t(k) = coef1 * coef2 * esicoef;
217  accrsc_t(k) = coef1 * esccoef;
218  coefice_t(k) = coef2;
219 
220  // evaporation of snow:
221  coef1 = (lsub/(tabs1d_t(k)*R_v)-one)*lsub/(therco*tabs1d_t(k));
222  coef2 = R_v * R_d / (diffelq * esti);
223  Prefactor = two * PI * nzeros / (rho1d_t(k) * (coef1 + coef2));
224  Prefactor *= (two/PI); // Shape factor snow
225  evaps1_t(k) = Prefactor * Real(0.65) * sqrt(rho1d_t(k) / (PI * rhos * nzeros));
226  evaps2_t(k) = Prefactor * Real(0.44) * sqrt(a_snow * rho1d_t(k) / muelq) * gams2
227  * sqrt(pratio) * pow(rho1d_t(k) / (PI * rhos * nzeros) , ((Real(5.0)+b_snow)/Real(8.0)));
228 
229  // accretion by graupel:
230  coef1 = fourth*PI*nzerog*a_grau*gamg1*pratio/pow((PI*rhog*nzerog/rho1d_t(k)) , ((three+b_grau)/Real(4.0)));
231  coef2 = exp(Real(0.025)*(tabs1d_t(k) - Real(273.15)));
232  accrgi_t(k) = coef1 * coef2 * egicoef;
233  accrgc_t(k) = coef1 * egccoef;
234 
235  // evaporation of graupel:
236  coef1 = (lsub/(tabs1d_t(k)*R_v)-one)*lsub/(therco*tabs1d_t(k));
237  coef2 = R_v * R_d / (diffelq * esti);
238  Prefactor = two * PI * nzerog / (rho1d_t(k) * (coef1 + coef2)); // Shape factor for graupel is 1
239  evapg1_t(k) = Prefactor * Real(0.78) * sqrt(rho1d_t(k) / (PI * rhog * nzerog));
240  evapg2_t(k) = Prefactor * Real(0.31) * sqrt(a_grau * rho1d_t(k) / muelq) * gamg2
241  * sqrt(pratio) * pow(rho1d_t(k) / (PI * rhog * nzerog) , ((Real(5.0)+b_grau)/Real(8.0)));
242 
243  // accretion by rain:
244  accrrc_t(k) = fourth * PI * nzeror * a_rain * gamr1 * pratio/pow((PI * rhor * nzeror / rho1d_t(k)) , ((three+b_rain)/Real(4.)))* erccoef;
245 
246  // evaporation of rain:
247  coef1 = (lcond/(tabs1d_t(k)*R_v)-one)*lcond/(therco*tabs1d_t(k));
248  coef2 = R_v * R_d / (diffelq * estw);
249  Prefactor = two * PI * nzeror / (rho1d_t(k) * (coef1 + coef2)); // Shape factor for rain is 1
250  evapr1_t(k) = Prefactor * Real(0.78) * sqrt(rho1d_t(k) / (PI * rhor * nzeror));
251  evapr2_t(k) = Prefactor * Real(0.31) * sqrt(a_rain * rho1d_t(k) / muelq) * gamr2
252  * sqrt(pratio) * pow(rho1d_t(k) / (PI * rhor * nzeror) , ((Real(5.0)+b_rain)/Real(8.0)));
253  });
254 }
rhog
constexpr amrex::Real rhog
Definition: ERF_Constants.H:40
R_v
constexpr amrex::Real R_v
Definition: ERF_Constants.H:21
three
constexpr amrex::Real three
Definition: ERF_Constants.H:9
muelq
constexpr amrex::Real muelq
Definition: ERF_Constants.H:85
nzerog
constexpr amrex::Real nzerog
Definition: ERF_Constants.H:69
two
constexpr amrex::Real two
Definition: ERF_Constants.H:8
a_grau
constexpr amrex::Real a_grau
Definition: ERF_Constants.H:52
lsub
constexpr amrex::Real lsub
Definition: ERF_Constants.H:78
esicoef
constexpr amrex::Real esicoef
Definition: ERF_Constants.H:63
diffelq
constexpr amrex::Real diffelq
Definition: ERF_Constants.H:83
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:12
therco
constexpr amrex::Real therco
Definition: ERF_Constants.H:84
b_grau
constexpr amrex::Real b_grau
Definition: ERF_Constants.H:53
egccoef
constexpr amrex::Real egccoef
Definition: ERF_Constants.H:64
egicoef
constexpr amrex::Real egicoef
Definition: ERF_Constants.H:65
b_rain
constexpr amrex::Real b_rain
Definition: ERF_Constants.H:49
lcond
constexpr amrex::Real lcond
Definition: ERF_Constants.H:76
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:16
nzeror
constexpr amrex::Real nzeror
Definition: ERF_Constants.H:67
rhos
constexpr amrex::Real rhos
Definition: ERF_Constants.H:39
esccoef
constexpr amrex::Real esccoef
Definition: ERF_Constants.H:62
rhor
constexpr amrex::Real rhor
Definition: ERF_Constants.H:38
a_rain
constexpr amrex::Real a_rain
Definition: ERF_Constants.H:48
nzeros
constexpr amrex::Real nzeros
Definition: ERF_Constants.H:68
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:20
a_snow
constexpr amrex::Real a_snow
Definition: ERF_Constants.H:50
b_snow
constexpr amrex::Real b_snow
Definition: ERF_Constants.H:51
erccoef
constexpr amrex::Real erccoef
Definition: ERF_Constants.H:61
ZDir
DirectionSelector< 2 > ZDir
Definition: ERF_DirectionSelector.H:38
getTgivenRandRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:46
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:81
get
pp get("wavelength", wavelength)
erf_esati
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_esati(amrex::Real t)
Definition: ERF_MicrophysicsUtils.H:26
erf_gammafff
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_gammafff(amrex::Real x)
Definition: ERF_MicrophysicsUtils.H:15
erf_esatw
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_esatw(amrex::Real t, bool use_empirical=false)
Definition: ERF_MicrophysicsUtils.H:68
rho_h
Vector< Real > rho_h(khi+1, zero)
rho_d
Gpu::DeviceVector< Real > rho_d(khi+1, zero)
PlaneAverage
Definition: ERF_PlaneAverage.H:14
SAM::rho1d
amrex::TableData< amrex::Real, 1 > rho1d
Definition: ERF_SAM.H:347
SAM::evapr2
amrex::TableData< amrex::Real, 1 > evapr2
Definition: ERF_SAM.H:344
SAM::evapg1
amrex::TableData< amrex::Real, 1 > evapg1
Definition: ERF_SAM.H:341
SAM::accrgi
amrex::TableData< amrex::Real, 1 > accrgi
Definition: ERF_SAM.H:339
SAM::evapg2
amrex::TableData< amrex::Real, 1 > evapg2
Definition: ERF_SAM.H:342
SAM::m_axis
int m_axis
Definition: ERF_SAM.H:313
SAM::evaps1
amrex::TableData< amrex::Real, 1 > evaps1
Definition: ERF_SAM.H:337
SAM::accrsi
amrex::TableData< amrex::Real, 1 > accrsi
Definition: ERF_SAM.H:334
SAM::accrrc
amrex::TableData< amrex::Real, 1 > accrrc
Definition: ERF_SAM.H:333
SAM::pres1d
amrex::TableData< amrex::Real, 1 > pres1d
Definition: ERF_SAM.H:348
SAM::nlev
int nlev
Definition: ERF_SAM.H:310
SAM::evapr1
amrex::TableData< amrex::Real, 1 > evapr1
Definition: ERF_SAM.H:343
SAM::evaps2
amrex::TableData< amrex::Real, 1 > evaps2
Definition: ERF_SAM.H:338
SAM::accrgc
amrex::TableData< amrex::Real, 1 > accrgc
Definition: ERF_SAM.H:340
SAM::accrsc
amrex::TableData< amrex::Real, 1 > accrsc
Definition: ERF_SAM.H:335
SAM::coefice
amrex::TableData< amrex::Real, 1 > coefice
Definition: ERF_SAM.H:336
SAM::tabs1d
amrex::TableData< amrex::Real, 1 > tabs1d
Definition: ERF_SAM.H:349
MicVar::rho
@ rho
Definition: ERF_SAM.H:30

Referenced by Update_Micro_Vars().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Copy_Micro_to_State()

void SAM::Copy_Micro_to_State ( amrex::MultiFab &  cons_in)
overridevirtual

Updates conserved and microphysics variables in the provided MultiFabs from the internal MultiFabs that store Microphysics module data.

Parameters
[out]consConserved variables
[out]qmoistqv, qc, qi, qr, qs, qg

Reimplemented from NullMoist.

16 {
17  // Get the temperature, density, theta, qt and qp from input
18  for ( MFIter mfi(cons,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
19  const auto& box3d = mfi.tilebox();
20 
21  auto states_arr = cons.array(mfi);
22 
23  auto rho_arr = mic_fab_vars[MicVar::rho]->array(mfi);
24  auto theta_arr = mic_fab_vars[MicVar::theta]->array(mfi);
25 
26  auto qv_arr = mic_fab_vars[MicVar::qv]->array(mfi);
27  auto qc_arr = mic_fab_vars[MicVar::qcl]->array(mfi);
28  auto qi_arr = mic_fab_vars[MicVar::qci]->array(mfi);
29 
30  auto qpr_arr = mic_fab_vars[MicVar::qpr]->array(mfi);
31  auto qps_arr = mic_fab_vars[MicVar::qps]->array(mfi);
32  auto qpg_arr = mic_fab_vars[MicVar::qpg]->array(mfi);
33 
34  // get potential total density, temperature, qt, qp
35  ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
36  {
37  states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k);
38 
39  states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*std::max(Real(0),qv_arr(i,j,k));
40  states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*std::max(Real(0),qc_arr(i,j,k));
41  states_arr(i,j,k,RhoQ3_comp) = rho_arr(i,j,k)*std::max(Real(0),qi_arr(i,j,k));
42 
43  states_arr(i,j,k,RhoQ4_comp) = rho_arr(i,j,k)*std::max(Real(0),qpr_arr(i,j,k));
44  states_arr(i,j,k,RhoQ5_comp) = rho_arr(i,j,k)*std::max(Real(0),qps_arr(i,j,k));
45  states_arr(i,j,k,RhoQ6_comp) = rho_arr(i,j,k)*std::max(Real(0),qpg_arr(i,j,k));
46  });
47  }
48 
49  // Fill interior ghost cells and periodic boundaries
50  cons.FillBoundary(m_geom.periodicity());
51 }
RhoQ4_comp
#define RhoQ4_comp
Definition: ERF_IndexDefines.H:45
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
RhoQ2_comp
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:43
RhoQ3_comp
#define RhoQ3_comp
Definition: ERF_IndexDefines.H:44
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
RhoQ6_comp
#define RhoQ6_comp
Definition: ERF_IndexDefines.H:47
RhoQ5_comp
#define RhoQ5_comp
Definition: ERF_IndexDefines.H:46
qv_arr
auto qv_arr
Definition: ERF_InitCustomPert_MultiSpeciesBubble.H:210
rho_arr
auto rho_arr
Definition: ERF_UpdateWSubsidence_SineMassFlux.H:3
MicVar::qpr
@ qpr
Definition: ERF_SAM.H:42
MicVar::qpg
@ qpg
Definition: ERF_SAM.H:44
MicVar::qps
@ qps
Definition: ERF_SAM.H:43
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:158

Referenced by Update_State_Vars().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Copy_State_to_Micro()

void SAM::Copy_State_to_Micro ( const amrex::MultiFab &  cons_in)
overridevirtual

Initializes the Microphysics module.

Parameters
[in]cons_inConserved variables input

Reimplemented from NullMoist.

80 {
81  // Get the temperature, density, theta, qt and qp from input
82  for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
83  const auto& box3d = mfi.growntilebox();
84 
85  auto states_array = cons_in.array(mfi);
86 
87  // Non-precipitating
88  auto qv_array = mic_fab_vars[MicVar::qv]->array(mfi);
89  auto qc_array = mic_fab_vars[MicVar::qcl]->array(mfi);
90  auto qi_array = mic_fab_vars[MicVar::qci]->array(mfi);
91  auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi);
92  auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi);
93 
94  // Precipitating
95  auto qpr_array = mic_fab_vars[MicVar::qpr]->array(mfi);
96  auto qps_array = mic_fab_vars[MicVar::qps]->array(mfi);
97  auto qpg_array = mic_fab_vars[MicVar::qpg]->array(mfi);
98  auto qp_array = mic_fab_vars[MicVar::qp]->array(mfi);
99 
100  auto rho_array = mic_fab_vars[MicVar::rho]->array(mfi);
101  auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi);
102  auto tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi);
103  auto pres_array = mic_fab_vars[MicVar::pres]->array(mfi);
104 
105  // Get pressure, theta, temperature, density, and qt, qp
106  ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
107  {
108  rho_array(i,j,k) = states_array(i,j,k,Rho_comp);
109  theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp);
110 
111  qv_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp));
112  qc_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp));
113  qi_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ3_comp)/states_array(i,j,k,Rho_comp));
114  qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k);
115  qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k);
116 
117  qpr_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ4_comp)/states_array(i,j,k,Rho_comp));
118  qps_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ5_comp)/states_array(i,j,k,Rho_comp));
119  qpg_array(i,j,k) = std::max(Real(0),states_array(i,j,k,RhoQ6_comp)/states_array(i,j,k,Rho_comp));
120  qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
121 
122  tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),
123  states_array(i,j,k,RhoTheta_comp),
124  qv_array(i,j,k));
125  pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)) * Real(0.01);
126  });
127  }
128 }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
MicVar::qp
@ qp
Definition: ERF_SAM.H:41

Referenced by Update_Micro_Vars().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Define()

void SAM::Define ( SolverChoice &  sc)
inlineoverridevirtual

Reimplemented from NullMoist.

79  {
80  m_fac_cond = lcond / sc.c_p;
81  m_fac_fus = lfus / sc.c_p;
82  m_fac_sub = lsub / sc.c_p;
83  m_axis = sc.ave_plane;
84  m_rdOcp = sc.rdOcp;
85  m_do_cond = (!sc.use_shoc);
86  }
lfus
constexpr amrex::Real lfus
Definition: ERF_Constants.H:77
SolverChoice::rdOcp
amrex::Real rdOcp
Definition: ERF_DataStruct.H:1156
SolverChoice::c_p
amrex::Real c_p
Definition: ERF_DataStruct.H:1155
SolverChoice::use_shoc
bool use_shoc
Definition: ERF_DataStruct.H:1187
SolverChoice::ave_plane
int ave_plane
Definition: ERF_DataStruct.H:1220

◆ IceFall()

void SAM::IceFall ( const SolverChoice &  sc)

Sedimentation of cloud ice (A32)

11  {
12 
13  if(sc.moisture_type == MoistureType::SAM_NoIce ||
14  sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce)
15  return;
16 
17  Real dtn = dt;
18  Real coef = dtn/m_dzmin;
19 
20  auto domain = m_geom.Domain();
21  int k_lo = domain.smallEnd(2);
22  int k_hi = domain.bigEnd(2);
23 
24  auto qcl = mic_fab_vars[MicVar::qcl];
25  auto qci = mic_fab_vars[MicVar::qci];
26  auto qn = mic_fab_vars[MicVar::qn];
27  auto qt = mic_fab_vars[MicVar::qt];
28  auto rho = mic_fab_vars[MicVar::rho];
29  auto tabs = mic_fab_vars[MicVar::tabs];
30 
31  MultiFab fz;
32  IntVect ng = qcl->nGrowVect();
33  BoxArray ba = qcl->boxArray();
34  DistributionMapping dm = qcl->DistributionMap();
35  fz.define(convert(ba, IntVect(0,0,1)), dm, 1, ng);
36  fz.setVal(0.);
37 
38  for (MFIter mfi(fz, TileNoZ()); mfi.isValid(); ++mfi) {
39  auto qci_array = qci->array(mfi);
40  auto rho_array = rho->array(mfi);
41  auto fz_array = fz.array(mfi);
42 
43  const auto& box3d = mfi.tilebox();
44 
45  ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
46  {
47  Real rho_avg, qci_avg;
48  if (k==k_lo) {
49  rho_avg = rho_array(i,j,k);
50  qci_avg = qci_array(i,j,k);
51  } else if (k==k_hi+1) {
52  rho_avg = rho_array(i,j,k-1);
53  qci_avg = qci_array(i,j,k-1);
54  } else {
55  rho_avg = myhalf*(rho_array(i,j,k-1) + rho_array(i,j,k));
56  qci_avg = myhalf*(qci_array(i,j,k-1) + qci_array(i,j,k));
57  }
58  Real vt_ice = min( Real(0.4),
59  Real(8.66)*std::pow( (amrex::max(Real(0),qci_avg)+Real(1.e-10)) , Real(0.24)) );
60 
61  // NOTE: Fz is the sedimentation flux from the advective operator.
62  // In the terrain-following coordinate system, the z-deriv in
63  // the divergence uses the normal velocity (Omega). However,
64  // there are no u/v components to the sedimentation velocity.
65  // Therefore, we simply end up with a division by detJ when
66  // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
67  fz_array(i,j,k) = rho_avg*vt_ice*qci_avg;
68  });
69  }
70 
71  // Compute number of substeps from maximum terminal velocity
72  Real wt_max;
73  int n_substep;
74  auto const& ma_fz_arr = fz.const_arrays();
75  GpuTuple<Real> max = ParReduce(TypeList<ReduceOpMax>{},
76  TypeList<Real>{},
77  fz, IntVect(0),
78  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
79  -> GpuTuple<Real>
80  {
81  return { ma_fz_arr[box_no](i,j,k) };
82  });
83  wt_max = get<0>(max) + std::numeric_limits<Real>::epsilon();
84  n_substep = int( std::ceil(wt_max * coef / CFL_MAX) );
85  AMREX_ALWAYS_ASSERT(n_substep >= 1);
86  coef /= Real(n_substep);
87  dtn /= Real(n_substep);
88 
89  // Substep the vertical advection
90  for (int nsub(0); nsub<n_substep; ++nsub) {
91  for (MFIter mfi(*qci, TileNoZ()); mfi.isValid(); ++mfi) {
92  auto qci_array = qci->array(mfi);
93  auto qn_array = qn->array(mfi);
94  auto qt_array = qt->array(mfi);
95  auto rho_array = rho->array(mfi);
96  auto fz_array = fz.array(mfi);
97 
98  const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4<const Real>{};
99 
100  const auto& tbx = mfi.tilebox();
101  const auto& tbz = mfi.tilebox(IntVect(0,0,1),IntVect(0));
102 
103  // Update vertical flux every substep
104  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
105  {
106  Real rho_avg, qci_avg;
107  if (k==k_lo) {
108  rho_avg = rho_array(i,j,k);
109  qci_avg = qci_array(i,j,k);
110  } else if (k==k_hi+1) {
111  rho_avg = rho_array(i,j,k-1);
112  qci_avg = qci_array(i,j,k-1);
113  } else {
114  rho_avg = myhalf*(rho_array(i,j,k-1) + rho_array(i,j,k));
115  qci_avg = myhalf*(qci_array(i,j,k-1) + qci_array(i,j,k));
116  }
117  Real vt_ice = min( Real(0.4),
118  Real(8.66)*std::pow( (amrex::max(Real(0),qci_avg)+Real(1.e-10)) , Real(0.24)) );
119 
120  // NOTE: Fz is the sedimentation flux from the advective operator.
121  // In the terrain-following coordinate system, the z-deriv in
122  // the divergence uses the normal velocity (Omega). However,
123  // there are no u/v components to the sedimentation velocity.
124  // Therefore, we simply end up with a division by detJ when
125  // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
126  fz_array(i,j,k) = rho_avg*vt_ice*qci_avg;
127  });
128 
129  // Update precip every substep
130  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
131  {
132  // Jacobian determinant
133  Real dJinv = (dJ_array) ? one/dJ_array(i,j,k) : one;
134 
135  //==================================================
136  // Cloud ice sedimentation (A32)
137  //==================================================
138  Real dqi = dJinv * (one/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef;
139  dqi = std::max(-qci_array(i,j,k), dqi);
140 
141  // Add this increment to both non-precipitating and total water.
142  qci_array(i,j,k) += dqi;
143  qn_array(i,j,k) += dqi;
144  qt_array(i,j,k) += dqi;
145 
146  // NOTE: Sedimentation does not affect the potential temperature,
147  // but it does affect the liquid/ice static energy.
148  // No source to Theta occurs here.
149  });
150  } // mfi
151  } // nsub
152 }
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
AMREX_ALWAYS_ASSERT
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
rho
rho
Definition: ERF_InitCustomPert_Bubble.H:106
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
SAM::m_detJ_cc
amrex::MultiFab * m_detJ_cc
Definition: ERF_SAM.H:327
SAM::CFL_MAX
static constexpr amrex::Real CFL_MAX
Definition: ERF_SAM.H:295
SAM::m_dzmin
amrex::Real m_dzmin
Definition: ERF_SAM.H:323
MicVar_Kess::qcl
@ qcl
Definition: ERF_Kessler.H:29
MicVar_Kess::tabs
@ tabs
Definition: ERF_Kessler.H:24
MicVar_Kess::qt
@ qt
Definition: ERF_Kessler.H:27
MicVar_Morr::ng
@ ng
Definition: ERF_Morrison.H:48
MicVar_Morr::qn
@ qn
Definition: ERF_Morrison.H:33
MicVar_Morr::qci
@ qci
Definition: ERF_Morrison.H:36
module_model_constants::epsilon
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12

Referenced by Advance().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Init()

void SAM::Init ( const amrex::MultiFab &  cons_in,
const amrex::BoxArray &  grids,
const amrex::Geometry &  geom,
const amrex::Real &  dt_advance,
std::unique_ptr< amrex::MultiFab > &  z_phys_nd,
std::unique_ptr< amrex::MultiFab > &  detJ_cc 
)
overridevirtual

Initializes the Microphysics module.

Parameters
[in]cons_inConserved variables input
[in]qc_inCloud variables input
[in,out]qv_inVapor variables input
[in]qi_inIce variables input
[in]gridsThe boxes on which we will evolve the solution
[in]geomGeometry associated with these MultiFabs and grids
[in]dt_advanceTimestep for the advance

Reimplemented from NullMoist.

28 {
29  dt = dt_advance;
30  m_geom = geom;
31 
32  m_z_phys_nd = z_phys_nd.get();
33  m_detJ_cc = detJ_cc.get();
34 
35  MicVarMap.resize(m_qmoist_size);
36  MicVarMap = {MicVar::rain_accum, MicVar::snow_accum, MicVar::graup_accum};
37 
38  // initialize microphysics variables
39  for (auto ivar = 0; ivar < MicVar::NumVars; ++ivar) {
40  mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
41  1, cons_in.nGrowVect());
42  mic_fab_vars[ivar]->setVal(0.);
43  }
44 
45  // NOTE: For multi-level not all ranks will own a box.
46  // Furthermore, the plane average allocates space
47  // for the entire domain. We make this consistent.
48  nlev = m_geom.Domain().length(2);
49  zlo = m_geom.Domain().smallEnd(2);
50  zhi = m_geom.Domain().bigEnd(2);
51 
52  // parameters
53  accrrc.resize({zlo}, {zhi});
54  accrsi.resize({zlo}, {zhi});
55  accrsc.resize({zlo}, {zhi});
56  coefice.resize({zlo}, {zhi});
57  evaps1.resize({zlo}, {zhi});
58  evaps2.resize({zlo}, {zhi});
59  accrgi.resize({zlo}, {zhi});
60  accrgc.resize({zlo}, {zhi});
61  evapg1.resize({zlo}, {zhi});
62  evapg2.resize({zlo}, {zhi});
63  evapr1.resize({zlo}, {zhi});
64  evapr2.resize({zlo}, {zhi});
65 
66  // data (input)
67  rho1d.resize({zlo}, {zhi});
68  pres1d.resize({zlo}, {zhi});
69  tabs1d.resize({zlo}, {zhi});
70 }
SAM::zlo
int zlo
Definition: ERF_SAM.H:310
SAM::m_qmoist_size
int m_qmoist_size
Definition: ERF_SAM.H:289
SAM::m_z_phys_nd
amrex::MultiFab * m_z_phys_nd
Definition: ERF_SAM.H:326
SAM::MicVarMap
amrex::Vector< int > MicVarMap
Definition: ERF_SAM.H:298
SAM::zhi
int zhi
Definition: ERF_SAM.H:310
MicVar::NumVars
@ NumVars
Definition: ERF_SAM.H:49
MicVar::snow_accum
@ snow_accum
Definition: ERF_SAM.H:47
MicVar::rain_accum
@ rain_accum
Definition: ERF_SAM.H:46
MicVar::graup_accum
@ graup_accum
Definition: ERF_SAM.H:48

◆ NewtonIterSat()

AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real SAM::NewtonIterSat ( int &  i,
int &  j,
int &  k,
const int &  SAM_moisture_type,
const amrex::Real &  fac_cond,
const amrex::Real &  fac_fus,
const amrex::Real &  fac_sub,
const amrex::Real &  an,
const amrex::Real &  bn,
const amrex::Array4< amrex::Real > &  tabs_array,
const amrex::Array4< amrex::Real > &  pres_array,
const amrex::Array4< amrex::Real > &  qv_array,
const amrex::Array4< amrex::Real > &  qc_array,
const amrex::Array4< amrex::Real > &  qi_array,
const amrex::Array4< amrex::Real > &  qn_array,
const amrex::Array4< amrex::Real > &  qt_array 
)
inlinestatic
174  {
175  // Solution tolerance
176  amrex::Real tol = amrex::Real(1.0e-4);
177 
178  // Saturation moisture fractions
179  amrex::Real omn, domn;
180  amrex::Real qsat, dqsat;
181  amrex::Real qsatw, dqsatw;
182  amrex::Real qsati, dqsati;
183 
184  // Newton iteration vars
185  int niter;
186  amrex::Real fff, dfff, dtabs;
187  amrex::Real lstar, dlstar;
188  amrex::Real lstarw, lstari;
189  amrex::Real delta_qv, delta_qc, delta_qi;
190 
191  // Initial guess for temperature & pressure
192  amrex::Real tabs = tabs_array(i,j,k);
193  amrex::Real pres = pres_array(i,j,k);
194 
195  niter = 0;
196  dtabs = 1;
197  //==================================================
198  // Newton iteration to qv=qsat (cloud phase only)
199  //==================================================
200  do {
201  // Latent heats and their derivatives wrt to T
202  lstarw = fac_cond;
203  lstari = fac_sub;
204  domn = zero;
205 
206  // Saturation moisture fractions
207  erf_qsatw(tabs, pres, qsatw);
208  erf_qsati(tabs, pres, qsati);
209  erf_dtqsatw(tabs, pres, dqsatw);
210  erf_dtqsati(tabs, pres, dqsati);
211 
212  if (SAM_moisture_type == 1) {
213  // Cloud ice not permitted (condensation)
214  if(tabs >= tbgmax) {
215  omn = one;
216  }
217  // Cloud water not permitted (sublimation)
218  else if(tabs <= tbgmin) {
219  omn = zero;
220  }
221  // Mixed cloud phase (condensation & fusion)
222  else {
223  lstari = fac_fus;
224  omn = an*tabs-bn;
225  domn = an;
226  }
227  } else if (SAM_moisture_type == 2) {
228  omn = one;
229  domn = zero;
230  }
231 
232  // Linear combination of each component
233  qsat = omn * qsatw + (one-omn) * qsati;
234  dqsat = omn * dqsatw + (one-omn) * dqsati
235  + domn * qsatw - domn * qsati;
236  lstar = omn * lstarw + (one-omn) * lstari;
237  dlstar = domn * lstarw - domn * lstari;
238 
239  // Function for root finding:
240  // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat)
241  fff = -tabs + tabs_array(i,j,k) + lstar*(qv_array(i,j,k) - qsat);
242 
243  // Derivative of function (T_new iterated on)
244  dfff = -one + dlstar*(qv_array(i,j,k) - qsat) - lstar*dqsat;
245 
246  // Update the temperature
247  dtabs = -fff/dfff;
248  tabs += dtabs;
249 
250  // Update iteration
251  niter = niter+1;
252  } while(std::abs(dtabs) > tol && niter < 20);
253 
254  // Update qsat from last iteration (dq = dq/dt * dt)
255  qsat += dqsat*dtabs;
256 
257  // Changes in each component
258  delta_qv = qv_array(i,j,k) - qsat;
259  delta_qc = std::max(-qc_array(i,j,k), delta_qv * omn);
260  delta_qi = std::max(-qi_array(i,j,k), delta_qv * (one-omn));
261 
262  // Partition the change in non-precipitating q
263  qv_array(i,j,k) = qsat;
264  qc_array(i,j,k) += delta_qc;
265  qi_array(i,j,k) += delta_qi;
266  qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k);
267  qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k);
268 
269  // Return to temperature
270  return tabs;
271  }
erf_dtqsatw
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_dtqsatw(amrex::Real t, amrex::Real p, amrex::Real &dtqsatw)
Definition: ERF_MicrophysicsUtils.H:186
erf_dtqsati
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_dtqsati(amrex::Real t, amrex::Real p, amrex::Real &dtqsati)
Definition: ERF_MicrophysicsUtils.H:178
Here is the call graph for this function:

◆ Precip()

void SAM::Precip ( const SolverChoice &  sc)

Autoconversion (A30), Accretion (A28), Evaporation (A24)

11 {
12 
13  if (sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) return;
14 
15  Real powr1 = (three + b_rain) / Real(4.0);
16  Real powr2 = (Real(5.0) + b_rain) / Real(8.0);
17  Real pows1 = (three + b_snow) / Real(4.0);
18  Real pows2 = (Real(5.0) + b_snow) / Real(8.0);
19  Real powg1 = (three + b_grau) / Real(4.0);
20  Real powg2 = (Real(5.0) + b_grau) / Real(8.0);
21 
22  auto accrrc_t = accrrc.table();
23  auto accrsc_t = accrsc.table();
24  auto accrsi_t = accrsi.table();
25  auto accrgc_t = accrgc.table();
26  auto accrgi_t = accrgi.table();
27  auto coefice_t = coefice.table();
28  auto evapr1_t = evapr1.table();
29  auto evapr2_t = evapr2.table();
30  auto evaps1_t = evaps1.table();
31  auto evaps2_t = evaps2.table();
32  auto evapg1_t = evapg1.table();
33  auto evapg2_t = evapg2.table();
34 
35  Real fac_cond = m_fac_cond;
36  Real fac_sub = m_fac_sub;
37  Real fac_fus = m_fac_fus;
38  Real rdOcp = m_rdOcp;
39 
40  Real eps = std::numeric_limits<Real>::epsilon();
41 
42  Real dtn = dt;
43 
44  int SAM_moisture_type = 1;
45  if (sc.moisture_type == MoistureType::SAM_NoIce) {
46  SAM_moisture_type = 2;
47  }
48 
49  auto domain = m_geom.Domain();
50  int i_lo = domain.smallEnd(0);
51  int i_hi = domain.bigEnd(0);
52  int j_lo = domain.smallEnd(1);
53  int j_hi = domain.bigEnd(1);
54 
55  // get the temperature, density, theta, qt and qp from input
56  for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]),TilingIfNotGPU()); mfi.isValid(); ++mfi) {
57  auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi);
58  auto tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi);
59  auto pres_array = mic_fab_vars[MicVar::pres]->array(mfi);
60 
61  // Non-precipitating
62  auto qv_array = mic_fab_vars[MicVar::qv]->array(mfi);
63  auto qcl_array = mic_fab_vars[MicVar::qcl]->array(mfi);
64  auto qci_array = mic_fab_vars[MicVar::qci]->array(mfi);
65  auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi);
66  auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi);
67 
68  // Precipitating
69  auto qpr_array = mic_fab_vars[MicVar::qpr]->array(mfi);
70  auto qps_array = mic_fab_vars[MicVar::qps]->array(mfi);
71  auto qpg_array = mic_fab_vars[MicVar::qpg]->array(mfi);
72  auto qp_array = mic_fab_vars[MicVar::qp]->array(mfi);
73 
74  auto tbx = mfi.tilebox();
75  if (tbx.smallEnd(0) == i_lo) { tbx.growLo(0,-m_real_width); }
76  if (tbx.bigEnd(0) == i_hi) { tbx.growHi(0,-m_real_width); }
77  if (tbx.smallEnd(1) == j_lo) { tbx.growLo(1,-m_real_width); }
78  if (tbx.bigEnd(1) == j_hi) { tbx.growHi(1,-m_real_width); }
79 
80  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
81  {
82  //------- Autoconversion/accretion
83  Real omn, omp, omg;
84  Real qsat, qsatw, qsati;
85 
86  Real qcc, qii, qpr, qps, qpg;
87  Real dprc, dpsc, dpgc;
88  Real dpsi, dpgi;
89 
90  Real dqc, dqca, dqi, dqia, dqp;
91  Real dqpr, dqps, dqpg;
92 
93  Real auto_r, autos;
94  Real accrcr, accrcs, accris, accrcg, accrig;
95 
96  // Work to be done for autoc/accr or evap
97  if (qn_array(i,j,k)+qp_array(i,j,k) > zero) {
98  if (SAM_moisture_type == 2) {
99  omn = Real(1);
100  omp = Real(1);
101  omg = Real(0);
102  } else {
103  omn = std::max(Real(0),std::min(Real(1),(tabs_array(i,j,k)-tbgmin)*a_bg));
104  omp = std::max(Real(0),std::min(Real(1),(tabs_array(i,j,k)-tprmin)*a_pr));
105  omg = std::max(Real(0),std::min(Real(1),(tabs_array(i,j,k)-tgrmin)*a_gr));
106  }
107 
108  qcc = qcl_array(i,j,k);
109  qii = qci_array(i,j,k);
110 
111  qpr = qpr_array(i,j,k);
112  qps = qps_array(i,j,k);
113  qpg = qpg_array(i,j,k);
114 
115  //==================================================
116  // Autoconversion (A30/A31) and accretion (A27)
117  //==================================================
118  if (qn_array(i,j,k) > zero) {
119  accrcr = zero;
120  accrcs = zero;
121  accris = zero;
122  accrcg = zero;
123  accrig = zero;
124 
125  if (qcc > qcw0) {
126  auto_r = alphaelq;
127  } else {
128  auto_r = zero;
129  }
130 
131  if (qii > qci0) {
132  autos = betaelq*coefice_t(k);
133  } else {
134  autos = zero;
135  }
136 
137  if (omp > Real(0.001)) {
138  accrcr = accrrc_t(k);
139  }
140 
141  if (omp < Real(0.999) && omg < Real(0.999)) {
142  accrcs = accrsc_t(k);
143  accris = accrsi_t(k);
144  }
145 
146  if (omp < Real(0.999) && omg > Real(0.001)) {
147  accrcg = accrgc_t(k);
148  accrig = accrgi_t(k);
149  }
150 
151  // Autoconversion & accretion (sink for cloud comps)
152  dqca = dtn * auto_r * (qcc-qcw0);
153  dprc = dtn * accrcr * qcc * std::pow(qpr, powr1);
154  dpsc = dtn * accrcs * qcc * std::pow(qps, pows1);
155  dpgc = dtn * accrcg * qcc * std::pow(qpg, powg1);
156 
157  dqia = dtn * autos * (qii-qci0);
158  dpsi = dtn * accris * qii * std::pow(qps, pows1);
159  dpgi = dtn * accrig * qii * std::pow(qpg, powg1);
160 
161  // Rescale sinks to avoid negative cloud fractions
162  dqc = dqca + dprc + dpsc + dpgc;
163  dqi = dqia + dpsi + dpgi;
164  Real scalec = std::min(qcl_array(i,j,k),dqc) / (dqc + eps);
165  Real scalei = std::min(qci_array(i,j,k),dqi) / (dqi + eps);
166  dqca *= scalec; dprc *= scalec; dpsc *= scalec; dpgc *= scalec;
167  dqia *= scalei; dpsi *= scalei; dpgi *= scalei;
168  dqc = dqca + dprc + dpsc + dpgc;
169  dqi = dqia + dpsi + dpgi;
170 
171  // NOTE: Autoconversion of cloud water and ice are sources
172  // to qp, while accretion is a source to an individual
173  // precipitating component (e.g., qpr/qps/qpg). So we
174  // only split autoconversion with omega. The omega
175  // splitting does imply a latent heat source.
176 
177  // Partition formed precip componentss
178  dqpr = (dqca + dqia) * omp + dprc;
179  dqps = (dqca + dqia) * (one - omp) * (one - omg) + dpsc + dpsi;
180  dqpg = (dqca + dqia) * (one - omp) * omg + dpgc + dpgi;
181 
182  // Update the primitive state variables
183  qcl_array(i,j,k) -= dqc;
184  qci_array(i,j,k) -= dqi;
185  qpr_array(i,j,k) += dqpr;
186  qps_array(i,j,k) += dqps;
187  qpg_array(i,j,k) += dqpg;
188 
189  // Update the primitive derived vars
190  qn_array(i,j,k) = qcl_array(i,j,k) + qci_array(i,j,k);
191  qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k);
192  qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
193 
194  // Update temperature
195  tabs_array(i,j,k) += fac_fus * ( dqca * (one - omp) - dqia * omp );
196 
197  // Update theta
198  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), Real(100.0)*pres_array(i,j,k), rdOcp);
199  }
200 
201  //==================================================
202  // Evaporation (A24)
203  //==================================================
204  erf_qsatw(tabs_array(i,j,k),pres_array(i,j,k),qsatw);
205  erf_qsati(tabs_array(i,j,k),pres_array(i,j,k),qsati);
206  qsat = qsatw * omn + qsati * (one-omn);
207  if((qp_array(i,j,k) > zero) && (qv_array(i,j,k) < qsat)) {
208 
209  dqpr = evapr1_t(k)*sqrt(qpr) + evapr2_t(k)*pow(qpr,powr2);
210  dqps = evaps1_t(k)*sqrt(qps) + evaps2_t(k)*pow(qps,pows2);
211  dqpg = evapg1_t(k)*sqrt(qpg) + evapg2_t(k)*pow(qpg,powg2);
212 
213  // NOTE: This is always a sink for precipitating comps
214  // since qv<qsat and thus (1 - qv/qsat)>zero If we are
215  // in a super-saturated state (qv>qsat) the Newton
216  // iterations in Cloud() will have handled condensation.
217  dqpr *= dtn * (one - qv_array(i,j,k)/qsat);
218  dqps *= dtn * (one - qv_array(i,j,k)/qsat);
219  dqpg *= dtn * (one - qv_array(i,j,k)/qsat);
220 
221  // Limit to avoid negative moisture fractions
222  dqpr = std::min(qpr_array(i,j,k),dqpr);
223  dqps = std::min(qps_array(i,j,k),dqps);
224  dqpg = std::min(qpg_array(i,j,k),dqpg);
225  dqp = dqpr + dqps + dqpg;
226 
227  // Update the primitive state variables
228  qv_array(i,j,k) += dqp;
229  qpr_array(i,j,k) -= dqpr;
230  qps_array(i,j,k) -= dqps;
231  qpg_array(i,j,k) -= dqpg;
232 
233  // Update the primitive derived vars
234  qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k);
235  qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
236 
237  // Update temperature
238  tabs_array(i,j,k) -= fac_cond * dqpr + fac_sub * (dqps + dqpg);
239 
240  // Update theta
241  theta_array(i,j,k) = getThgivenTandP(tabs_array(i,j,k), Real(100.0)*pres_array(i,j,k), rdOcp);
242  }
243  }
244  });
245  }
246 }
a_bg
constexpr amrex::Real a_bg
Definition: ERF_Constants.H:87
a_gr
constexpr amrex::Real a_gr
Definition: ERF_Constants.H:89
qci0
constexpr amrex::Real qci0
Definition: ERF_Constants.H:57
betaelq
constexpr amrex::Real betaelq
Definition: ERF_Constants.H:59
alphaelq
constexpr amrex::Real alphaelq
Definition: ERF_Constants.H:58
qcw0
constexpr amrex::Real qcw0
Definition: ERF_Constants.H:56
tprmin
constexpr amrex::Real tprmin
Definition: ERF_Constants.H:43
tgrmin
constexpr amrex::Real tgrmin
Definition: ERF_Constants.H:45
a_pr
constexpr amrex::Real a_pr
Definition: ERF_Constants.H:88
MicVar_Morr::qpg
@ qpg
Definition: ERF_Morrison.H:41
MicVar_Morr::qps
@ qps
Definition: ERF_Morrison.H:40
MicVar_Morr::qpr
@ qpr
Definition: ERF_Morrison.H:39

Referenced by Advance().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ PrecipFall()

void SAM::PrecipFall ( const SolverChoice &  sc)

Precipitation fluxes P_{r/s/g} (A19)

Code modified from SAMXX, the C++ version of the SAM code.

Parameters
[in]hydro_typeType selection for the precipitation advection hydrodynamics scheme (0-3)
17 {
18  if(sc.moisture_type == MoistureType::SAM_NoPrecip_NoIce) return;
19 
20  Real rho_0 = Real(1.29);
21 
22  Real gamr3 = erf_gammafff(Real(4.0)+b_rain);
23  Real gams3 = erf_gammafff(Real(4.0)+b_snow);
24  Real gamg3 = erf_gammafff(Real(4.0)+b_grau);
25 
26  Real vrain = (a_rain*gamr3/Real(6.0))*pow((PI*rhor*nzeror),-crain);
27  Real vsnow = (a_snow*gams3/Real(6.0))*pow((PI*rhos*nzeros),-csnow);
28  Real vgrau = (a_grau*gamg3/Real(6.0))*pow((PI*rhog*nzerog),-cgrau);
29 
30  Real dtn = dt;
31  Real coef = dtn/m_dzmin;
32 
33  auto domain = m_geom.Domain();
34  int k_lo = domain.smallEnd(2);
35  int k_hi = domain.bigEnd(2);
36 
37  auto qpr = mic_fab_vars[MicVar::qpr];
38  auto qps = mic_fab_vars[MicVar::qps];
39  auto qpg = mic_fab_vars[MicVar::qpg];
40  auto qp = mic_fab_vars[MicVar::qp];
41  auto rho = mic_fab_vars[MicVar::rho];
42  auto tabs = mic_fab_vars[MicVar::tabs];
43  auto theta = mic_fab_vars[MicVar::theta];
44  auto rain_accum = mic_fab_vars[MicVar::rain_accum];
45  auto snow_accum = mic_fab_vars[MicVar::snow_accum];
46  auto graup_accum = mic_fab_vars[MicVar::graup_accum];
47 
48  auto ba = tabs->boxArray();
49  auto dm = tabs->DistributionMap();
50  auto ngrow = tabs->nGrowVect();
51 
52  MultiFab fz;
53  fz.define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow);
54 
55  int SAM_moisture_type = 1;
56  if (sc.moisture_type == MoistureType::SAM_NoIce) {
57  SAM_moisture_type = 2;
58  }
59 
60  // Precompute the vertical fluxes for CFL constraint
61  for (MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
62  auto qp_array = qp->array(mfi);
63  auto rho_array = rho->array(mfi);
64  auto tabs_array = tabs->array(mfi);
65  auto fz_array = fz.array(mfi);
66 
67  const auto& box3d = mfi.tilebox();
68 
69  ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
70  {
71  Real rho_avg, tab_avg, qp_avg;
72  if (k==k_lo) {
73  rho_avg = rho_array(i,j,k);
74  tab_avg = tabs_array(i,j,k);
75  qp_avg = qp_array(i,j,k);
76  } else if (k==k_hi+1) {
77  rho_avg = rho_array(i,j,k-1);
78  tab_avg = tabs_array(i,j,k-1);
79  qp_avg = qp_array(i,j,k-1);
80  } else {
81  rho_avg = myhalf*( rho_array(i,j,k-1) + rho_array(i,j,k));
82  tab_avg = myhalf*(tabs_array(i,j,k-1) + tabs_array(i,j,k));
83  qp_avg = myhalf*( qp_array(i,j,k-1) + qp_array(i,j,k));
84  }
85 
86  Real Pprecip = zero;
87  if(qp_avg > qp_threshold) {
88  Real omp, omg;
89  if (SAM_moisture_type == 2) {
90  omp = one;
91  omg = zero;
92  } else {
93  omp = std::max(Real(0),std::min(Real(1),(tab_avg-tprmin)*a_pr));
94  omg = std::max(Real(0),std::min(Real(1),(tab_avg-tgrmin)*a_gr));
95  }
96  Real qrr = omp*qp_avg;
97  Real qss = (one-omp)*(one-omg)*qp_avg;
98  Real qgg = (one-omp)*(omg)*qp_avg;
99  Pprecip = omp*vrain*std::pow(rho_avg*qrr,one+crain)
100  + (one-omp)*( (one-omg)*vsnow*std::pow(rho_avg*qss,one+csnow)
101  + omg *vgrau*std::pow(rho_avg*qgg,one+cgrau) );
102  }
103 
104  // NOTE: Fz is the sedimentation flux from the advective operator.
105  // In the terrain-following coordinate system, the z-deriv in
106  // the divergence uses the normal velocity (Omega). However,
107  // there are no u/v components to the sedimentation velocity.
108  // Therefore, we simply end up with a division by detJ when
109  // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
110  fz_array(i,j,k) = Pprecip * std::sqrt(rho_0/rho_avg);
111  });
112  }
113 
114  // Compute number of substeps from maximum terminal velocity
115  Real wt_max;
116  int n_substep;
117  auto const& ma_fz_arr = fz.const_arrays();
118  GpuTuple<Real> max = ParReduce(TypeList<ReduceOpMax>{},
119  TypeList<Real>{},
120  fz, IntVect(0),
121  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
122  -> GpuTuple<Real>
123  {
124  return { ma_fz_arr[box_no](i,j,k) };
125  });
126  wt_max = get<0>(max) + std::numeric_limits<Real>::epsilon();
127  n_substep = int( std::ceil(wt_max * coef / CFL_MAX) );
128  AMREX_ALWAYS_ASSERT(n_substep >= 1);
129  coef /= Real(n_substep);
130  dtn /= Real(n_substep);
131 
132  // Substep the vertical advection
133  for (int nsub(0); nsub<n_substep; ++nsub) {
134  for (MFIter mfi(*qp, TileNoZ()); mfi.isValid(); ++mfi) {
135  auto qpr_array = qpr->array(mfi);
136  auto qps_array = qps->array(mfi);
137  auto qpg_array = qpg->array(mfi);
138  auto qp_array = qp->array(mfi);
139  auto rho_array = rho->array(mfi);
140  auto tabs_array = tabs->array(mfi);
141  auto fz_array = fz.array(mfi);
142 
143  auto rain_accum_array = rain_accum->array(mfi);
144  auto snow_accum_array = snow_accum->array(mfi);
145  auto graup_accum_array = graup_accum->array(mfi);
146 
147  const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4<const Real>{};
148 
149  const auto& tbx = mfi.tilebox();
150  const auto& tbz = mfi.tilebox(IntVect(0,0,1),IntVect(0));
151 
152  // Update vertical flux every substep
153  ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
154  {
155  Real rho_avg, tab_avg, qp_avg;
156  if (k==k_lo) {
157  rho_avg = rho_array(i,j,k);
158  tab_avg = tabs_array(i,j,k);
159  qp_avg = qp_array(i,j,k);
160  } else if (k==k_hi+1) {
161  rho_avg = rho_array(i,j,k-1);
162  tab_avg = tabs_array(i,j,k-1);
163  qp_avg = qp_array(i,j,k-1);
164  } else {
165  rho_avg = myhalf*( rho_array(i,j,k-1) + rho_array(i,j,k));
166  tab_avg = myhalf*(tabs_array(i,j,k-1) + tabs_array(i,j,k));
167  qp_avg = myhalf*( qp_array(i,j,k-1) + qp_array(i,j,k));
168  }
169 
170  Real Pprecip = zero;
171  if(qp_avg > qp_threshold) {
172  Real omp, omg;
173  if (SAM_moisture_type == 2) {
174  omp = one;
175  omg = zero;
176  } else {
177  omp = std::max(Real(0),std::min(Real(1),(tab_avg-tprmin)*a_pr));
178  omg = std::max(Real(0),std::min(Real(1),(tab_avg-tgrmin)*a_gr));
179  }
180  Real qrr = omp*qp_avg;
181  Real qss = (one-omp)*(one-omg)*qp_avg;
182  Real qgg = (one-omp)*(omg)*qp_avg;
183  Pprecip = omp*vrain*std::pow(rho_avg*qrr,one+crain)
184  + (one-omp)*( (one-omg)*vsnow*std::pow(rho_avg*qss,one+csnow)
185  + omg *vgrau*std::pow(rho_avg*qgg,one+cgrau) );
186  }
187 
188  // NOTE: Fz is the sedimentation flux from the advective operator.
189  // In the terrain-following coordinate system, the z-deriv in
190  // the divergence uses the normal velocity (Omega). However,
191  // there are no u/v components to the sedimentation velocity.
192  // Therefore, we simply end up with a division by detJ when
193  // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
194  fz_array(i,j,k) = Pprecip * std::sqrt(rho_0/rho_avg);
195 
196  if(k==k_lo){
197  Real omp, omg;
198  if (SAM_moisture_type == 2) {
199  omp = one;
200  omg = zero;
201  } else {
202  omp = std::max(Real(0),std::min(Real(1),(tab_avg-tprmin)*a_pr));
203  omg = std::max(Real(0),std::min(Real(1),(tab_avg-tgrmin)*a_gr));
204  }
205  rain_accum_array(i,j,k) = rain_accum_array(i,j,k) + rho_avg*(omp*qp_avg)*vrain*dtn/rhor*Real(1000.0); // Divide by rho_water and convert to mm
206  snow_accum_array(i,j,k) = snow_accum_array(i,j,k) + rho_avg*(one-omp)*(one-omg)*qp_avg*vsnow*dtn/rhos*Real(1000.0); // Divide by rho_snow and convert to mm
207  graup_accum_array(i,j,k) = graup_accum_array(i,j,k) + rho_avg*(one-omp)*(omg)*qp_avg*vgrau*dtn/rhog*Real(1000.0); // Divide by rho_graupel and convert to mm
208  }
209  });
210 
211  // Update precip every substep
212  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
213  {
214  // Jacobian determinant
215  Real dJinv = (dJ_array) ? one/dJ_array(i,j,k) : one;
216 
217  //==================================================
218  // Precipitating sedimentation (A19)
219  //==================================================
220  Real dqp = dJinv * (one/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef;
221  Real omp, omg;
222  if (SAM_moisture_type == 2) {
223  omp = one;
224  omg = zero;
225  } else {
226  omp = std::max(Real(0),std::min(Real(1),(tabs_array(i,j,k)-tprmin)*a_pr));
227  omg = std::max(Real(0),std::min(Real(1),(tabs_array(i,j,k)-tgrmin)*a_gr));
228  }
229 
230  qpr_array(i,j,k) = std::max(Real(0), qpr_array(i,j,k) + dqp*omp);
231  qps_array(i,j,k) = std::max(Real(0), qps_array(i,j,k) + dqp*(one-omp)*(one-omg));
232  qpg_array(i,j,k) = std::max(Real(0), qpg_array(i,j,k) + dqp*(one-omp)*omg);
233  qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
234 
235  // NOTE: Sedimentation does not affect the potential temperature,
236  // but it does affect the liquid/ice static energy.
237  // No source to Theta occurs here.
238  });
239  } // mfi
240  } // nsub
241 }
csnow
constexpr amrex::Real csnow
Definition: ERF_Constants.H:92
cgrau
constexpr amrex::Real cgrau
Definition: ERF_Constants.H:93
qp_threshold
constexpr amrex::Real qp_threshold
Definition: ERF_Constants.H:70
crain
constexpr amrex::Real crain
Definition: ERF_Constants.H:91
rho_0
Real rho_0
Definition: ERF_InitCustomPert_ABL.H:4
LsmVar_MM5::theta
@ theta
Definition: ERF_MM5.H:20
MicVar_Kess::qp
@ qp
Definition: ERF_Kessler.H:31
MicVar_Kess::rain_accum
@ rain_accum
Definition: ERF_Kessler.H:33
MicVar_Morr::graup_accum
@ graup_accum
Definition: ERF_Morrison.H:52
MicVar_Morr::snow_accum
@ snow_accum
Definition: ERF_Morrison.H:51

Referenced by Advance().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Qmoist_Ptr()

amrex::MultiFab* SAM::Qmoist_Ptr ( const int &  varIdx)
inlineoverridevirtual

Reimplemented from NullMoist.

140  {
141  AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
142  return mic_fab_vars[MicVarMap[varIdx]].get();
143  }
Here is the call graph for this function:

◆ Qmoist_Restart_Vars()

void SAM::Qmoist_Restart_Vars ( const SolverChoice &  ,
std::vector< int > &  a_idx,
std::vector< std::string > &  a_names 
) const
inlineoverridevirtual

Reimplemented from NullMoist.

277  {
278  a_idx.clear();
279  a_names.clear();
280 
281  // NOTE: These are the indices to access into qmoist (not mic_fab_vars)
282  a_idx.push_back(0); a_names.push_back("RainAccum");
283  a_idx.push_back(1); a_names.push_back("SnowAccum");
284  a_idx.push_back(2); a_names.push_back("GraupAccum");
285  }

◆ Qmoist_Size()

int SAM::Qmoist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

149 { return SAM::m_qmoist_size; }

◆ Qstate_Moist_Size()

int SAM::Qstate_Moist_Size ( )
inlineoverridevirtual

Reimplemented from NullMoist.

152 { return SAM::n_qstate_moist_size; }
SAM::n_qstate_moist_size
int n_qstate_moist_size
Definition: ERF_SAM.H:292

◆ Set_dzmin()

void SAM::Set_dzmin ( const amrex::Real  dz_min)
inlineoverridevirtual

Reimplemented from NullMoist.

100  {
101  m_dzmin = dz_min;
102  }

◆ Set_RealWidth()

void SAM::Set_RealWidth ( const int  real_width)
inlineoverridevirtual

Reimplemented from NullMoist.

155 { m_real_width = real_width; }

◆ Update_Micro_Vars()

void SAM::Update_Micro_Vars ( amrex::MultiFab &  cons_in)
inlineoverridevirtual

Reimplemented from NullMoist.

114  {
115  this->Copy_State_to_Micro(cons_in);
116  this->Compute_Coefficients();
117  }
SAM::Copy_State_to_Micro
void Copy_State_to_Micro(const amrex::MultiFab &cons_in) override
Definition: ERF_InitSAM.cpp:79
SAM::Compute_Coefficients
void Compute_Coefficients()
Definition: ERF_InitSAM.cpp:131
Here is the call graph for this function:

◆ Update_State_Vars()

void SAM::Update_State_Vars ( amrex::MultiFab &  cons_in)
inlineoverridevirtual

Reimplemented from NullMoist.

121  {
122  this->Copy_Micro_to_State(cons_in);
123  }
SAM::Copy_Micro_to_State
void Copy_Micro_to_State(amrex::MultiFab &cons_in) override
Definition: ERF_UpdateSAM.cpp:15
Here is the call graph for this function:

Member Data Documentation

◆ accrgc

amrex::TableData<amrex::Real, 1> SAM::accrgc
private

◆ accrgi

amrex::TableData<amrex::Real, 1> SAM::accrgi
private

◆ accrrc

amrex::TableData<amrex::Real, 1> SAM::accrrc
private

◆ accrsc

amrex::TableData<amrex::Real, 1> SAM::accrsc
private

◆ accrsi

amrex::TableData<amrex::Real, 1> SAM::accrsi
private

◆ CFL_MAX

constexpr amrex::Real SAM::CFL_MAX = myhalf
staticconstexprprivate

◆ coefice

amrex::TableData<amrex::Real, 1> SAM::coefice
private

◆ dt

amrex::Real SAM::dt
private

Referenced by Advance().

◆ evapg1

amrex::TableData<amrex::Real, 1> SAM::evapg1
private

◆ evapg2

amrex::TableData<amrex::Real, 1> SAM::evapg2
private

◆ evapr1

amrex::TableData<amrex::Real, 1> SAM::evapr1
private

◆ evapr2

amrex::TableData<amrex::Real, 1> SAM::evapr2
private

◆ evaps1

amrex::TableData<amrex::Real, 1> SAM::evaps1
private

◆ evaps2

amrex::TableData<amrex::Real, 1> SAM::evaps2
private

◆ m_axis

int SAM::m_axis
private

Referenced by Define().

◆ m_detJ_cc

amrex::MultiFab* SAM::m_detJ_cc
private

◆ m_do_cond

bool SAM::m_do_cond
private

Referenced by Define().

◆ m_dzmin

amrex::Real SAM::m_dzmin
private

Referenced by Set_dzmin().

◆ m_fac_cond

amrex::Real SAM::m_fac_cond
private

Referenced by Define().

◆ m_fac_fus

amrex::Real SAM::m_fac_fus
private

Referenced by Define().

◆ m_fac_sub

amrex::Real SAM::m_fac_sub
private

Referenced by Define().

◆ m_geom

amrex::Geometry SAM::m_geom
private

◆ m_qmoist_size

int SAM::m_qmoist_size = 3
private

Referenced by Qmoist_Ptr(), and Qmoist_Size().

◆ m_rdOcp

amrex::Real SAM::m_rdOcp
private

Referenced by Define().

◆ m_real_width

int SAM::m_real_width {0}
private

Referenced by Set_RealWidth().

◆ m_z_phys_nd

amrex::MultiFab* SAM::m_z_phys_nd
private

◆ mic_fab_vars

amrex::Array<FabPtr, MicVar::NumVars> SAM::mic_fab_vars
private

Referenced by Qmoist_Ptr().

◆ MicVarMap

amrex::Vector<int> SAM::MicVarMap
private

Referenced by Qmoist_Ptr().

◆ n_qstate_moist_size

int SAM::n_qstate_moist_size = 6
private

Referenced by Qstate_Moist_Size().

◆ nlev

int SAM::nlev
private

◆ pres1d

amrex::TableData<amrex::Real, 1> SAM::pres1d
private

◆ qn1d

amrex::TableData<amrex::Real, 1> SAM::qn1d
private

◆ qt1d

amrex::TableData<amrex::Real, 1> SAM::qt1d
private

◆ qv1d

amrex::TableData<amrex::Real, 1> SAM::qv1d
private

◆ rho1d

amrex::TableData<amrex::Real, 1> SAM::rho1d
private

◆ tabs1d

amrex::TableData<amrex::Real, 1> SAM::tabs1d
private

◆ zhi

int SAM::zhi
private

◆ zlo

int SAM::zlo
private
The documentation for this class was generated from the following files: