Skip to content

Commit

Permalink
More Radiation (#1587)
Browse files Browse the repository at this point in the history 
* Fix heating rate calculation and add qv to the EOS calls.

* Fix units and pressure/temp extrap. LW fluxes look weird still.
  • Loading branch information
@AMLattanzi
AMLattanzi authored Apr 23, 2024
1 parent 0062ec9 commit 94815bb
Show file tree
Hide file tree
Showing 4 changed files with 74 additions and 59 deletions.
2 changes: 1 addition & 1 deletion Source/ERF_make_new_arrays.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -246,7 +246,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
//*********************************************************
// Radiation heating source terms
//*********************************************************
qheating_rates[lev] = std::make_unique<MultiFab>(ba, dm, 1, ngrow_state);
qheating_rates[lev] = std::make_unique<MultiFab>(ba, dm, 2, ngrow_state);
qheating_rates[lev]->setVal(0.);
#endif
}
Expand Down
68 changes: 39 additions & 29 deletions Source/Radiation/Radiation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,6 +177,7 @@ void Radiation::initialize (const MultiFab& cons_in,
for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) {
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> {};

Expand All @@ -188,26 +189,32 @@ void Radiation::initialize (const MultiFab& cons_in,
{
auto icol = j*nx+i+1;
auto ilev = k+1;
Real qv = (qv_array) ? qv_array(i,j,k): 0.0;
qt(icol,ilev) = (qt_array) ? qt_array(i,j,k): 0.0;
qc(icol,ilev) = (qc_array) ? qc_array(i,j,k): 0.0;
qi(icol,ilev) = (qi_array) ? qi_array(i,j,k): 0.0;
qn(icol,ilev) = qc(icol,ilev) + qi(icol,ilev);
tmid(icol,ilev) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp));
pmid(icol,ilev) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/1.0e3;
tmid(icol,ilev) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp),qv);
// NOTE: RRTMGP code expects pressure in mb so we convert it here
pmid(icol,ilev) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp),qv) * 1.0e-2;
});
}

parallel_for(SimpleBounds<2>(ncol, nlev+1), YAKL_LAMBDA (int icol, int ilev)
{
if (ilev == 1) {
pint(icol, 1) = 2.*pmid(icol, 2) - pmid(icol, 1);
tint(icol, 1) = 2.*tmid(icol, 2) - tmid(icol, 1);
//pint(icol, 1) = 2.*pmid(icol, 2) - pmid(icol, 1);
//tint(icol, 1) = 2.*tmid(icol, 2) - tmid(icol, 1);
pint(icol, 1) = -0.5*pmid(icol, 2) + 1.5*pmid(icol, 1);
tint(icol, 1) = -0.5*tmid(icol, 2) + 1.5*tmid(icol, 1);
} else if (ilev <= nlev) {
pint(icol, ilev) = 0.5*(pmid(icol, ilev-1) + pmid(icol, ilev));
tint(icol, ilev) = 0.5*(tmid(icol, ilev-1) + tmid(icol, ilev));
} else {
pint(icol, nlev+1) = 2.*pmid(icol, nlev-1) - pmid(icol, nlev);
tint(icol, nlev+1) = 2.*tmid(icol, nlev-1) - tmid(icol, nlev);
//pint(icol, nlev+1) = 2.*pmid(icol, nlev-1) - pmid(icol, nlev);
//tint(icol, nlev+1) = 2.*tmid(icol, nlev-1) - tmid(icol, nlev);
pint(icol, nlev+1) = -0.5*pmid(icol, nlev-1) + 1.5*pmid(icol, nlev);
tint(icol, nlev+1) = -0.5*tmid(icol, nlev-1) + 1.5*tmid(icol, nlev);
}
});

Expand Down Expand Up @@ -252,9 +259,6 @@ void Radiation::initialize (const MultiFab& cons_in,
// run radiation model
void Radiation::run ()
{
// Temporary variable for heating rate output
real2d hr("hr", ncol, nlev);

// Cosine solar zenith angle for all columns in chunk
real1d coszrs("coszrs", ncol);

Expand Down Expand Up @@ -512,32 +516,20 @@ void Radiation::run ()
}
} // dolw

