3 #ifndef ERF_MODAL_AERO_WATERUPTAKE_H_
4 #define ERF_MODAL_AERO_WATERUPTAKE_H_
9 #include <AMReX_GpuComplex.H>
15 using yakl::fortran::parallel_for;
16 using yakl::fortran::SimpleBounds;
20 constexpr
static real
third = 1./3.;
21 constexpr
static real
pi43 =
PI*4.0/3.0;
22 constexpr
static real
huge_real = std::numeric_limits<real>::max();
23 constexpr
static int imax = 200;
27 MamConstituents& consti,
const real2d& h2ommr,
const real2d& t,
const real2d& pmid,
28 const real3d& dgnumdry_m,
const real3d& dgnumwet_m,
const real3d& qaerwat_m,
29 const real3d& wetdens_m,
const real2d& clear_rh_in)
31 real2d raer(
"raer",ncol,nlev);
32 real2d cldn(
"cldn",ncol,nlev);
33 real3d dgncur_a(
"dgncur",ncol,nlev,nmodes);
34 real3d dgncur_awet(
"dgncur_awet",ncol,nlev,nmodes);
35 real3d wetdens(
"wetdens",ncol,nlev,nmodes);
36 real3d qaerwat(
"qaerwat",ncol,nlev,nmodes);
38 real2d dryvolmr(
"dryvolmr",ncol,nlev);
40 real spechygro, spechygro_1;
46 real2d rh(
"rh",ncol,nlev);
50 real2d aerosol_water(
"aerosol_water",ncol,nlev);
55 real3d naer(
"naer",ncol,nlev,nmodes);
56 real3d dryvol(
"dryvol",ncol,nlev,nmodes);
57 real3d drymass(
"drymass",ncol,nlev,nmodes);
58 real3d dryrad(
"dryrad",ncol,nlev,nmodes);
59 real3d wetrad(
"wetrad",ncol,nlev,nmodes);
60 real3d wetvol(
"wetvol",ncol,nlev,nmodes);
61 real3d wtrvol(
"wtrvol",ncol,nlev,nmodes);
62 real1d rhcrystal(
"rhcrystal",nmodes);
63 real1d rhdeliques(
"rhdeliques",nmodes);
64 real1d specdens_1(
"specdens_1",nmodes);
65 real3d maer(
"maer",ncol,nlev,nmodes);
66 real3d hygro(
"hygro",ncol,nlev,nmodes);
80 yakl::memset(maer, 0.0);
81 yakl::memset(hygro, 0.);
84 compute_wetdens =
true;
86 dgncur_a = dgnumdry_m;
87 dgncur_awet = dgnumwet_m;
92 for (
auto m = 1; m <= nmodes; ++m) {
93 yakl::memset(dryvolmr, 0.);
96 consti.
get_mode_props(list_idx, m-1, sigmag, rhcrystal(m), rhdeliques(m));
101 for(
auto l = 0; l < nspec; ++l) {
109 specdens_1(m) = specdens;
110 spechygro_1 = spechygro;
114 for(
auto k = top_lev; k <= nlev; ++k) {
115 for(
auto i = 1; i <= ncol; ++i) {
117 maer(i,k,m) = maer(i,k,m) + duma;
118 dumb = duma/specdens;
119 dryvolmr(i,k) = dryvolmr(i,k) + dumb;
120 hygro(i,k,m) = hygro(i,k,m) + dumb*spechygro;
127 for(
auto k = top_lev; k <= nlev; ++k) {
128 for(
auto i = 1; i <= ncol; ++i) {
130 if (dryvolmr(i,k) > 1.0e-30) {
131 hygro(i,k,m) = hygro(i,k,m)/dryvolmr(i,k);
133 hygro(i,k,m) = spechygro_1;
137 v2ncur_a = 1. / ( (
PI/6.)*(std::pow(dgncur_a(i,k,m),3.))*std::exp(4.5*std::pow(alnsg, 2.)) );
139 naer(i,k,m) = dryvolmr(i,k)*v2ncur_a;
146 if (maer(i,k,m) > 1.0e-31) {
147 drydens = maer(i,k,m)/dryvolmr(i,k);
152 dryvol(i,k,m) = 1.0/v2ncur_a;
153 drymass(i,k,m) = drydens*dryvol(i,k,m);
154 dryrad(i,k,m) = std::pow((dryvol(i,k,m)/
pi43),
third);
160 bool has_clear_rh =
false;
164 for(
auto k = top_lev; k <= nlev; --k) {
165 for(
auto i = 1; i <= ncol; ++i) {
166 rh(i,k) = clear_rh_in(i,k);
168 amrex::Print() <<
"modal_aero_wateruptake_dr: clear_rh_in is negative - rh:" << rh(i,k);
172 rh(i,k) = std::min(rh(i,k), 0.98);
177 for(
auto k = top_lev; k <= nlev; ++k) {
178 for(
auto i = 1; i <= ncol; ++i) {
180 if (qs(i) > h2ommr(i,k)) {
181 rh(i,k) = h2ommr(i,k)/qs(i);
186 rh(i,k) = std::max(rh(i,k), 0.0);
187 rh(i,k) = std::min(rh(i,k), 0.98);
188 if (cldn(i,k) < 0.998) {
189 rh(i,k) = (rh(i,k) - cldn(i,k)) / (1.0 - cldn(i,k));
191 rh(i,k) = std::max(rh(i,k), 0.0);
198 rhcrystal, rhdeliques, dryrad,
200 wetrad, wetvol, wtrvol);
202 for(
auto m = 1; m <= nmodes; ++m) {
203 for(
auto k = top_lev; k <= nlev; ++k) {
204 for(
auto i = 1; i <= ncol; ++i) {
206 dgncur_awet(i,k,m) = dgncur_a(i,k,m) * (wetrad(i,k,m)/dryrad(i,k,m));
207 qaerwat(i,k,m) =
rhoh2o*naer(i,k,m)*wtrvol(i,k,m);
210 if(compute_wetdens) {
211 if (wetvol(i,k,m) > 1.0e-30) {
212 wetdens(i,k,m) = (drymass(i,k,m) +
rhoh2o*wtrvol(i,k,m))/wetvol(i,k,m);
215 wetdens(i,k,m) = specdens_1(m);
226 const real1d& rhcrystal,
const real1d& rhdeliques,
const real3d& dryrad,
227 const real3d& hygro,
const real2d& rh,
const real3d& dryvol,
228 const real3d& wetrad,
const real3d& wetvol,
const real3d& wtrvol)
231 for(
auto m = 1; m <= nmodes; ++m) {
232 real hystfac = 1.0 / std::max(1.0e-5, (rhdeliques(m) - rhcrystal(m)));
233 for(
auto k = top_lev; k <= nlev; ++k) {
234 for(
auto i = 1; i <= ncol; ++i) {
237 modal_aero_kohler(dryrad(i,k,m), hygro(i,k,m), rh(i,k), wetrad(i,k,m));
239 wetrad(i,k,m) = std::max(wetrad(i,k,m), dryrad(i,k,m));
240 wetvol(i,k,m) =
pi43*std::pow(wetrad(i,k,m), 3.);
241 wetvol(i,k,m) = std::max(wetvol(i,k,m), dryvol(i,k,m));
242 wtrvol(i,k,m) = wetvol(i,k,m) - dryvol(i,k,m);
243 wtrvol(i,k,m) = std::max(wtrvol(i,k,m), 0.0);
248 if (rh(i,k) < rhcrystal(m)) {
249 wetrad(i,k,m) = dryrad(i,k,m);
250 wetvol(i,k,m) = dryvol(i,k,m);
253 else if (rh(i,k) < rhdeliques(m)) {
254 wtrvol(i,k,m) = wtrvol(i,k,m)*hystfac*(rh(i,k) - rhcrystal(m));
255 wtrvol(i,k,m) = std::max(wtrvol(i,k,m), 0.0);
256 wetvol(i,k,m) = dryvol(i,k,m) + wtrvol(i,k,m);
257 wetrad(i,k,m) = std::pow(wetvol(i,k,m)/
pi43, 1./3.);
271 const real& hygro,
const real& s,
274 const real eps = 1.e-4;
276 const real rhow = 1.;
277 const real surften = 76.;
278 const real tair = 273.;
279 const real
third = 1./3.;
280 const real ugascon = 8.3e7;
283 auto a=2.e4*mw*surften/(ugascon*tair*rhow);
285 auto rdry = rdry_in*1.0e6;
286 auto vol = std::pow(rdry, 3.);
290 auto ss = std::max(std::min(s,1.-eps), 1.e-10);
291 auto slog = std::log(ss);
294 auto p41 = b/slog-vol;
295 auto p40 = a*vol/slog;
308 auto p = std::abs(p31)/(rdry*rdry);
311 r=rdry*(1.+p*
third/(1.-slog*rdry/a));
314 amrex::GpuComplex<real> cx4[4];
319 for(
int n=0; n<4; ++n) {
320 auto xr=cx4[n].real();
321 auto xi=cx4[n].imag();
322 if(abs(xi) > abs(xr)*eps)
continue;
324 if(xr < rdry*(1.-eps))
continue;
325 if(xr != xr)
continue;
330 amrex::Print() <<
"kohlerc: no real solution found(quartic), nsol= " << nsol <<
"\n";
339 p=abs(p31)/(rdry*rdry);
344 amrex::GpuComplex<real> cx3[3];
349 for(
auto n=0; n<3; ++n) {
350 auto xr = cx3[n].real();
351 auto xi = cx3[n].imag();
352 if(abs(xi) > abs(xr)*eps)
continue;
354 if(xr < rdry*(1.-eps))
continue;
355 if(xr != xr)
continue;
360 amrex::Print() <<
"kohlerc: no real solution found (cubic), nsol= " << nsol <<
"\n";
366 r=(r4*(1.-s)+r3*(s-1.+eps))/eps;
381 const real eps = 1.0e-20;
382 const real
third=1./3.;
384 auto ci = amrex::GpuComplex<real>(0., 1.);
385 auto sqrt3 = amrex::GpuComplex<real>(std::sqrt(3.), 0.);
386 auto cw = 0.5*(-1.+ci*sqrt3);
387 auto cwsq = 0.5*(-1.-ci*sqrt3);
391 cx[0] = amrex::GpuComplex<real>(std::pow(-p0,
third), 0.);
396 auto q = amrex::GpuComplex<real>(p1/3., 0.);
397 auto r = amrex::GpuComplex<real>(p0/2., 0.);
398 auto crad = r*r+q*q*q;
399 crad = amrex::sqrt(crad);
402 if (amrex::abs(cy) > eps) cy = amrex::pow(cy,
third);
407 cx[1] = -cw*cy-cwsq*cz;
408 cx[2] = -cwsq*cy-cw*cz;
421 const real
third=1./3.;
422 auto q = -p2*p2/36.+(p3*p1-4*p0)/12.;
423 auto r = -std::pow((p2/6), 3.0) + p2*(p3*p1-4.*p0)/48.
424 + (4.*p0*p2-p0*p3*p3-p1*p1)/16.;
426 auto crad = amrex::sqrt(amrex::GpuComplex<real>(r*r+q*q*q, 0.));
429 if (cb.real() == 0. && cb.imag() == 0.) {
431 cx[0] = amrex::pow(amrex::GpuComplex<real>(-p1, 0.0),
third);
437 cb = amrex::pow(cb,
third);
438 auto cy = -cb+q/cb+p2/6.;
439 auto cb0 = amrex::sqrt(cy*cy-p0);
440 auto cb1 = (p3*cy-p1)/(2.*cb0);
443 crad = cb*cb-4.*(cy+cb0);
444 crad = amrex::sqrt(crad);
445 cx[1] = (-cb+crad)/2.;
446 cx[2] = (-cb-crad)/2.;
448 crad = cb*cb-4.*(cy-cb0);
449 crad = amrex::sqrt(crad);
450 cx[3] = (-cb+crad)/2.;
451 cx[4] = (-cb-crad)/2.;
constexpr amrex::Real rhoh2o
Definition: ERF_Constants.H:94
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
Definition: ERF_Mam4_constituents.H:18
void get_mode_nspec(int list_idx, int m_idx, int &nspec) const
Definition: ERF_Mam4_constituents.H:394
void rad_cnst_get_mam_mmr_by_idx(int list_idx, int mode_idx, int spec_idx, const std::string &phase, real2d &mmr) const
Definition: ERF_Mam4_constituents.H:639
void get_mam_density_aer(int list_idx, int mode_idx, int spec_idx, real &density_aer) const
Definition: ERF_Mam4_constituents.H:889
void get_mam_hygro_aer(int list_idx, int mode_idx, int spec_idx, real &hygro_aer) const
Definition: ERF_Mam4_constituents.H:903
void get_mode_props(int list_idx, int mode_idx, real &sigmag, real &rhcrystal, real &rhdeliques) const
Definition: ERF_Mam4_constituents.H:808
Definition: ERF_Modal_aero_wateruptake.H:18
static void makoh_cubic(amrex::GpuComplex< real > cx[], const real &p2, const real &p1, const real &p0)
Definition: ERF_Modal_aero_wateruptake.H:376
constexpr static real pi43
Definition: ERF_Modal_aero_wateruptake.H:21
constexpr static int imax
Definition: ERF_Modal_aero_wateruptake.H:23
static void makoh_quartic(amrex::GpuComplex< real > cx[], const real &p3, const real &p2, const real &p1, const real &p0)
Definition: ERF_Modal_aero_wateruptake.H:415
constexpr static real third
Definition: ERF_Modal_aero_wateruptake.H:20
constexpr static real huge_real
Definition: ERF_Modal_aero_wateruptake.H:22
static void modal_aero_wateruptake_sub(int ncol, int nlev, int nmodes, int top_lev, const real1d &rhcrystal, const real1d &rhdeliques, const real3d &dryrad, const real3d &hygro, const real2d &rh, const real3d &dryvol, const real3d &wetrad, const real3d &wetvol, const real3d &wtrvol)
Definition: ERF_Modal_aero_wateruptake.H:225
static void modal_aero_wateruptake_dr(int list_idx, int ncol, int nlev, int nmodes, int top_lev, MamConstituents &consti, const real2d &h2ommr, const real2d &t, const real2d &pmid, const real3d &dgnumdry_m, const real3d &dgnumwet_m, const real3d &qaerwat_m, const real3d &wetdens_m, const real2d &clear_rh_in)
Definition: ERF_Modal_aero_wateruptake.H:26
static YAKL_INLINE void modal_aero_kohler(const real &rdry_in, const real &hygro, const real &s, real &rwet_out)
Definition: ERF_Modal_aero_wateruptake.H:270
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE void qsat(const amrex::Real &t, const amrex::Real &p, amrex::Real &es, amrex::Real &qs, amrex::Real *gam=nullptr, amrex::Real *dqsdt=nullptr, amrex::Real *enthalpy=nullptr)
Definition: ERF_Water_vapor_saturation.H:152