15 auto dz =
m_geom.CellSize(2);
16 auto domain =
m_geom.Domain();
17 int k_lo = domain.smallEnd(2);
18 int k_hi = domain.bigEnd(2);
21 auto ba =
tabs->boxArray();
22 auto dm =
tabs->DistributionMap();
23 fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0);
29 for ( MFIter mfi(*
tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
39 const auto& tbx = mfi.tilebox();
44 ParallelFor(tbx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
46 qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
47 qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
48 qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
51 Real qcc, auto_r, accrr;
53 Real dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater;
55 Real pressure = pres_array(i,j,k);
56 erf_qsatw(tabs_array(i,j,k), pressure, qsat);
60 amrex::Warning(
"qsat computed as non-positive; setting to 0.!");
68 dq_clwater_to_rain = 0.0;
69 dq_rain_to_vapor = 0.0;
70 dq_vapor_to_clwater = 0.0;
71 dq_clwater_to_vapor = 0.0;
75 Real fac = 1.0 + (
L_v/
Cp_d)*dtqsat;
78 if (qv_array(i,j,k) > qsat) {
79 dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac));
84 if (qv_array(i,j,k) < qsat && qc_array(i,j,k) > 0.0) {
85 dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
88 if (qp_array(i,j,k) > 0.0 && qv_array(i,j,k) < qsat) {
89 Real C = 1.6 + 124.9*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.2046);
90 dq_rain_to_vapor = 1.0/(0.001*rho_array(i,j,k))*(1.0 - qv_array(i,j,k)/qsat)*C*std::pow(0.001*rho_array(i,j,k)*qp_array(i,j,k),0.525)/
91 (5.4e5 + 2.55e6/(pressure*qsat))*dtn;
94 dq_rain_to_vapor = std::min({qp_array(i,j,k), dq_rain_to_vapor});
100 if (qc_array(i,j,k) > 0.0) {
101 qcc = qc_array(i,j,k);
109 accrr = 2.2 * std::pow(qp_array(i,j,k) , 0.875);
110 dq_clwater_to_rain = dtn *(accrr*qcc + auto_r*(qcc -
qcw0));
113 dq_clwater_to_rain = std::min(dq_clwater_to_rain, qc_array(i,j,k));
116 qv_array(i,j,k) += -dq_vapor_to_clwater + dq_clwater_to_vapor + dq_rain_to_vapor;
117 qc_array(i,j,k) += dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain;
118 qp_array(i,j,k) += dq_clwater_to_rain - dq_rain_to_vapor;
120 Real theta_over_T = theta_array(i,j,k)/tabs_array(i,j,k);
121 theta_array(i,j,k) += theta_over_T * d_fac_cond * (dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor);
123 qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
124 qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
125 qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
127 qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k);
132 for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){
136 auto fz_array = fz.array(mfi);
137 const Box& tbz = mfi.tilebox();
139 ParallelFor(tbz, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
141 Real rho_avg, qp_avg;
144 rho_avg = rho_array(i,j,k);
145 qp_avg = qp_array(i,j,k);
146 }
else if (k==k_hi+1) {
147 rho_avg = rho_array(i,j,k-1);
148 qp_avg = qp_array(i,j,k-1);
150 rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k));
151 qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k));
154 qp_avg = std::max(0.0, qp_avg);
156 Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5);
164 fz_array(i,j,k) = rho_avg*V_terminal*qp_avg;
171 auto const& ma_fz_arr = fz.const_arrays();
172 GpuTuple<Real> max = ParReduce(TypeList<ReduceOpMax>{},
175 [=] AMREX_GPU_DEVICE (
int box_no,
int i,
int j,
int k) noexcept
178 return { ma_fz_arr[box_no](i,j,k) };
180 wt_max = get<0>(max) + std::numeric_limits<Real>::epsilon();
181 n_substep = int( std::ceil(wt_max * coef /
CFL_MAX) );
182 AMREX_ALWAYS_ASSERT(n_substep >= 1);
183 coef /= Real(n_substep);
184 dtn /= Real(n_substep);
187 for (
int nsub(0); nsub<n_substep; ++nsub) {
188 for ( MFIter mfi(*
tabs, TilingIfNotGPU()); mfi.isValid(); ++mfi ){
192 auto fz_array = fz.array(mfi);
194 const auto dJ_array = (
m_detJ_cc) ?
m_detJ_cc->const_array(mfi) : Array4<const Real>{};
196 const Box& tbx = mfi.tilebox();
197 const Box& tbz = mfi.tilebox(IntVect(0,0,1),IntVect(0));
200 ParallelFor(tbz, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
202 Real rho_avg, qp_avg;
204 rho_avg = rho_array(i,j,k);
205 qp_avg = qp_array(i,j,k);
206 }
else if (k==k_hi+1) {
207 rho_avg = rho_array(i,j,k-1);
208 qp_avg = qp_array(i,j,k-1);
210 rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k));
211 qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k));
214 qp_avg = std::max(0.0, qp_avg);
216 Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5);
224 fz_array(i,j,k) = rho_avg*V_terminal*qp_avg;
227 rain_accum_array(i,j,k) = rain_accum_array(i,j,k) + rho_avg*qp_avg*V_terminal*dtn/1000.0*1000.0;
232 ParallelFor(tbx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
235 Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0;
237 if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0;
238 if(std::fabs(fz_array(i,j,k )) < 1e-14) fz_array(i,j,k ) = 0.0;
239 Real dq_sed = dJinv * (1.0/rho_array(i,j,k)) * (fz_array(i,j,k+1) - fz_array(i,j,k)) * coef;
240 if(std::fabs(dq_sed) < 1e-14) dq_sed = 0.0;
242 qp_array(i,j,k) += dq_sed;
243 qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
250 if (solverChoice.
moisture_type == MoistureType::Kessler_NoRain){
253 for ( MFIter mfi(*
tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
261 const auto& box3d = mfi.tilebox();
266 ParallelFor(box3d, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept
268 qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
272 Real dq_clwater_to_vapor, dq_vapor_to_clwater;
274 Real pressure = pres_array(i,j,k);
275 erf_qsatw(tabs_array(i,j,k), pressure, qsat);
282 dq_vapor_to_clwater = 0.0;
283 dq_clwater_to_vapor = 0.0;
287 Real fac = 1.0 + (
L_v/
Cp_d)*dtqsat;
290 if (qv_array(i,j,k) > qsat){
291 dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac));
295 if (qv_array(i,j,k) < qsat && qc_array(i,j,k) > 0.0){
296 dq_clwater_to_vapor = std::min(qc_array(i,j,k), (qsat - qv_array(i,j,k))/(1.0 + fac));
299 qv_array(i,j,k) += -dq_vapor_to_clwater + dq_clwater_to_vapor;
300 qc_array(i,j,k) += dq_vapor_to_clwater - dq_clwater_to_vapor;
302 Real theta_over_T = theta_array(i,j,k)/tabs_array(i,j,k);
304 theta_array(i,j,k) += theta_over_T * d_fac_cond * (dq_vapor_to_clwater - dq_clwater_to_vapor);
306 qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
307 qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
309 qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k);
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:12
constexpr amrex::Real alphaelq
Definition: ERF_Constants.H:48
constexpr amrex::Real qcw0
Definition: ERF_Constants.H:46
constexpr amrex::Real L_v
Definition: ERF_Constants.H:16
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_dtqsatw(amrex::Real t, amrex::Real p, amrex::Real &dtqsatw)
Definition: ERF_MicrophysicsUtils.H:170
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw(amrex::Real t, amrex::Real p, amrex::Real &qsatw)
Definition: ERF_MicrophysicsUtils.H:158
amrex::Geometry m_geom
Definition: ERF_Kessler.H:153
amrex::MultiFab * m_detJ_cc
Definition: ERF_Kessler.H:178
amrex::Real m_fac_cond
Definition: ERF_Kessler.H:171
amrex::Array< FabPtr, MicVar_Kess::NumVars > mic_fab_vars
Definition: ERF_Kessler.H:181
static constexpr amrex::Real CFL_MAX
Definition: ERF_Kessler.H:147
@ qp
Definition: ERF_Kessler.H:39
@ qcl
Definition: ERF_Kessler.H:37
@ tabs
Definition: ERF_Kessler.H:32
@ pres
Definition: ERF_Kessler.H:33
@ rho
Definition: ERF_Kessler.H:30
@ theta
Definition: ERF_Kessler.H:31
@ qt
Definition: ERF_Kessler.H:35
@ rain_accum
Definition: ERF_Kessler.H:41
@ qv
Definition: ERF_Kessler.H:36
MoistureType moisture_type
Definition: ERF_DataStruct.H:664