// Compute heating rate for dtheta/dt
parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
{
// TODO: CHECK THIS UNIT CONVERSION!!
hr(icol,ilev) = (qrs(icol,ilev) + qrl(icol,ilev)) / Cp_d * (1.e5 / std::pow(pmid(icol,ilev), R_d/Cp_d));
});

// Populate theta source from hr
// Populate source term for theta dycore variable
for (MFIter mfi(*(qrad_src)); mfi.isValid(); ++mfi) {
auto qrad_src_array = qrad_src->array(mfi);
const auto& box3d = mfi.tilebox();
auto nx = box3d.length(0);
amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
// Map (col,lev) to (i,j,k)
auto icol = j*nx+i+1;
auto ilev = k+1;
qrad_src_array(i,j,k) = hr(icol,ilev);
});
}

// convert radiative heating rates to Q*dp for energy conservation
if (conserve_energy) {
parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
{
qrs(icol,ilev) = qrs(icol,ilev)*pdel(icol,ilev);
qrl(icol,ilev) = qrl(icol,ilev)*pdel(icol,ilev);
// SW and LW sources
qrad_src_array(i,j,k,0) = qrs(icol,ilev);
qrad_src_array(i,j,k,1) = qrl(icol,ilev);
});
}
}
Expand Down Expand Up @@ -894,11 +886,29 @@ void Radiation::calculate_heating_rate (const real2d& flux_up,
const real2d& pint,
const real2d& heating_rate)
{
// NOTE: The pressure is in [mb] for RRTMGP to use.
// The fluxes are in [W/m^2] and gravity is [m/s^2].
// We need to convert pressure from [mb] -> [Pa] and divide by Cp [J/kg*K]
// The heating rate is {dF/dP * g / Cp} with units [K/s]
real1d heatfac("heatfac",1);
yakl::memset(heatfac, 1.0e-2/Cp_d);
parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
{
heating_rate(icol,ilev) = ( (flux_up(icol,ilev+1) - flux_dn(icol,ilev+1))
- (flux_up(icol,ilev ) - flux_dn(icol,ilev )) )
* CONST_GRAV/(pint(icol,ilev+1)-pint(icol,ilev));
heating_rate(icol,ilev) = heatfac(1) * ( (flux_up(icol,ilev+1) - flux_dn(icol,ilev+1))
- (flux_up(icol,ilev ) - flux_dn(icol,ilev )) )
* CONST_GRAV/(pint(icol,ilev+1)-pint(icol,ilev));
/*
if (icol==1) {
amrex::Print() << "HR: " << ilev << ' '
<< heating_rate(icol,ilev) << ' '
<< (flux_up(icol,ilev+1) - flux_dn(icol,ilev+1)) << ' '
<< (flux_up(icol,ilev ) - flux_dn(icol,ilev )) << ' '
<< (flux_up(icol,ilev+1) - flux_dn(icol,ilev+1))
- (flux_up(icol,ilev ) - flux_dn(icol,ilev )) << ' '
<< (pint(icol,ilev+1)-pint(icol,ilev)) << ' '
<< CONST_GRAV << "\n";
}
*/
});
}

Expand Down
60 changes: 32 additions & 28 deletions Source/Radiation/Run_shortwave_rrtmgp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,16 +34,18 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay,
bool1d top_at_1_g("top_at_1_g",1);
boolHost1d top_at_1_h("top_at_1_h",1);
bool top_at_1;
parallel_for(SimpleBounds<1>(1), YAKL_LAMBDA (int icol) { // HACK: Single loop kernel is not efficient
top_at_1_g(1) = pmid(1, 1) < pmid (1, 2);
});
top_at_1_g.deep_copy_to(top_at_1_h);
top_at_1 = top_at_1_h(1);
real2d toa_flux("toa_flux", ncol, nswgpts);
k_dist_sw.gas_optics(ncol, nlay, &top_at_1, pmid, pint, tmid, gas_concs, combined_optics, toa_flux); // TODO: top_at_1 is not a valid address and this doesn't match longwave call
k_dist_sw.gas_optics(ncol, nlay, top_at_1, pmid, pint, tmid, gas_concs, combined_optics, toa_flux);

