Skip to content

Commit

Permalink
Radiation SW w/ moisture (#1597)
Browse files Browse the repository at this point in the history 
* Getting close with moisture.

* This runs a step.
  • Loading branch information
@AMLattanzi
AMLattanzi authored May 2, 2024
1 parent 5da4c3f commit 571acd4
Show file tree
Hide file tree
Showing 6 changed files with 46 additions and 57 deletions.
5 changes: 3 additions & 2 deletions Source/ERF_Constants.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,8 +66,9 @@ constexpr amrex::Real mwwv = 18.016;
constexpr amrex::Real lcond = 2.501e6;
constexpr amrex::Real lfus = 2.11727e3;
constexpr amrex::Real lsub = lcond+lfus;
constexpr amrex::Real rgas = boltz*avogadro/mwdair;
constexpr amrex::Real rv = rgas/mwwv;
constexpr amrex::Real rair = boltz*avogadro/mwdair;
constexpr amrex::Real rh20 = rair/mwwv;
constexpr amrex::Real rga = 1.0/CONST_GRAV;

constexpr amrex::Real diffelq = 2.21e-05; // Diffusivity of water vapor, m2/s
constexpr amrex::Real therco = 2.40e-02; // Thermal conductivity of air, J/m/s/K
Expand Down
12 changes: 6 additions & 6 deletions Source/Microphysics/SAM/Init_SAM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -230,8 +230,8 @@ void SAM::Compute_Coefficients ()
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);
coef1 = (lsub/(tabs1d_t(k)*rh20)-1.0)*lsub/(therco*rrr1*tabs1d_t(k));
coef2 = rh20*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);
Expand All @@ -243,8 +243,8 @@ void SAM::Compute_Coefficients ()
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);
coef1 = (lsub/(tabs1d_t(k)*rh20)-1.0)*lsub/(therco*rrr1*tabs1d_t(k));
coef2 = rh20*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);
Expand All @@ -253,8 +253,8 @@ void SAM::Compute_Coefficients ()
accrrc_t(k) = 0.25 * PI * nzeror * a_rain * gamr1 * pratio/pow((PI * rhor * nzeror / rho1d_t(k)) , ((3.0+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);
coef1 = (lcond/(tabs1d_t(k)*rh20)-1.0)*lcond/(therco*rrr1*tabs1d_t(k));
coef2 = rh20*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))/
Expand Down
12 changes: 6 additions & 6 deletions Source/Radiation/Aero_rad_props.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,8 @@ void AerRadProps::initialize (int num_gas, int num_modes, int naeroes,
nrh = num_rh;
top_lev = top_levels;

pmid = pmiddle;
// TODO: Interpretation for pdeldry (it is NOT absolute pressure)!
pmid = pmiddle;
pdeldry = pmiddle;
temp = temperature;
qt = qtotal;
Expand All @@ -52,12 +53,12 @@ void AerRadProps::aer_rad_props_sw (const int& list_idx, const real& dt, const i
real2d mmr_to_mass("mmr_to_mass", ncol, nlev); // conversion factor for mmr to mass
parallel_for(SimpleBounds<2>(ncol, nlev) , YAKL_LAMBDA (int icol, int ilev)
{
mmr_to_mass(icol,ilev) = rgas * pdeldry(icol,ilev);
mmr_to_mass(icol,ilev) = rga * pdeldry(icol,ilev);
});

// initialize to conditions that would cause failure
yakl::memset(tau, -100.);
yakl::memset(tau_w, -100.);
yakl::memset(tau , -100.);
yakl::memset(tau_w , -100.);
yakl::memset(tau_w_g, -100.);
yakl::memset(tau_w_f, -100.);
yakl::memset(ext_cmip6_sw_inv_m, 1.0e-3);
Expand Down Expand Up @@ -125,7 +126,6 @@ void AerRadProps::aer_rad_props_sw (const int& list_idx, const real& dt, const i
if (tlev != -1) trop_level(i) = tlev;
});

//
//Quit if tropopause is not found
if (yakl::intrinsics::any(trop_level) == -1) {
amrex::Print() << "aer_rad_props.F90: subr aer_rad_props_sw: tropopause not found\n";
Expand Down Expand Up @@ -350,7 +350,7 @@ void AerRadProps::aer_rad_props_lw (const bool& is_cmip6_volc,
// compute mixing ratio to mass conversion
parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
{
mmr_to_mass(icol,ilev) = rgas * pdeldry(icol,ilev);
mmr_to_mass(icol,ilev) = rga * pdeldry(icol,ilev);
});

// calculate relative humidity for table lookup into rh grid
Expand Down
26 changes: 14 additions & 12 deletions Source/Radiation/Mam4_aero.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -387,8 +387,8 @@ class Mam4_aer {

parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k)
{
mass(i,k) = pdeldry(i,k)*rgas;
air_density(i,k) = pmid(i,k)/(rv*temperature(i,k));
mass(i,k) = pdeldry(i,k)*rga;
air_density(i,k) = pmid(i,k)/(rair*temperature(i,k));
});

// Calculate aerosol size distribution parameters and aerosol water uptake
Expand Down Expand Up @@ -429,8 +429,8 @@ class Mam4_aer {
for(auto k = top_lev; k <= nlev; ++k) {
// form bulk refractive index
yakl::memset(crefin_real, 0.);
yakl::memset(crefin_im, 0.);
yakl::memset(dryvol, 0.);
yakl::memset(crefin_im , 0.);
yakl::memset(dryvol , 0.);

// aerosol species loop
for(auto l = 0; l < nspec; ++l) {
Expand Down Expand Up @@ -603,7 +603,7 @@ class Mam4_aer {
// dry mass in each cell
parallel_for (SimpleBounds<2> (ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
{
//mass(icol,ilev) = pdeldry(icol,ilev)*rga;
mass(icol,ilev) = pdeldry(icol,ilev)*rga;
});

// Calculate aerosol size distribution parameters and aerosol water uptake
Expand Down Expand Up @@ -638,7 +638,7 @@ class Mam4_aer {
// some intermediate results are saved and the chebyshev polynomials are stored
// in a array with different index order. Could be unified.
top_lev = 1;
for(auto k = top_lev; k > nlev; ++k) {
for(auto k = top_lev; k <= nlev; ++k) {
for(auto i = 1; i <= ncol; ++i) {
alnsg_amode = std::log(sigma_logr_aer);
// convert from number diameter to surface area
Expand All @@ -650,15 +650,14 @@ class Mam4_aer {
// chebyshev polynomials
cheby(1,i,k) = 1.0;
cheby(2,i,k) = xrad(i);
for(auto nc = 3; nc < ncoef; ++nc) {
for(auto nc = 3; nc <= ncoef; ++nc) {
cheby(nc,i,k) = 2.0*xrad(i)*cheby(nc-1,i,k)-cheby(nc-2,i,k);
}
}
}

//int nlwbands; // add by xyuan to be compiled
for(auto ilw = 1; ilw <= nlwbands; ++ilw) {
for(auto k = top_lev; k > nlev; --k) {
for(auto k = top_lev; k <= nlev; ++k) {
// form bulk refractive index. Use volume mixing for infrared
yakl::memset(crefinr, 0.0);
yakl::memset(crefini, 0.0);
Expand All @@ -669,7 +668,7 @@ class Mam4_aer {
mam_consti.get_mam_props_lw(list_idx, m, l, specdens, specrefrindex, specrefiindex);
for(auto i = 1; i <= ncol; ++i) {
vol(i) = specmmr(i,k)/specdens;
dryvol(i) = dryvol(i) + vol(i);
dryvol(i) = dryvol(i) + vol(i);
crefinr(i) = crefinr(i) + vol(i)*specrefrindex(ilw);
crefini(i) = crefini(i) + vol(i)*specrefiindex(ilw);
}
Expand Down Expand Up @@ -715,7 +714,7 @@ class Mam4_aer {
// parameterized optical properties
for(auto i = 1; i <= ncol; ++i) {
pabs(i) = 0.5*cabs(i,1);
for(auto nc = 2; nc < ncoef; ++nc) {
for(auto nc = 2; nc <= ncoef; ++nc) {
pabs(i) = pabs(i) + cheby(nc,i,k)*cabs(i,nc);
}
pabs(i) = pabs(i)*wetvol(i)*rhoh2o;
Expand Down Expand Up @@ -743,7 +742,10 @@ class Mam4_aer {
}
}
}
for(auto i = 1; i <= ncol; ++i) tauxar(i,k,ilw) = tauxar(i,k,ilw) + dopaer(i);

for(auto i = 1; i <= ncol; ++i) {
tauxar(i,k,ilw) = tauxar(i,k,ilw) + dopaer(i);
}
} // k = top_lev, nlev
} // nlwbands
} // m = 1, nmodes
Expand Down
9 changes: 0 additions & 9 deletions Source/Radiation/Radiation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,19 +187,10 @@ void Radiation::initialize (const MultiFab& cons_in,
auto nx = box3d.length(0);

auto states_array = cons_in.array(mfi);
/*
auto qt_array = (qmoist[0]) ? qmoist[0]->array(mfi) : Array4<Real> {};
auto qv_array = (qmoist[1]) ? qmoist[1]->array(mfi) : Array4<Real> {};
auto qc_array = (qmoist[2]) ? qmoist[2]->array(mfi) : Array4<Real> {};
auto qi_array = (qmoist.size()>=8) ? qmoist[3]->array(mfi) : Array4<Real> {};
*/

// DEBUG: delete and uncomment when done
// NOTE : SW crashes with moisture but runs without it (pressure unit issue?!)
auto qt_array = (false) ? qmoist[0]->array(mfi) : Array4<Real> {};
auto qv_array = (false) ? qmoist[1]->array(mfi) : Array4<Real> {};
auto qc_array = (false) ? qmoist[2]->array(mfi) : Array4<Real> {};
auto qi_array = (false) ? qmoist[3]->array(mfi) : Array4<Real> {};

// Get pressure, theta, temperature, density, and qt, qp
ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
Expand Down
39 changes: 17 additions & 22 deletions Source/Radiation/Run_shortwave_rrtmgp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,31 +47,23 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay,
toa_flux(icol, igpt) = tsi_scaling * toa_flux(icol,igpt);
});

// Add in aerosol
// TODO: should we avoid allocating here?
// Add in aerosol, allocate on gpts and map aer_*_bnd from bands
OpticalProps2str aerosol_optics;
auto &aerosol_optics_tau = aerosol_optics.tau;
auto &aerosol_optics_ssa = aerosol_optics.ssa;
auto &aerosol_optics_g = aerosol_optics.g ;
if (true) {
aerosol_optics.alloc_2str(ncol, nlay, k_dist_sw);
auto gpt_bnd = aerosol_optics.get_gpoint_bands();
parallel_for(SimpleBounds<3>(nswgpts,nlay,ncol) , YAKL_LAMBDA (int igpt, int ilay, int icol) {
aerosol_optics_tau(icol,ilay,igpt) = aer_tau_bnd(icol,ilay,gpt_bnd(igpt));
aerosol_optics_ssa(icol,ilay,igpt) = aer_ssa_bnd(icol,ilay,gpt_bnd(igpt));
aerosol_optics_g (icol,ilay,igpt) = aer_asm_bnd(icol,ilay,gpt_bnd(igpt));
});
} else {
aerosol_optics.alloc_2str(ncol, nlay, k_dist_sw.get_band_lims_wavenumber());
parallel_for(SimpleBounds<3>(nswbands,nlay,ncol), YAKL_LAMBDA (int ibnd, int ilay, int icol) {
aerosol_optics_tau(icol,ilay,ibnd) = aer_tau_bnd(icol,ilay,ibnd);
aerosol_optics_ssa(icol,ilay,ibnd) = aer_ssa_bnd(icol,ilay,ibnd);
aerosol_optics_g (icol,ilay,ibnd) = aer_asm_bnd(icol,ilay,ibnd);
});
}
aerosol_optics.alloc_2str(ncol, nlay, k_dist_sw);
auto gpt_bnd = aerosol_optics.get_gpoint_bands();
parallel_for(SimpleBounds<3>(nswgpts,nlay,ncol) , YAKL_LAMBDA (int igpt, int ilay, int icol)
{
aerosol_optics.tau(icol,ilay,igpt) = aer_tau_bnd(icol,ilay,gpt_bnd(igpt));
aerosol_optics.ssa(icol,ilay,igpt) = aer_ssa_bnd(icol,ilay,gpt_bnd(igpt));
aerosol_optics.g (icol,ilay,igpt) = aer_asm_bnd(icol,ilay,gpt_bnd(igpt));
});

aerosol_optics.delta_scale();

// NOTE: aero_optics is allocated on nswgpts and combined_optics is on nswgpts
// The `increment` call below can handle matching or differing (nswgpts/nswbnds)
// sizes; see calls in `mo_optical_props_kernels.cpp`. Since we have matching
// gpt sizes, we will call `increment_2stream_by_2stream`
aerosol_optics.increment(combined_optics);

// Do the clearsky calculation before adding in clouds
Expand All @@ -97,7 +89,7 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay,
fluxes_clrsky.bnd_flux_net.deep_copy_to(clrsky_bnd_flux_net);
fluxes_clrsky.bnd_flux_dn_dir.deep_copy_to(clrsky_bnd_flux_dn_dir);

// Add in clouds
// Add in clouds, which are already on gpts
OpticalProps2str cloud_optics;
cloud_optics.alloc_2str(ncol, nlay, k_dist_sw);
auto &cloud_optics_tau = cloud_optics.tau;
Expand All @@ -108,7 +100,10 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay,
cloud_optics_ssa(icol,ilay,igpt) = cld_ssa_gpt(icol,ilay,igpt);
cloud_optics_g (icol,ilay,igpt) = cld_asm_gpt(icol,ilay,igpt);
});

cloud_optics.delta_scale();

// NOTE: See above and here we again have matching gpt sizes
cloud_optics.increment(combined_optics);

// Call SW flux driver
Expand Down

0 comments on commit 571acd4

Please sign in to comment.