Skip to content

Commit

Permalink
Merge pull request #1404 from xyuan/xyuan/rrtmgp
Browse files Browse the repository at this point in the history 
fix some bugs in the aerosol tracer's mass fraction, and cloud species mass fraction to match e3sm
  • Loading branch information
@xyuan
xyuan authored Jan 26, 2024
2 parents 8609801 + dd79be6 commit 02493b8
Show file tree
Hide file tree
Showing 8 changed files with 281 additions and 381 deletions.
10 changes: 5 additions & 5 deletions Source/Radiation/Aero_rad_props.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,12 @@ class AerRadProps {
~AerRadProps() = default;

// initialization
void initialize(int ngas_, int nmodes_, int num_aeros_,
void initialize(int num_gas, int num_modes, int naeroes,
int nswbands_, int nlwbands_,
int ncol_, int nlev_, int nrh_, int top_lev_,
const std::vector<std::string>& aero_names_,
const real2d& zi_, const real2d& pmid_, const real2d& temp_,
const real2d& qi_, const real2d& geom_radius);
int ncoloum, int nlevel, int num_rh, int top_levels,
const std::vector<std::string>& aerosol_names,
const real2d& zint, const real2d& pmiddle, const real2d& temperature,
const real2d& qtotal, const real2d& geom_rad);

void aer_rad_props_sw(const int& list_idx,
const real& dt,
Expand Down
393 changes: 154 additions & 239 deletions Source/Radiation/Aero_rad_props.cpp
View file

Large diffs are not rendered by default.

70 changes: 28 additions & 42 deletions Source/Radiation/Mam4_aero.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,12 +53,13 @@ class Mam4_aer {
public:
// initialization
inline
void initialize(int ncol_, int nlev_, int top_lev_, int nswbands_, int nlwbands_) {
ncol = ncol_;
nlev = nlev_;
top_lev = top_lev_;
nswbands = nswbands_;
nlwbands = nlwbands_;
void initialize(int ncols, int nlevs, int top_levs, int nsw_bands, int nlw_bands)
{
ncol = ncols;
nlev = nlevs;
top_lev = top_levs;
nswbands = nsw_bands;
nlwbands = nlw_bands;
clim_modal_aero = true;

const real rmmin = 0.01e-6;
Expand All @@ -81,8 +82,7 @@ class Mam4_aer {
water_refindex_file = erf_rad_data_dir+"water_refindex_rrtmg_c080910.nc";
read_water_refindex(water_refindex_file);

// Allocate dry and wet size variables in an OMP PARRALLEL region as these
// arrays are private for each thread and needs to be allocated for each OMP thread
// Allocate dry and wet size variables
dgnumdry_m = real3d("dgnumdry", ncol, nlev, nmodes);
dgnumwet_m = real3d("dgnumwet", ncol, nlev, nmodes);
qaerwat_m = real3d("qaerwat" , ncol, nlev, nmodes);
Expand Down Expand Up @@ -239,22 +239,13 @@ class Mam4_aer {
real2d specmmr("specmmr",ncol,nlev); // specie mmr
real2d dryvol_a("dryvol_a",ncol,nlev); // interstital aerosol dry volume (cm^3/mol_air)

real dgnum, dgnumhi, dgnumlo;
real dgnyy, dgnxx; // dgnumlo/hi of current mode
real drv_a; // dry volume (cm3/mol_air)
real dumfac, dummwdens; // work variables
real num_a0; // initial number (#/mol_air)
real num_a; // final number (#/mol_air)
real voltonumbhi, voltonumblo;
real v2nmax, v2nmin; // voltonumblo/hi of current mode
real sigmag, alnsg;

for (auto n = 1; n <= nmodes; ++n) {
parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k) {
dgncur_a(i,k) = dgnum_m(i,k,n);
});

// get mode properties
real dgnum, dgnumhi, dgnumlo, sigmag;
mam_consti.get_mode_props(list_idx, n-1, dgnum, dgnumhi, dgnumlo, sigmag=sigmag);

// get mode number mixing ratio
Expand All @@ -272,7 +263,7 @@ class Mam4_aer {
mam_consti.get_mam_props(list_idx, n-1, l1, specdens);

// need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
dummwdens = 1.0 / specdens;
real dummwdens = 1.0 / specdens;
top_lev = 1;

for (auto k=top_lev; k <= nlev; ++k) {
Expand All @@ -282,21 +273,21 @@ class Mam4_aer {
}
}

alnsg = std::log(sigmag);
dumfac = std::exp(4.5*std::pow(alnsg, 2))*PI/6.0;
voltonumblo = 1./((PI/6.)*(std::pow(dgnumlo, 3))*std::exp(4.5*std::pow(alnsg, 2)));
voltonumbhi = 1./((PI/6.)*(std::pow(dgnumhi, 3))*std::exp(4.5*std::pow(alnsg, 2)));
v2nmin = voltonumbhi;
v2nmax = voltonumblo;
dgnxx = dgnumhi;
dgnyy = dgnumlo;
auto alnsg = std::log(sigmag);
auto dumfac = std::exp(4.5*std::pow(alnsg, 2))*PI/6.0;
auto voltonumblo = 1./((PI/6.)*(std::pow(dgnumlo, 3))*std::exp(4.5*std::pow(alnsg, 2)));
auto voltonumbhi = 1./((PI/6.)*(std::pow(dgnumhi, 3))*std::exp(4.5*std::pow(alnsg, 2)));
auto v2nmin = voltonumbhi;
auto v2nmax = voltonumblo;
auto dgnxx = dgnumhi;
auto dgnyy = dgnumlo;

top_lev = 1;
for (auto k = top_lev; k <= nlev; ++k) {
for (auto i = 1; i <= ncol; ++i) {
drv_a = dryvol_a(i,k);
num_a0 = mode_num(i,k);
num_a = std::max(0.0, num_a0);
auto drv_a = dryvol_a(i,k);
auto num_a0 = mode_num(i,k);
auto num_a = std::max(0.0, num_a0);

if (drv_a > 0.0) {
if (num_a <= drv_a*v2nmin)
Expand Down Expand Up @@ -371,21 +362,18 @@ class Mam4_aer {
real1d pasm("pasm",ncol); // parameterized asymmetry factor
real1d palb("palb",ncol); // parameterized single scattering albedo

// error code
int nerr_dopaer;

// initialize output variables
yakl::memset(tauxar, 0.);
yakl::memset(wa, 0.);
yakl::memset(ga, 0.);
yakl::memset(fa, 0.);

// zero'th layer does not contain aerosol
parallel_for(SimpleBounds<2>(ncol, nswbands), YAKL_LAMBDA (int i, int ibnd) {
wa(i,1,ibnd) = 0.925;
ga(i,1,ibnd) = 0.850;
fa(i,1,ibnd) = 0.7225;
});
// parallel_for(SimpleBounds<2>(ncol, nswbands), YAKL_LAMBDA (int i, int ibnd) {
// wa(i,1,ibnd) = 0.925;
// ga(i,1,ibnd) = 0.850;
// fa(i,1,ibnd) = 0.7225;
// });

parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k) {
mass(i,k) = pdeldry(i,k)*rgas;
Expand Down Expand Up @@ -464,7 +452,6 @@ class Mam4_aer {
refr(i) = crefin_real(i);
refi(i) = std::abs(crefin_im(i));
}

// interpolate coefficients linear in refractive index
// first call calcs itab,jtab,ttab,utab
real3d extpswr("extpswr", ncoef, prefr, prefi);
Expand Down Expand Up @@ -533,9 +520,8 @@ class Mam4_aer {
mam_consti.get_mam_props_sw(list_idx, m-1, l, specdens, specrefindex_real, specrefindex_im);
volf = specmmr(i,k)/specdens;
}
nerr_dopaer = nerr_dopaer + 1;
if (dopaer(i) < -1.e-10) {
amrex::Print() << "*** halting after nerr_dopaer\n";
amrex::Print() << "*** halting with error!\n";
exit(0);
}
}
Expand Down Expand Up @@ -637,7 +623,7 @@ class Mam4_aer {
top_lev = 1;
for(auto k = top_lev; k > nlev; ++k) {
for(auto i = 1; i <= ncol; ++i) {
alnsg_amode = std::log( sigma_logr_aer );
alnsg_amode = std::log(sigma_logr_aer);
// convert from number diameter to surface area
xrad(i) = std::log(0.5*dgnumwet(i,k)) + 2.0*alnsg_amode*alnsg_amode;
// normalize size parameter
Expand Down
55 changes: 28 additions & 27 deletions Source/Radiation/Radiation.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
#include <AMReX_TableData.H>
#include <AMReX_MultiFabUtil.H>

#include "ERF_Config.H"
#include "ERF_Constants.H"
#include "Rad_constants.H"
#include "Rrtmgp.H"
Expand All @@ -34,12 +35,12 @@
// Radiation code interface class
class Radiation {
public:
Radiation () {
Radiation() {
// First, make sure yakl has been initialized
if (!yakl::isInitialized()) yakl::init();
}

~Radiation () = default;
~Radiation() = default;

// init
void initialize(const amrex::MultiFab& cons_in,
Expand All @@ -54,37 +55,37 @@ class Radiation {
const bool& is_cmip6_volcano);

// run radiation model
void run ();
void run();

// call back
void on_complete ();

void radiation_driver_lw (int ncol, int nlev,
const real3d& gas_vmr,
const real2d& pmid, const real2d& pint, const real2d& tmid, const real2d& tint,
const real3d& cld_tau_gpt, const real3d& aer_tau_bnd,
FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
const real2d& qrl, const real2d& qrlc);

void radiation_driver_sw (int ncol,
const real3d& gas_vmr, const real2d& pmid, const real2d& pint, const real2d& tmid,
const real2d& albedo_dir, const real2d& albedo_dif, const real1d& coszrs,
const real3d& cld_tau_gpt, const real3d& cld_ssa_gpt, const real3d& cld_asm_gpt,
const real3d& aer_tau_bnd, const real3d& aer_ssa_bnd, const real3d& aer_asm_bnd,
FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
void on_complete();

void radiation_driver_lw(int ncol, int nlev,
const real3d& gas_vmr,
const real2d& pmid, const real2d& pint, const real2d& tmid, const real2d& tint,
const real3d& cld_tau_gpt, const real3d& aer_tau_bnd,
FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
const real2d& qrl, const real2d& qrlc);

void radiation_driver_sw(int ncol,
const real3d& gas_vmr, const real2d& pmid, const real2d& pint, const real2d& tmid,
const real2d& albedo_dir, const real2d& albedo_dif, const real1d& coszrs,
const real3d& cld_tau_gpt, const real3d& cld_ssa_gpt, const real3d& cld_asm_gpt,
const real3d& aer_tau_bnd, const real3d& aer_ssa_bnd, const real3d& aer_asm_bnd,
FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
const real2d& qrs, const real2d& qrsc);

void set_daynight_indices (const real1d& coszrs,
const int1d& day_indices,
const int1d& night_indices);
void set_daynight_indices(const real1d& coszrs,
const int1d& day_indices,
const int1d& night_indices);

void get_gas_vmr (const std::vector<std::string>& gas_names,
const real3d& gas_vmr);
void get_gas_vmr(const std::vector<std::string>& gas_names,
const real3d& gas_vmr);

void calculate_heating_rate (const real2d& flux_up,
const real2d& flux_dn,
const real2d& pint,
const real2d& heating_rate);
void calculate_heating_rate(const real2d& flux_up,
const real2d& flux_dn,
const real2d& pint,
const real2d& heating_rate);
private:
// geometry
amrex::Geometry m_geom;
Expand Down
Loading

0 comments on commit 02493b8

Please sign in to comment.