// Apply TOA flux scaling
parallel_for(SimpleBounds<2>(nswgpts,ncol), YAKL_LAMBDA (int igpt, int icol) {
toa_flux(icol, igpt) = tsi_scaling * toa_flux(icol,igpt);
top_at_1_g(1) = pmid(1, 1) < pmid (1, 2);
});
top_at_1_g.deep_copy_to(top_at_1_h);
top_at_1 = top_at_1_h(1);

// Add in aerosol
// TODO: should we avoid allocating here?
Expand Down Expand Up @@ -74,25 +76,26 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay,

// Do the clearsky calculation before adding in clouds
FluxesByband fluxes_clrsky;
fluxes_clrsky.flux_up = real2d("clrsky_flux_up", ncol, nlay+1); // clrsky_flux_up;
fluxes_clrsky.flux_dn = real2d("clrsky_flux_up", ncol, nlay+1); //clrsky_flux_dn;
fluxes_clrsky.flux_dn_dir = real2d("clrsky_flux_up", ncol, nlay+1); //clrsky_flux_dn_dir;
fluxes_clrsky.flux_net = real2d("clrsky_flux_up", ncol, nlay+1); //clrsky_flux_net;
fluxes_clrsky.bnd_flux_up = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_up;
fluxes_clrsky.bnd_flux_dn = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn;
fluxes_clrsky.bnd_flux_dn_dir = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn_dir;
fluxes_clrsky.bnd_flux_net = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_net;
fluxes_clrsky.flux_up = real2d("clrsky_flux_up" , ncol, nlay+1); // clrsky_flux_up;
fluxes_clrsky.flux_dn = real2d("clrsky_flux_nd" , ncol, nlay+1); //clrsky_flux_dn;
fluxes_clrsky.flux_net = real2d("clrsky_flux_net", ncol, nlay+1); //clrsky_flux_net;
fluxes_clrsky.flux_dn_dir = real2d("clrsky_flux_dn_dir", ncol, nlay+1); //clrsky_flux_dn_dir;
fluxes_clrsky.bnd_flux_up = real3d("clrsky_bnd_flux_up" , ncol, nlay+1, nswbands); //clrsky_bnd_flux_up;
fluxes_clrsky.bnd_flux_dn = real3d("clrsky_bnd_flux_dn" , ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn;
fluxes_clrsky.bnd_flux_net = real3d("clrsky_bnd_flux_net", ncol, nlay+1, nswbands); //clrsky_bnd_flux_net;
fluxes_clrsky.bnd_flux_dn_dir = real3d("clrsky_bnd_flux_dn_dir", ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn_dir;

rte_sw(combined_optics, top_at_1, coszrs, toa_flux, albedo_dir, albedo_dif, fluxes_clrsky);

// Copy fluxes back out of FluxesByband object
fluxes_clrsky.flux_up.deep_copy_to(clrsky_flux_up);
fluxes_clrsky.flux_dn.deep_copy_to(clrsky_flux_dn);
fluxes_clrsky.flux_dn_dir.deep_copy_to(clrsky_flux_dn_dir);
fluxes_clrsky.flux_up.deep_copy_to (clrsky_flux_up);
fluxes_clrsky.flux_dn.deep_copy_to (clrsky_flux_dn);
fluxes_clrsky.flux_net.deep_copy_to(clrsky_flux_net);
fluxes_clrsky.bnd_flux_up.deep_copy_to(clrsky_bnd_flux_up);
fluxes_clrsky.bnd_flux_dn.deep_copy_to(clrsky_bnd_flux_dn);
fluxes_clrsky.bnd_flux_dn_dir.deep_copy_to(clrsky_bnd_flux_dn_dir);
fluxes_clrsky.flux_dn_dir.deep_copy_to(clrsky_flux_dn_dir);
fluxes_clrsky.bnd_flux_up.deep_copy_to (clrsky_bnd_flux_up);
fluxes_clrsky.bnd_flux_dn.deep_copy_to (clrsky_bnd_flux_dn);
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
OpticalProps2str cloud_optics;
Expand All @@ -110,25 +113,26 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay,

// Call SW flux driver
FluxesByband fluxes_allsky;
fluxes_allsky.flux_up = real2d("allsky_flux_up", ncol, nlay+1); // allsky_flux_up;
fluxes_allsky.flux_dn = real2d("allsky_flux_up", ncol, nlay+1); //allsky_flux_dn;
fluxes_allsky.flux_dn_dir = real2d("allsky_flux_up", ncol, nlay+1); //allsky_flux_dn_dir;
fluxes_allsky.flux_net = real2d("allsky_flux_up", ncol, nlay+1); //allsky_flux_net;
fluxes_allsky.bnd_flux_up = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_up;
fluxes_allsky.bnd_flux_dn = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_dn;
fluxes_allsky.bnd_flux_dn_dir = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_dn_dir;
fluxes_allsky.bnd_flux_net = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_net;
fluxes_allsky.flux_up = real2d("allsky_flux_up" , ncol, nlay+1); //allsky_flux_up;
fluxes_allsky.flux_dn = real2d("allsky_flux_dn" , ncol, nlay+1); //allsky_flux_dn;
fluxes_allsky.flux_net = real2d("allsky_flux_net", ncol, nlay+1); //allsky_flux_net;
fluxes_allsky.flux_dn_dir = real2d("allsky_flux_dn_dir", ncol, nlay+1); //allsky_flux_dn_dir;
fluxes_allsky.bnd_flux_up = real3d("allsky_bnd_flux_up" , ncol, nlay+1, nswbands); //allsky_bnd_flux_up;
fluxes_allsky.bnd_flux_dn = real3d("allsky_bnd_flux_dn" , ncol, nlay+1, nswbands); //allsky_bnd_flux_dn;
fluxes_allsky.bnd_flux_net = real3d("allsky_bnd_flux_net", ncol, nlay+1, nswbands); //allsky_bnd_flux_net;
fluxes_allsky.bnd_flux_dn_dir = real3d("allsky_bnd_flux_dn_dir", ncol, nlay+1, nswbands); //allsky_bnd_flux_dn_dir;

rte_sw(combined_optics, top_at_1, coszrs, toa_flux, albedo_dir, albedo_dif, fluxes_allsky);

// Copy fluxes back out of FluxesByband object
fluxes_allsky.flux_up.deep_copy_to(allsky_flux_up);
fluxes_allsky.flux_dn.deep_copy_to(allsky_flux_dn);
fluxes_allsky.flux_dn_dir.deep_copy_to(allsky_flux_dn_dir);
fluxes_allsky.flux_net.deep_copy_to(allsky_flux_net);
fluxes_allsky.flux_dn_dir.deep_copy_to(allsky_flux_dn_dir);
fluxes_allsky.bnd_flux_up.deep_copy_to(allsky_bnd_flux_up);
fluxes_allsky.bnd_flux_dn.deep_copy_to(allsky_bnd_flux_dn);
fluxes_allsky.bnd_flux_dn_dir.deep_copy_to(allsky_bnd_flux_dn_dir);
fluxes_allsky.bnd_flux_net.deep_copy_to(allsky_bnd_flux_net);
fluxes_allsky.bnd_flux_dn_dir.deep_copy_to(allsky_bnd_flux_dn_dir);
yakl::fence();
}

3 changes: 2 additions & 1 deletion Source/SourceTerms/ERF_make_sources.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -153,7 +153,8 @@ void make_sources (int level,
auto const& qheating_arr = qheating_rates->const_array(mfi);
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
cell_src(i,j,k,RhoTheta_comp) += qheating_arr(i,j,k);
// Short-wavelength and long-wavelength radiation source terms
cell_src(i,j,k,RhoTheta_comp) += qheating_arr(i,j,k,0) + qheating_arr(i,j,k,1);
});
}

Expand Down

0 comments on commit 94815bb

Please sign in to comment.