Radiation Class Reference

#include <ERF_Radiation.H>

Collaboration diagram for Radiation:

Public Member Functions
	Radiation ()
	~Radiation ()=default
void	initialize (const amrex::MultiFab &cons_in, amrex::MultiFab *lsm_fluxes, amrex::MultiFab *lsm_zenith, amrex::MultiFab *qheating_rates, amrex::MultiFab *lat, amrex::MultiFab *lon, amrex::Vector< amrex::MultiFab * > qmoist, const amrex::BoxArray &grids, const amrex::Geometry &geom, const amrex::Real &dt_advance, const bool &do_sw_rad, const bool &do_lw_rad, const bool &do_aero_rad, const bool &do_snow_opt, const bool &is_cmip6_volcano)
void	run ()
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	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	export_surface_fluxes (FluxesByband &fluxes, std::string band)
void	yakl_to_mf (const real2d &data, amrex::MultiFab &mf)
void	expand_yakl1d_to_mf (const real1d &data, amrex::MultiFab &mf)
void	writePlotfile (const std::string &plot_prefix, const amrex::Real time, const int level_step)

Private Attributes
amrex::Geometry	m_geom
amrex::BoxArray	m_box
amrex::MultiFab *	qrad_src = nullptr
amrex::MultiFab *	m_lat = nullptr
amrex::MultiFab *	m_lon = nullptr
amrex::MultiFab *	m_lsm_fluxes = nullptr
amrex::MultiFab *	m_lsm_zenith = nullptr
std::string	moisture_type = "None"
bool	has_qmoist
amrex::Vector< int >	rank_offsets
amrex::Real	uniform_angle = 78.463
int	nlev
int	zlo
int	zhi
int	ncol
int	nlwgpts
int	nswgpts
int	nlwbands
int	nswbands
bool	do_short_wave_rad
bool	do_long_wave_rad
bool	do_snow_optics
bool	do_aerosol_rad = true
Rrtmgp	radiation
Optics	optics
AerRadProps	aer_rad
real1d	net_flux
bool	is_cmip6_volc
real	dt
real1d	fsns
real1d	fsnt
real1d	flns
real1d	flnt
real1d	fsds
const std::vector< std::string >	active_gases
bool	spectralflux = false
bool	use_rad_dt_cosz = false
real	fixed_total_solar_irradiance = -1.
bool	rrtmgp_enable_temperature_warnings = true
bool	dohirs = false
int	ihirsfq = 1
real	dt_avg = 0.0
std::string	rrtmgp_data_path
std::string	rrtmgp_coefficients_file_sw
std::string	rrtmgp_coefficients_file_lw
std::string	rrtmgp_coefficients_file_name_sw = "rrtmgp_coefficients_sw_20181204.nc"
std::string	rrtmgp_coefficients_file_name_lw = "rrtmgp_coefficients_lw_20181204.nc"
real1d	sw_band_midpoints
real1d	lw_band_midpoints
int	ngas
int	naer
std::vector< std::string >	gasnames
std::vector< std::string >	aernames
int1d	rrtmg_to_rrtmgp
real2d	qrs
real2d	qrl
real2d	zi
real2d	clear_rh
real2d	qrsc
real2d	qrlc
real2d	qt
real2d	qi
real2d	qc
real2d	qn
real2d	tmid
real2d	pmid
real2d	pdel
real2d	pint
real2d	tint
real2d	albedo_dir
real2d	albedo_dif
real1d	coszrs
real2d	cld
real2d	cldfsnow
real2d	iclwp
real2d	iciwp
real2d	icswp
real2d	dei
real2d	des
real2d	lambdac
real2d	mu
real2d	rei
real2d	rel
real3d	cld_tau_gpt_sw
real3d	cld_ssa_gpt_sw
real3d	cld_asm_gpt_sw
real3d	cld_tau_bnd_sw
real3d	cld_ssa_bnd_sw
real3d	cld_asm_bnd_sw
real3d	aer_tau_bnd_sw
real3d	aer_ssa_bnd_sw
real3d	aer_asm_bnd_sw
real3d	cld_tau_bnd_lw
real3d	aer_tau_bnd_lw
real3d	cld_tau_gpt_lw
real3d	liq_tau_bnd_sw
real3d	ice_tau_bnd_sw
real3d	snw_tau_bnd_sw
real3d	liq_tau_bnd_lw
real3d	ice_tau_bnd_lw
real3d	snw_tau_bnd_lw
real3d	gas_vmr
int1d	gpoint_bands_sw
int1d	gpoint_bands_lw
FluxesByband	sw_fluxes_allsky
FluxesByband	sw_fluxes_clrsky
real1d	cld_tau_bnd_sw_1d
real1d	cld_ssa_bnd_sw_1d
real1d	cld_asm_bnd_sw_1d
real1d	cld_tau_bnd_sw_o_1d
real1d	cld_ssa_bnd_sw_o_1d
real1d	cld_asm_bnd_sw_o_1d
real1d	aer_tau_bnd_sw_1d
real1d	aer_ssa_bnd_sw_1d
real1d	aer_asm_bnd_sw_1d
real1d	aer_tau_bnd_sw_o_1d
real1d	aer_ssa_bnd_sw_o_1d
real1d	aer_asm_bnd_sw_o_1d
FluxesByband	lw_fluxes_allsky
FluxesByband	lw_fluxes_clrsky

Static Private Attributes
static constexpr amrex::Real	eccen = 0.0
static constexpr amrex::Real	obliqr = 23.0 * PI / 180.0
static constexpr amrex::Real	mvelpp = 103.0 * PI / 180.0 + PI
static constexpr amrex::Real	lambm0 = -3.2503635878519378e-2

Constructor & Destructor Documentation

Radiation()

Radiation::Radiation ( )

inline

                 {
        // First, make sure yakl has been initialized
        if (!yakl::isInitialized()) yakl::init();
    }

~Radiation()

Radiation::~Radiation ( )

default

Member Function Documentation

calculate_heating_rate()

void Radiation::calculate_heating_rate	(	const real2d &	flux_up,
		const real2d &	flux_dn,
		const real2d &	pint,
		const real2d &	heating_rate
	)

{
    // NOTE: The pressure is in [pa] for RRTMGP to use.
    //       The fluxes are in [W/m^2] and gravity is [m/s^2].
    //       The heating rate is {dF/dP * g / Cp} with units [K/s]
    real1d heatfac("heatfac",1);
    yakl::memset(heatfac, 1.0/Cp_d);
    parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int 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+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_yakl1d_to_mf()

void Radiation::expand_yakl1d_to_mf	(	const real1d &	data,
		amrex::MultiFab &	mf
	)

{
    // copies the 1D yakl data to a 3D MF
    AMREX_ASSERT(data.get_dimensions()(1) == ncol);
    mf = amrex::MultiFab(m_box, qrad_src->DistributionMap(), 1, 0);
    if (!data.initialized())
    {
        // yakl array hasn't been created yet, so create an empty MF
        mf.setVal(0.0);
        return;
    }
 
    for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
        auto mf_arr = mf.array(mfi);
        const auto& box3d = mfi.tilebox();
        const int nx = box3d.length(0);
        const int offset = rank_offsets[mfi.LocalIndex()];
        amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            // map [i,j,k] 0-based to [icol, ilev] 1-based
            const int icol = (j-box3d.smallEnd(1))*nx + (i-box3d.smallEnd(0)) + 1 + offset;
            AMREX_ASSERT(icol <= static_cast<int>(data.get_dimensions()(1)));
            mf_arr(i, j, k) = data(icol);
        });
    }
}

Here is the call graph for this function:

export_surface_fluxes()

void Radiation::export_surface_fluxes	(	FluxesByband &	fluxes,
		std::string	band
	)

{
    // No work to be done if we don't have valid pointers
    if (!m_lsm_fluxes) return;
 
    if (band == "shortwave") {
        real3d flux_dn_diffuse("flux_dn_diffuse", ncol, nlev+1, nswbands);
 
        // Calculate diffuse flux from total and direct
        // This only occurs at the bottom level (k index)
        parallel_for (SimpleBounds<3>(ncol, 1, nswbands), YAKL_LAMBDA (int icol, int ilev, int ibnd)
        {
            flux_dn_diffuse(icol,ilev,ibnd) = fluxes.bnd_flux_dn(icol,ilev,ibnd)
                                            - fluxes.bnd_flux_dn_dir(icol,ilev,ibnd);
        });
 
        // Populate the LSM data structure (this is a 2D MF)
        for (MFIter mfi(*(m_lsm_fluxes)); mfi.isValid(); ++mfi) {
            auto lsm_array = m_lsm_fluxes->array(mfi);
            const auto& box3d = mfi.tilebox();
            auto nx = box3d.length(0);
            const int offset = rank_offsets[mfi.LocalIndex()];
            amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                // Map (col,lev) to (i,j,k)
                auto icol = (j-box3d.smallEnd(1))*nx + (i-box3d.smallEnd(0)) + 1 + offset;
                auto ilev = k+1;
 
                // Direct fluxes
                Real sum1(0.0), sum2(0.0);
                for (int ibnd(1); ibnd<=9; ++ibnd) {
                    sum1 += fluxes.bnd_flux_dn_dir(icol,ilev,ibnd);
                }
                for (int ibnd(11); ibnd<=14; ++ibnd) {
                    sum2 += fluxes.bnd_flux_dn_dir(icol,ilev,ibnd);
                }
                sum1 += 0.5 * fluxes.bnd_flux_dn_dir(icol,ilev,10);
                sum2 += 0.5 * fluxes.bnd_flux_dn_dir(icol,ilev,10);
                lsm_array(i,j,k,0) = sum1;
                lsm_array(i,j,k,1) = sum2;
 
                // Diffuse fluxes
                sum1=0.0; sum2=0.0;
                for (int ibnd(1); ibnd<=9; ++ibnd) {
                    sum1 += flux_dn_diffuse(icol,ilev,ibnd);
                }
                for (int ibnd(11); ibnd<=14; ++ibnd) {
                    sum2 += flux_dn_diffuse(icol,ilev,ibnd);
                }
                sum1 += 0.5 * flux_dn_diffuse(icol,ilev,10);
                sum2 += 0.5 * flux_dn_diffuse(icol,ilev,10);
                lsm_array(i,j,k,2) = sum1;
                lsm_array(i,j,k,3) = sum2;
 
                // Net fluxes
                lsm_array(i,j,k,4) = fluxes.flux_net(icol,ilev);
            });
        }
    } else if (band == "longwave") {
        // Populate the LSM data structure (this is a 2D MF)
        for (MFIter mfi(*(m_lsm_fluxes)); mfi.isValid(); ++mfi) {
            auto lsm_array = m_lsm_fluxes->array(mfi);
            const auto& box3d = mfi.tilebox();
            auto nx = box3d.length(0);
            const int offset = rank_offsets[mfi.LocalIndex()];
            amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
            {
                // Map (col,lev) to (i,j,k)
                auto icol = (j-box3d.smallEnd(1))*nx + (i-box3d.smallEnd(0)) + 1 + offset;
                auto ilev = k+1;
 
                // Net fluxes
                lsm_array(i,j,k,5) = fluxes.flux_dn(icol,ilev);
            });
        }
    } else {
         amrex::Abort("Unknown radiation band type!");
    }
}

Here is the call graph for this function:

get_gas_vmr()

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

{
    // Mass mixing ratio
    real2d mmr("mmr", ncol, nlev);
 
    // Gases and molecular weights. Note that we do NOT have CFCs yet (I think
    // this is coming soon in RRTMGP). RRTMGP also allows for absorption due to
    // CO and N2, which RRTMG did not have.
    const std::vector<std::string> gas_species = {"H2O", "CO2", "O3", "N2O",
                                                  "CO" , "CH4", "O2", "N2"};
    const std::vector<real> mol_weight_gas = {18.01528, 44.00950, 47.9982, 44.0128,
                                              28.01010, 16.04246, 31.9980, 28.0134}; // g/mol
    // Molar weight of air
    //const real mol_weight_air = 28.97; // g/mol
    // Defaults for gases that are not available (TODO: is this still accurate?)
    const real co_vol_mix_ratio = 1.0e-7;
    const real n2_vol_mix_ratio = 0.7906;
 
    // initialize
    yakl::memset(gas_vmr, 0.);
 
    // For each gas species needed for RRTMGP, read the mass mixing ratio from the
    // CAM rad_constituents interface, convert to volume mixing ratios, and
    // subset for daytime-only indices if needed.
    for (auto igas = 0; igas < gas_names.size(); ++igas) {
 
        std::cout << "gas_name: " << gas_names[igas] << "; igas: " << igas << std::endl;
 
        if (gas_names[igas] == "CO"){
            // CO not available, use default
            parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
            {
                gas_vmr(igas+1,icol,ilev) = co_vol_mix_ratio;
            });
        } else if (gas_names[igas] == "N2") {
            // N2 not available, use default
            parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
            {
                gas_vmr(igas+1,icol,ilev) = n2_vol_mix_ratio;
            });
        } else if (gas_names[igas] == "H2O") {
            // Water vapor is represented as specific humidity in CAM, so we
            // need to handle water a little differently
            //        rad_cnst_get_gas(icall, gas_species[igas], mmr);
 
            // Convert to volume mixing ratio by multiplying by the ratio of
            // molecular weight of dry air to molecular weight of gas. Note that
            // first specific humidity (held in the mass_mix_ratio array read
            // from rad_constituents) is converted to an actual mass mixing ratio.
            parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
            {
                gas_vmr(igas+1,icol,ilev) = qt(icol,ilev); //mmr(icol,ilev) / (
                //                  1. - mmr(icol,ilev))*mol_weight_air / mol_weight_gas[igas];
            });
        } else if (gas_names[igas] == "CO2") {
            parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
            {
                gas_vmr(igas+1,icol,ilev) = 3.8868676125307193E-4;
            });
        } else if (gas_names[igas] == "O3") {
            parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
            {
                gas_vmr(igas+1,icol,ilev) = 1.8868676125307193E-7;
            });
        } else if (gas_names[igas] == "N2O") {
            parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
            {
                gas_vmr(igas+1,icol,ilev) = 3.8868676125307193E-7;
            });
        } else if (gas_names[igas] == "CH4") {
            parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
            {
                gas_vmr(igas+1,icol,ilev) = 1.8868676125307193E-6;
            });
        } else if (gas_names[igas] == "O2") {
            parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
            {
                gas_vmr(igas+1,icol,ilev) = 0.2095;
            });
        } else {
            // Get mass mixing ratio from the rad_constituents interface
            //        rad_cnst_get_gas(icall, gas_species[igas], mmr);
 
            // Convert to volume mixing ratio by multiplying by the ratio of
            // molecular weight of dry air to molecular weight of gas
            parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
            {
                gas_vmr(igas+1,icol,ilev) = 1.0e-6; //mmr(icol,ilev)
                //                                     * mol_weight_air / mol_weight_gas[igas];
            });
        }
    } // igas
}

initialize()

void Radiation::initialize	(	const amrex::MultiFab &	cons_in,
		amrex::MultiFab *	lsm_fluxes,
		amrex::MultiFab *	lsm_zenith,
		amrex::MultiFab *	qheating_rates,
		amrex::MultiFab *	lat,
		amrex::MultiFab *	lon,
		amrex::Vector< amrex::MultiFab * >	qmoist,
		const amrex::BoxArray &	grids,
		const amrex::Geometry &	geom,
		const amrex::Real &	dt_advance,
		const bool &	do_sw_rad,
		const bool &	do_lw_rad,
		const bool &	do_aero_rad,
		const bool &	do_snow_opt,
		const bool &	is_cmip6_volcano
	)

{
    m_geom = geom;
    m_box = grids;
 
    qrad_src = qheating_rates;
 
    auto dz   = m_geom.CellSize(2);
    auto lowz = m_geom.ProbLo(2);
 
    dt = dt_advance;
 
    do_short_wave_rad = do_sw_rad;
    do_long_wave_rad  = do_lw_rad;
    do_aerosol_rad    = do_aero_rad;
    do_snow_optics    = do_snow_opt;
    is_cmip6_volc     = is_cmip6_volcano;
 
    m_lat = lat;
    m_lon = lon;
 
    m_lsm_fluxes = lsm_fluxes;
    m_lsm_zenith = lsm_zenith;
 
    rrtmgp_data_path = getRadiationDataDir() + "/";
    rrtmgp_coefficients_file_sw = rrtmgp_data_path + rrtmgp_coefficients_file_name_sw;
    rrtmgp_coefficients_file_lw = rrtmgp_data_path + rrtmgp_coefficients_file_name_lw;
 
    ParmParse pp("erf");
    pp.query("fixed_total_solar_irradiance", fixed_total_solar_irradiance);
    pp.query("radiation_uniform_angle"     , uniform_angle);
    pp.query("moisture_model", moisture_type); // TODO: get from SolverChoice?
    has_qmoist = (moisture_type != "None");
 
    nlev = geom.Domain().length(2);
    ncol = 0;
    rank_offsets.resize(cons_in.local_size());
    for (MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) {
        const auto& box3d = mfi.tilebox();
        int nx = box3d.length(0);
        int ny = box3d.length(1);
        rank_offsets[mfi.LocalIndex()] = ncol;
        ncol += nx * ny;
    }
 
    ngas = active_gases.size();
 
    // initialize cloud, aerosol, and radiation
    radiation.initialize(ngas, active_gases,
                         rrtmgp_coefficients_file_sw.c_str(),
                         rrtmgp_coefficients_file_lw.c_str());
 
    // initialize the radiation data
    nswbands = radiation.get_nband_sw();
    nswgpts  = radiation.get_ngpt_sw();
    nlwbands = radiation.get_nband_lw();
    nlwgpts  = radiation.get_ngpt_lw();
 
    rrtmg_to_rrtmgp = int1d("rrtmg_to_rrtmgp",14);
    parallel_for(14, YAKL_LAMBDA (int i)
    {
        if (i == 1) {
            rrtmg_to_rrtmgp(i) = 13;
        } else {
            rrtmg_to_rrtmgp(i) = i - 1;
        }
    });
 
    tmid = real2d("tmid", ncol, nlev);
    pmid = real2d("pmid", ncol, nlev);
    pdel = real2d("pdel", ncol, nlev);
 
    pint = real2d("pint", ncol, nlev+1);
    tint = real2d("tint", ncol, nlev+1);
 
    qt   = real2d("qt", ncol, nlev);
    qc   = real2d("qc", ncol, nlev);
    qi   = real2d("qi", ncol, nlev);
    qn   = real2d("qn", ncol, nlev);
    zi   = real2d("zi", ncol, nlev);
 
    // Get the temperature, density, theta, qt and qp from input
    for (MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) {
        const auto& box3d = mfi.tilebox();
        auto nx = box3d.length(0);
 
        auto states_array = cons_in.array(mfi);
        auto qt_array = (has_qmoist) ? qmoist[0]->array(mfi) : Array4<Real> {};
        auto qv_array = (has_qmoist) ? qmoist[1]->array(mfi) : Array4<Real> {};
        auto qc_array = (has_qmoist) ? qmoist[2]->array(mfi) : Array4<Real> {};
        auto qi_array = (has_qmoist && qmoist.size()>=8) ? qmoist[3]->array(mfi) : Array4<Real> {};
        const int offset = rank_offsets[mfi.LocalIndex()];
 
        // Get pressure, theta, temperature, density, and qt, qp
        ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            auto icol = (j-box3d.smallEnd(1))*nx + (i-box3d.smallEnd(0)) + 1 + offset;
            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),qv);
            // NOTE: RRTMGP code expects pressure in pa
            pmid(icol,ilev) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp),qv);
        });
    }
 
    parallel_for(SimpleBounds<2>(ncol, nlev+1), YAKL_LAMBDA (int icol, int ilev)
    {
        if (ilev == 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) = -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);
        }
    });
 
    parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
    {
        zi(icol, ilev)  = lowz + (ilev+0.5)*dz;
        pdel(icol,ilev) = pint(icol,ilev+1) - pint(icol,ilev);
    });
 
    albedo_dir = real2d("albedo_dir", nswbands, ncol);
    albedo_dif = real2d("albedo_dif", nswbands, ncol);
 
    qrs = real2d("qrs", ncol, nlev);   // shortwave radiative heating rate
    qrl = real2d("qrl", ncol, nlev);   // longwave  radiative heating rate
 
    // Clear-sky heating rates are not on the physics buffer, and we have no
    // reason to put them there, so declare these are regular arrays here
    qrsc = real2d("qrsc", ncol, nlev);
    qrlc = real2d("qrlc", ncol, nlev);
 
    int nmodes = 3;
    int nrh = 1;
    int top_lev = 1;
    naer = 4;
    std::vector<std::string> aero_names {"H2O", "N2", "O2", "O3"};
    auto geom_radius = real2d("geom_radius", ncol, nlev);
    yakl::memset(geom_radius, 0.1);
 
    optics.initialize(ngas, nmodes, naer, nswbands, nlwbands,
                      ncol, nlev, nrh, top_lev, aero_names, zi,
                      pmid, pdel, tmid, qt, geom_radius);
 
    amrex::Print() << "LW coefficients file: " << rrtmgp_coefficients_file_lw
                   << "\nSW coefficients file: " << rrtmgp_coefficients_file_sw
                   << "\nFrequency (timesteps) of Shortwave Radiation calc: " << dt
                   << "\nFrequency (timesteps) of Longwave Radiation calc:  " << dt
                   << "\nDo aerosol radiative calculations: " << do_aerosol_rad << std::endl;
 
}

Here is the call graph for this function:

on_complete()

void Radiation::on_complete

(

)

{ }

radiation_driver_lw()

void Radiation::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
	)

{
    real3d cld_tau_gpt_rad("cld_tau_gpt_rad", ncol, nlev+1, nlwgpts);
    real3d aer_tau_bnd_rad("aer_tau_bnd_rad", ncol, nlev+1, nlwgpts);
 
    // Surface emissivity needed for longwave
    real2d surface_emissivity("surface_emissivity", nlwbands, ncol);
 
    // Temporary heating rates on radiation vertical grid
    real3d gas_vmr_rad("gas_vmr_rad", ngas, ncol, nlev);
 
    // Set surface emissivity to 1 here. There is a note in the RRTMG
    // implementation that this is treated in the land model, but the old
    // RRTMG implementation also sets this to 1. This probably does not make
    // a lot of difference either way, but if a more intelligent value
    // exists or is assumed in the model we should use it here as well.
    // TODO: set this more intelligently?
    yakl::memset(surface_emissivity, 1.0);
 
    // Add an empty level above model top
    yakl::memset(cld_tau_gpt_rad, 0.);
    yakl::memset(aer_tau_bnd_rad, 0.);
 
    parallel_for(SimpleBounds<3>(ncol, nlev, nlwgpts), YAKL_LAMBDA (int icol, int ilev, int igpt)
    {
        cld_tau_gpt_rad(icol,ilev,igpt) = cld_tau_gpt(icol,ilev,igpt);
        aer_tau_bnd_rad(icol,ilev,igpt) = aer_tau_bnd(icol,ilev,igpt);
    });
 
    parallel_for(SimpleBounds<3>(ncol, nlev, ngas), YAKL_LAMBDA (int icol, int ilev, int igas)
    {
        gas_vmr_rad(igas,icol,ilev) = gas_vmr(igas,icol,ilev);
    });
 
    // Do longwave radiative transfer calculations
    radiation.run_longwave_rrtmgp(ngas, ncol, nlev,
                                  gas_vmr_rad, pmid, tmid, pint, tint,
                                  surface_emissivity, cld_tau_gpt_rad, aer_tau_bnd_rad,
                                  fluxes_allsky.flux_up    , fluxes_allsky.flux_dn    , fluxes_allsky.flux_net    ,
                                  fluxes_allsky.bnd_flux_up, fluxes_allsky.bnd_flux_dn, fluxes_allsky.bnd_flux_net,
                                  fluxes_clrsky.flux_up    , fluxes_clrsky.flux_dn    , fluxes_clrsky.flux_net    ,
                                  fluxes_clrsky.bnd_flux_up, fluxes_clrsky.bnd_flux_dn, fluxes_clrsky.bnd_flux_net);
 
    // Calculate heating rates
    calculate_heating_rate(fluxes_allsky.flux_up,
                           fluxes_allsky.flux_dn,
                           pint, qrl);
 
    calculate_heating_rate(fluxes_allsky.flux_up,
                           fluxes_allsky.flux_dn,
                           pint, qrlc);
}

radiation_driver_sw()

void Radiation::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
	)

{
    // Incoming solar radiation, scaled for solar zenith angle
    // and earth-sun distance
    real2d solar_irradiance_by_gpt("solar_irradiance_by_gpt",ncol,nswgpts);
 
    // Gathered indices of day and night columns
    // chunk_column_index = day_indices(daylight_column_index)
    int1d day_indices("day_indices",ncol), night_indices("night_indices", ncol);   // Indices of daylight coumns
 
    real1d coszrs_day("coszrs_day", ncol);
    real2d albedo_dir_day("albedo_dir_day", nswbands, ncol), albedo_dif_day("albedo_dif_day", nswbands, ncol);
    real2d pmid_day("pmid_day", ncol, nlev);
    real2d tmid_day("tmid_day", ncol, nlev);
    real2d pint_day("pint_day", ncol, nlev+1);
 
    real3d gas_vmr_day("gas_vmr_day", ngas, ncol, nlev);
 
    real3d cld_tau_gpt_day("cld_tau_gpt_day", ncol, nlev, nswgpts);
    real3d cld_ssa_gpt_day("cld_ssa_gpt_day", ncol, nlev, nswgpts);
    real3d cld_asm_gpt_day("cld_asm_gpt_day", ncol, nlev, nswgpts);
    real3d aer_tau_bnd_day("aer_tau_bnd_day", ncol, nlev, nswbands);
    real3d aer_ssa_bnd_day("aer_ssa_bnd_day", ncol, nlev, nswbands);
    real3d aer_asm_bnd_day("aer_asm_bnd_day", ncol, nlev, nswbands);
 
    real3d cld_tau_gpt_rad("cld_tau_gpt_rad", ncol, nlev+1, nswgpts);
    real3d cld_ssa_gpt_rad("cld_ssa_gpt_rad", ncol, nlev+1, nswgpts);
    real3d cld_asm_gpt_rad("cld_asm_gpt_rad", ncol, nlev+1, nswgpts);
    real3d aer_tau_bnd_rad("aer_tau_bnd_rad", ncol, nlev+1, nswgpts);
    real3d aer_ssa_bnd_rad("aer_ssa_bnd_rad", ncol, nlev+1, nswgpts);
    real3d aer_asm_bnd_rad("aer_asm_bnd_rad", ncol, nlev+1, nswgpts);
 
    // Scaling factor for total sky irradiance; used to account for orbital
    // eccentricity, and could be used to scale total sky irradiance for different
    // climates as well (i.e., paleoclimate simulations)
    real tsi_scaling;
    real solar_declination;
 
    if (fixed_total_solar_irradiance<0) {
        // TODO: Integrate calendar day computation
        int calday = 1;
        // Get orbital eccentricity factor to scale total sky irradiance
        shr_orb_decl(calday, eccen, mvelpp, lambm0, obliqr, solar_declination, tsi_scaling);
    } else {
        // For fixed TSI we divide by the default solar constant of 1360.9
        // At some point we will want to replace this with a method that
        // retrieves the solar constant
        tsi_scaling = fixed_total_solar_irradiance / 1360.9;
    }
 
    // Gather night/day column indices for subsetting SW inputs; we only want to
    // do the shortwave radiative transfer during the daytime to save
    // computational cost (and because RRTMGP will fail for cosine solar zenith
    // angles less than or equal to zero)
    set_daynight_indices(coszrs, day_indices, night_indices);
    int1d nday("nday",1);
    int1d nnight("nnight", 1);
    yakl::memset(nday, 0);
    yakl::memset(nnight, 0);
    parallel_for(SimpleBounds<1>(ncol), YAKL_LAMBDA (int icol)
    {
        if (day_indices(icol) > 0) nday(1)++;
        if (night_indices(icol) > 0) nnight(1)++;
    });
 
    AMREX_ASSERT(nday(1) + nnight(1) == ncol);
 
    intHost1d num_day("num_day",1);
    intHost1d num_night("num_night",1);
    nday.deep_copy_to(num_day);
    nnight.deep_copy_to(num_night);
 
    // If no daytime columns in this chunk, then we return zeros
    if (num_day(1) == 0) {
        //    reset_fluxes(fluxes_allsky)
        //    reset_fluxes(fluxes_clrsky)
        yakl::memset(qrs, 0.);
        yakl::memset(qrsc, 0.);
        return;
    }
 
    // Compress to daytime-only arrays
    parallel_for(SimpleBounds<2>(num_day(1), nlev), YAKL_LAMBDA (int iday, int ilev)
    {
        // 2D arrays
        auto icol = day_indices(iday);
        tmid_day(iday,ilev) = tmid(icol,ilev);
        pmid_day(iday,ilev) = pmid(icol,ilev);
        pint_day(iday,ilev) = pint(icol,ilev);
    });
    parallel_for(SimpleBounds<1>(num_day(1)), YAKL_LAMBDA (int iday)
    {
        // copy extra level for pmid
        auto icol = day_indices(iday);
        pint_day(iday,nlev+1) = pint(icol,nlev+1);
 
        coszrs_day(iday) = coszrs(icol);
        AMREX_ASSERT(coszrs_day(iday) > 0.0);
    });
    parallel_for(SimpleBounds<3>(num_day(1), nlev, nswgpts), YAKL_LAMBDA (int iday, int ilev, int igpt)
    {
        auto icol = day_indices(iday);
        cld_tau_gpt_day(iday,ilev,igpt) = cld_tau_gpt(icol,ilev,igpt);
        cld_ssa_gpt_day(iday,ilev,igpt) = cld_ssa_gpt(icol,ilev,igpt);
        cld_asm_gpt_day(iday,ilev,igpt) = cld_asm_gpt(icol,ilev,igpt);
    });
    parallel_for(SimpleBounds<2>(num_day(1), nswbands), YAKL_LAMBDA (int iday, int ibnd)
    {
        // albedo dims: [nswbands, ncol]
        auto icol = day_indices(iday);
        albedo_dir_day(ibnd,iday) = albedo_dir(ibnd,icol);
        albedo_dif_day(ibnd,iday) = albedo_dif(ibnd,icol);
    });
 
    parallel_for(SimpleBounds<3>(num_day(1), nlev, nswbands), YAKL_LAMBDA (int iday, int ilev, int ibnd)
    {
        auto icol = day_indices(iday);
        aer_tau_bnd_day(iday,ilev,ibnd) = aer_tau_bnd(icol,ilev,ibnd);
        aer_ssa_bnd_day(iday,ilev,ibnd) = aer_ssa_bnd(icol,ilev,ibnd);
        aer_asm_bnd_day(iday,ilev,ibnd) = aer_asm_bnd(icol,ilev,ibnd);
    });
 
    // Allocate shortwave fluxes (allsky and clearsky)
    // NOTE: fluxes defined at interfaces, so initialize to have vertical
    // dimension nlev_rad+1, while we initialized the RRTMGP input variables to
    // have vertical dimension nlev_rad (defined at midpoints).
    FluxesByband fluxes_clrsky_day, fluxes_allsky_day;
    internal::initial_fluxes(num_day(1), nlev+1, nswbands, fluxes_allsky_day);
    internal::initial_fluxes(num_day(1), nlev+1, nswbands, fluxes_clrsky_day);
 
    // Add an empty level above model top
    // TODO: combine with day compression above
    yakl::memset(cld_tau_gpt_rad, 0.);
    yakl::memset(cld_ssa_gpt_rad, 0.);
    yakl::memset(cld_asm_gpt_rad, 0.);
 
    yakl::memset(aer_tau_bnd_rad, 0.);
    yakl::memset(aer_ssa_bnd_rad, 0.);
    yakl::memset(aer_asm_bnd_rad, 0.);
 
    parallel_for(SimpleBounds<3>(num_day(1), nlev, nswgpts), YAKL_LAMBDA (int iday, int ilev, int igpt)
    {
        cld_tau_gpt_rad(iday,ilev,igpt) = cld_tau_gpt_day(iday,ilev,igpt);
        cld_ssa_gpt_rad(iday,ilev,igpt) = cld_ssa_gpt_day(iday,ilev,igpt);
        cld_asm_gpt_rad(iday,ilev,igpt) = cld_asm_gpt_day(iday,ilev,igpt);
    });
 
    parallel_for(SimpleBounds<3>(num_day(1), nlev, ngas), YAKL_LAMBDA (int iday, int ilev, int igas)
    {
        auto icol = day_indices(iday);
        gas_vmr_day(igas,iday,ilev) = gas_vmr(igas,icol,ilev);
    });
 
    parallel_for(SimpleBounds<3>(num_day(1), nlev, nswbands), YAKL_LAMBDA (int iday, int ilev, int ibnd)
    {
        aer_tau_bnd_rad(iday,ilev,ibnd) = aer_tau_bnd_day(iday,ilev,ibnd);
        aer_ssa_bnd_rad(iday,ilev,ibnd) = aer_ssa_bnd_day(iday,ilev,ibnd);
        aer_asm_bnd_rad(iday,ilev,ibnd) = aer_asm_bnd_day(iday,ilev,ibnd);
    });
 
    // Do shortwave radiative transfer calculations
    radiation.run_shortwave_rrtmgp(ngas, num_day(1), nlev, gas_vmr_day, pmid_day,
                                   tmid_day, pint_day, coszrs_day, albedo_dir_day, albedo_dif_day,
                                   cld_tau_gpt_rad, cld_ssa_gpt_rad, cld_asm_gpt_rad, aer_tau_bnd_rad, aer_ssa_bnd_rad, aer_asm_bnd_rad,
                                   fluxes_allsky_day.flux_up    , fluxes_allsky_day.flux_dn    , fluxes_allsky_day.flux_net    , fluxes_allsky_day.flux_dn_dir    ,
                                   fluxes_allsky_day.bnd_flux_up, fluxes_allsky_day.bnd_flux_dn, fluxes_allsky_day.bnd_flux_net, fluxes_allsky_day.bnd_flux_dn_dir,
                                   fluxes_clrsky_day.flux_up    , fluxes_clrsky_day.flux_dn    , fluxes_clrsky_day.flux_net    , fluxes_clrsky_day.flux_dn_dir    ,
                                   fluxes_clrsky_day.bnd_flux_up, fluxes_clrsky_day.bnd_flux_dn, fluxes_clrsky_day.bnd_flux_net, fluxes_clrsky_day.bnd_flux_dn_dir,
                                   tsi_scaling);
 
    // Expand fluxes from daytime-only arrays to full chunk arrays
    internal::expand_day_fluxes(fluxes_allsky_day, fluxes_allsky, day_indices);
    internal::expand_day_fluxes(fluxes_clrsky_day, fluxes_clrsky, day_indices);
 
    // Calculate heating rates
    calculate_heating_rate(fluxes_allsky.flux_up,
                           fluxes_allsky.flux_dn,
                           pint, qrs);
 
    calculate_heating_rate(fluxes_clrsky.flux_up,
                           fluxes_allsky.flux_dn,
                           pint, qrsc);
}

Here is the call graph for this function:

run()

void Radiation::run

(

)

{
    // Cosine solar zenith angle for all columns in chunk
    coszrs = real1d("coszrs", ncol);
 
    // Pointers to fields on the physics buffer
    cld = real2d("cld", ncol, nlev);
    cldfsnow = real2d("cldfsnow", ncol, nlev);
    iclwp = real2d("iclwp", ncol, nlev);
    iciwp = real2d("iciwp", ncol, nlev);
    icswp = real2d("icswp", ncol, nlev);
    dei = real2d("dei", ncol, nlev);
    des = real2d("des", ncol, nlev);
    lambdac = real2d("lambdac", ncol, nlev);
    mu = real2d("mu", ncol, nlev);
    rei = real2d("rei", ncol, nlev);
    rel = real2d("rel", ncol, nlev);
 
    // Cloud, snow, and aerosol optical properties
    cld_tau_gpt_sw = real3d("cld_tau_gpt_sw", ncol, nlev, nswgpts);
    cld_ssa_gpt_sw = real3d("cld_ssa_gpt_sw", ncol, nlev, nswgpts);
    cld_asm_gpt_sw = real3d("cld_asm_gpt_sw", ncol, nlev, nswgpts);
 
    cld_tau_bnd_sw = real3d("cld_tau_bnd_sw", ncol, nlev, nswbands);
    cld_ssa_bnd_sw = real3d("cld_ssa_bnd_sw", ncol, nlev, nswbands);
    cld_asm_bnd_sw = real3d("cld_asm_bnd_sw", ncol, nlev, nswbands);
 
    aer_tau_bnd_sw = real3d("aer_tau_bnd_sw", ncol, nlev, nswbands);
    aer_ssa_bnd_sw = real3d("aer_ssa_bnd_sw", ncol, nlev, nswbands);
    aer_asm_bnd_sw = real3d("aer_asm_bnd_sw", ncol, nlev, nswbands);
 
    cld_tau_bnd_lw = real3d("cld_tau_bnd_lw", ncol, nlev, nlwbands);
    aer_tau_bnd_lw = real3d("aer_tau_bnd_lw", ncol, nlev, nlwbands);
 
    cld_tau_gpt_lw = real3d("cld_tau_gpt_lw", ncol, nlev, nlwgpts);
 
    // NOTE: these are diagnostic only
    liq_tau_bnd_sw = real3d("liq_tau_bnd_sw", ncol, nlev, nswbands);
    ice_tau_bnd_sw = real3d("ice_tau_bnd_sw", ncol, nlev, nswbands);
    snw_tau_bnd_sw = real3d("snw_tau_bnd_sw", ncol, nlev, nswbands);
    liq_tau_bnd_lw = real3d("liq_tau_bnd_lw", ncol, nlev, nlwbands);
    ice_tau_bnd_lw = real3d("ice_tau_bnd_lw", ncol, nlev, nlwbands);
    snw_tau_bnd_lw = real3d("snw_tau_bnd_lw", ncol, nlev, nlwbands);
 
    // Gas volume mixing ratios
    gas_vmr = real3d("gas_vmr", ngas, ncol, nlev);
 
    // Needed for shortwave aerosol;
    //int nday, nnight;     // Number of daylight columns
    int1d day_indices("day_indices", ncol), night_indices("night_indices", ncol);   // Indices of daylight coumns
 
    // Flag to carry (QRS,QRL)*dp across time steps.
    // TODO: what does this mean?
    bool conserve_energy = true;
 
    // For loops over diagnostic calls
    //bool active_calls(0:N_DIAG)
 
    // Zero-array for cloud properties if not diagnosed by microphysics
    real2d zeros("zeros", ncol, nlev);
 
    gpoint_bands_sw = int1d("gpoint_bands_sw", nswgpts);
    gpoint_bands_lw = int1d("gpoint_bands_lw", nlwgpts);
 
    // Do shortwave stuff...
    if (do_short_wave_rad) {
        // Radiative fluxes
        internal::initial_fluxes(ncol, nlev+1, nswbands, sw_fluxes_allsky);
        internal::initial_fluxes(ncol, nlev+1, nswbands, sw_fluxes_clrsky);
 
        // TODO: Integrate calendar day computation
        int calday = 1;
        // Get cosine solar zenith angle for current time step.
        if (m_lat) {
            zenith(calday, m_lat, m_lon, rank_offsets, coszrs, ncol,
                   eccen,  mvelpp, lambm0, obliqr);
        } else {
            zenith(calday, m_lat, m_lon, rank_offsets, coszrs, ncol,
                   eccen,  mvelpp, lambm0, obliqr, uniform_angle);
        }
 
        // Get albedo. This uses CAM routines internally and just provides a
        // wrapper to improve readability of the code here.
        set_albedo(coszrs, albedo_dir, albedo_dif);
 
        // Do shortwave cloud optics calculations
        yakl::memset(cld_tau_gpt_sw, 0.);
        yakl::memset(cld_ssa_gpt_sw, 0.);
        yakl::memset(cld_asm_gpt_sw, 0.);
 
        // set cloud fraction to be 1, and snow fraction 0
        yakl::memset(cldfsnow, 0.0);
        yakl::memset(cld, 1.0);
 
        parallel_for (SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k)
        {
            iciwp(i,k) = std::min(qi(i,k)/std::max(1.0e-4,cld(i,k)),0.005)*pmid(i,k)/CONST_GRAV;
            iclwp(i,k) = std::min(qt(i,k)/std::max(1.0e-4,cld(i,k)),0.005)*pmid(i,k)/CONST_GRAV;
            icswp(i,k) = qn(i,k)/std::max(1.0e-4,cldfsnow(i,k))*pmid(i,k)/CONST_GRAV;
        });
 
        m2005_effradius(qc, qc, qi, qi, qt, qt, cld, pmid, tmid,
                        rel, rei, dei, lambdac, mu, des);
 
        // calculate the cloud radiation
        optics.get_cloud_optics_sw(ncol, nlev, nswbands, do_snow_optics, cld,
                                   cldfsnow, iclwp, iciwp, icswp,
                                   lambdac, mu, dei, des, rel, rei,
                                   cld_tau_bnd_sw, cld_ssa_bnd_sw, cld_asm_bnd_sw,
                                   liq_tau_bnd_sw, ice_tau_bnd_sw, snw_tau_bnd_sw);
 
        // Now reorder bands to be consistent with RRTMGP
        // We need to fix band ordering because the old input files assume RRTMG
        // band ordering, but this has changed in RRTMGP.
        // TODO: fix the input files themselves!
        cld_tau_bnd_sw_1d = real1d("cld_tau_bnd_sw_1d", nswbands);
        cld_ssa_bnd_sw_1d = real1d("cld_ssa_bnd_sw_1d", nswbands);
        cld_asm_bnd_sw_1d = real1d("cld_asm_bnd_sw_1d", nswbands);
        cld_tau_bnd_sw_o_1d = real1d("cld_tau_bnd_sw_1d", nswbands);
        cld_ssa_bnd_sw_o_1d = real1d("cld_ssa_bnd_sw_1d", nswbands);
        cld_asm_bnd_sw_o_1d = real1d("cld_asm_bnd_sw_1d", nswbands);
 
        parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilay)
        {
            for (auto ibnd = 1; ibnd <= nswbands; ++ibnd) {
                cld_tau_bnd_sw_1d(ibnd) = cld_tau_bnd_sw(icol,ilay,ibnd);
                cld_ssa_bnd_sw_1d(ibnd) = cld_ssa_bnd_sw(icol,ilay,ibnd);
                cld_asm_bnd_sw_1d(ibnd) = cld_asm_bnd_sw(icol,ilay,ibnd);
            }
            internal::reordered(cld_tau_bnd_sw_1d, rrtmg_to_rrtmgp, cld_tau_bnd_sw_o_1d);
            internal::reordered(cld_ssa_bnd_sw_1d, rrtmg_to_rrtmgp, cld_ssa_bnd_sw_o_1d);
            internal::reordered(cld_asm_bnd_sw_1d, rrtmg_to_rrtmgp, cld_asm_bnd_sw_o_1d);
            for (auto ibnd = 1; ibnd <= nswbands; ++ibnd) {
                cld_tau_bnd_sw(icol,ilay,ibnd) = cld_tau_bnd_sw_o_1d(ibnd);
                cld_ssa_bnd_sw(icol,ilay,ibnd) = cld_ssa_bnd_sw_o_1d(ibnd);
                cld_asm_bnd_sw(icol,ilay,ibnd) = cld_asm_bnd_sw_o_1d(ibnd);
            }
        });
 
        // And now do the MCICA sampling to get cloud optical properties by
        // gpoint/cloud state
        radiation.get_gpoint_bands_sw(gpoint_bands_sw);
 
        optics.sample_cloud_optics_sw(ncol, nlev, nswgpts, gpoint_bands_sw,
                                      pmid, cld, cldfsnow,
                                      cld_tau_bnd_sw, cld_ssa_bnd_sw, cld_asm_bnd_sw,
                                      cld_tau_gpt_sw, cld_ssa_gpt_sw, cld_asm_gpt_sw);
 
        // Aerosol needs night indices
        // TODO: remove this dependency, it's just used to mask aerosol outputs
        set_daynight_indices(coszrs, day_indices, night_indices);
        int1d nday("nday",1);
        int1d nnight("nnight",1);
        yakl::memset(nday, 0);
        yakl::memset(nnight, 0);
        for (auto icol=1; icol<=ncol; ++icol) {
            if (day_indices(icol) > 0) nday(1)++;
            if (night_indices(icol) > 0) nnight(1)++;
        }
 
        AMREX_ALWAYS_ASSERT(nday(1) + nnight(1) == ncol);
 
        // get aerosol optics
        do_aerosol_rad = false; // TODO: this causes issues if enabled
        {
            // Get gas concentrations
            get_gas_vmr(active_gases, gas_vmr);
 
            // Get aerosol optics
            if (do_aerosol_rad) {
                yakl::memset(aer_tau_bnd_sw, 0.);
                yakl::memset(aer_ssa_bnd_sw, 0.);
                yakl::memset(aer_asm_bnd_sw, 0.);
 
                clear_rh = real2d("clear_rh",ncol, nswbands);
                yakl::memset(clear_rh, 0.01);
 
                optics.set_aerosol_optics_sw(0, ncol, nlev, nswbands, dt, night_indices,
                                             is_cmip6_volc, aer_tau_bnd_sw, aer_ssa_bnd_sw, aer_asm_bnd_sw, clear_rh);
 
                // Now reorder bands to be consistent with RRTMGP
                // TODO: fix the input files themselves!
                aer_tau_bnd_sw_1d = real1d("cld_tau_bnd_sw_1d", nswbands);
                aer_ssa_bnd_sw_1d = real1d("cld_ssa_bnd_sw_1d", nswbands);
                aer_asm_bnd_sw_1d = real1d("cld_asm_bnd_sw_1d", nswbands);
                aer_tau_bnd_sw_o_1d = real1d("cld_tau_bnd_sw_1d", nswbands);
                aer_ssa_bnd_sw_o_1d = real1d("cld_ssa_bnd_sw_1d", nswbands);
                aer_asm_bnd_sw_o_1d = real1d("cld_asm_bnd_sw_1d", nswbands);
 
                parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilay)
                {
                    for (auto ibnd = 1; ibnd <= nswbands; ++ibnd) {
                        aer_tau_bnd_sw_1d(ibnd) = aer_tau_bnd_sw(icol,ilay,ibnd);
                        aer_ssa_bnd_sw_1d(ibnd) = aer_ssa_bnd_sw(icol,ilay,ibnd);
                        aer_asm_bnd_sw_1d(ibnd) = aer_asm_bnd_sw(icol,ilay,ibnd);
                    }
                    internal::reordered(aer_tau_bnd_sw_1d, rrtmg_to_rrtmgp, aer_tau_bnd_sw_o_1d);
                    internal::reordered(aer_ssa_bnd_sw_1d, rrtmg_to_rrtmgp, aer_ssa_bnd_sw_o_1d);
                    internal::reordered(aer_asm_bnd_sw_1d, rrtmg_to_rrtmgp, aer_asm_bnd_sw_o_1d);
                    for (auto ibnd = 1; ibnd <= nswbands; ++ibnd) {
                        aer_tau_bnd_sw(icol,ilay,ibnd) = aer_tau_bnd_sw_o_1d(ibnd);
                        aer_ssa_bnd_sw(icol,ilay,ibnd) = aer_ssa_bnd_sw_o_1d(ibnd);
                        aer_asm_bnd_sw(icol,ilay,ibnd) = aer_asm_bnd_sw_o_1d(ibnd);
                    }
                });
            } else {
                yakl::memset(aer_tau_bnd_sw, 0.);
                yakl::memset(aer_ssa_bnd_sw, 0.);
                yakl::memset(aer_asm_bnd_sw, 0.);
            }
 
         yakl::memset(cld_tau_gpt_sw, 0.);
         yakl::memset(cld_ssa_gpt_sw, 0.);
         yakl::memset(cld_asm_gpt_sw, 0.);
 
         // Call the shortwave radiation driver
         radiation_driver_sw(ncol, gas_vmr,
                             pmid, pint, tmid, albedo_dir, albedo_dif, coszrs,
                             cld_tau_gpt_sw, cld_ssa_gpt_sw, cld_asm_gpt_sw,
                             aer_tau_bnd_sw, aer_ssa_bnd_sw, aer_asm_bnd_sw,
                             sw_fluxes_allsky, sw_fluxes_clrsky, qrs, qrsc);
        }
 
        // Set surface fluxes that are used by the land model
        export_surface_fluxes(sw_fluxes_allsky, "shortwave");
 
    } else {
        // Conserve energy
        if (conserve_energy) {
            parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
            {
                qrs(icol,ilev) = qrs(icol,ilev)/pdel(icol,ilev);
            });
        }
    }  // dosw
 
    // Do longwave stuff...
    if (do_long_wave_rad) {
        // Allocate longwave outputs; why is this not part of the fluxes_t object?
        internal::initial_fluxes(ncol, nlev, nlwbands, lw_fluxes_allsky);
        internal::initial_fluxes(ncol, nlev, nlwbands, lw_fluxes_clrsky);
 
        // NOTE: fluxes defined at interfaces, so initialize to have vertical dimension nlev_rad+1
        yakl::memset(cld_tau_gpt_lw, 0.);
 
        optics.get_cloud_optics_lw(ncol, nlev, nlwbands, do_snow_optics, cld, cldfsnow, iclwp, iciwp, icswp,
                                   lambdac, mu, dei, des, rei,
                                   cld_tau_bnd_lw, liq_tau_bnd_lw, ice_tau_bnd_lw, snw_tau_bnd_lw);
 
        radiation.get_gpoint_bands_lw(gpoint_bands_lw);
 
        optics.sample_cloud_optics_lw(ncol, nlev, nlwgpts, gpoint_bands_lw,
                                      pmid, cld, cldfsnow,
                                      cld_tau_bnd_lw, cld_tau_gpt_lw);
 
        // Get gas concentrations
        get_gas_vmr(active_gases, gas_vmr);
 
        // Get aerosol optics
        yakl::memset(aer_tau_bnd_lw, 0.);
        if (do_aerosol_rad) {
            aer_rad.aer_rad_props_lw(is_cmip6_volc, 0, dt, zi, aer_tau_bnd_lw, clear_rh);
        }
 
        // Call the longwave radiation driver to calculate fluxes and heating rates
        radiation_driver_lw(ncol, nlev, gas_vmr, pmid, pint, tmid, tint, cld_tau_gpt_lw, aer_tau_bnd_lw,
                            lw_fluxes_allsky, lw_fluxes_clrsky, qrl, qrlc);
 
        // Set surface fluxes that are used by the land model
        export_surface_fluxes(lw_fluxes_allsky, "longwave");
 
    }
    else {
        // Conserve energy (what does this mean exactly?)
        if (conserve_energy) {
            parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
            {
                qrl(icol,ilev) = qrl(icol,ilev)/pdel(icol,ilev);
            });
        }
    } // dolw
 
    // 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);
        int const offset = rank_offsets[mfi.LocalIndex()];
        amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            // Map (col,lev) to (i,j,k)
            auto icol = (j-box3d.smallEnd(1))*nx + (i-box3d.smallEnd(0)) + 1 + offset;
            auto ilev = k+1;
 
            // TODO: We do not include the cloud source term qrsc/qrlc.
            //       Do these simply sum for a net source or do we pick one?
 
            // SW and LW sources
            qrad_src_array(i,j,k,0) = qrs(icol,ilev);
            qrad_src_array(i,j,k,1) = qrl(icol,ilev);
        });
    }
}

Here is the call graph for this function:

set_daynight_indices()

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

{
    // Loop over columns and identify daytime columns as those where the cosine
    // solar zenith angle exceeds zero. Note that we wrap the setting of
    // day_indices in an if-then to make sure we are not accessing day_indices out
    // of bounds, and stopping with an informative error message if we do for some reason.
    int1d iday("iday", 1);
    int1d inight("inight",1);
    yakl::memset(iday, 0);
    yakl::memset(inight, 0);
    yakl::memset(day_indices, 0);
    yakl::memset(night_indices, 0);
    parallel_for(SimpleBounds<1>(ncol), YAKL_LAMBDA (int icol)
    {
        if (coszrs(icol) > 0.) {
            iday(1) += 1;
            day_indices(iday(1)) = icol;
        } else {
            inight(1) += 1;
            night_indices(inight(1)) = icol;
        }
    });
}

writePlotfile()

void Radiation::writePlotfile	(	const std::string &	plot_prefix,
		const amrex::Real	time,
		const int	level_step
	)

{
    // Note: Radiation::initialize() is not called until the first timestep, so skip over the initial file write at t=0
    if (!qrad_src)
    {
        return;
    }
 
    std::string plotfilename = amrex::Concatenate(plot_prefix + "_rad", level_step, 5);
 
    // list of real2d (3D) variables to plot
    amrex::Vector<real2d*> plotvars_2d;
    plotvars_2d.push_back(&cld);
    plotvars_2d.push_back(&cldfsnow);
    plotvars_2d.push_back(&iclwp);
    plotvars_2d.push_back(&iciwp);
    plotvars_2d.push_back(&icswp);
    plotvars_2d.push_back(&dei);
    plotvars_2d.push_back(&des);
    plotvars_2d.push_back(&lambdac);
    plotvars_2d.push_back(&mu);
    plotvars_2d.push_back(&rei);
    plotvars_2d.push_back(&rel);
 
    //   SW allsky
    plotvars_2d.push_back(&sw_fluxes_allsky.flux_up);
    plotvars_2d.push_back(&sw_fluxes_allsky.flux_dn);
    plotvars_2d.push_back(&sw_fluxes_allsky.flux_net);
    plotvars_2d.push_back(&sw_fluxes_allsky.flux_dn_dir);
    //   SW clearsky
    plotvars_2d.push_back(&sw_fluxes_clrsky.flux_up);
    plotvars_2d.push_back(&sw_fluxes_clrsky.flux_dn);
    plotvars_2d.push_back(&sw_fluxes_clrsky.flux_net);
    plotvars_2d.push_back(&sw_fluxes_clrsky.flux_dn_dir);
 
    //   LW allsky
    plotvars_2d.push_back(&lw_fluxes_allsky.flux_up);
    plotvars_2d.push_back(&lw_fluxes_allsky.flux_dn);
    plotvars_2d.push_back(&lw_fluxes_allsky.flux_net);
    //   LW clearsky
    plotvars_2d.push_back(&lw_fluxes_clrsky.flux_up);
    plotvars_2d.push_back(&lw_fluxes_clrsky.flux_dn);
    plotvars_2d.push_back(&lw_fluxes_clrsky.flux_net);
 
    plotvars_2d.push_back(&qrs);
    plotvars_2d.push_back(&qrl);
    plotvars_2d.push_back(&zi);
    plotvars_2d.push_back(&clear_rh);
    plotvars_2d.push_back(&qrsc);
    plotvars_2d.push_back(&qrlc);
    plotvars_2d.push_back(&qt);
    plotvars_2d.push_back(&qi);
    plotvars_2d.push_back(&qc);
    plotvars_2d.push_back(&qn);
    plotvars_2d.push_back(&tmid);
    plotvars_2d.push_back(&pmid);
    plotvars_2d.push_back(&pdel);
 
    //plotvars_2d.push_back(&pint); // NOTE these have nlev + 1
    //plotvars_2d.push_back(&tint);
 
    //plotvars_2d.push_back(&albedo_dir); // [nswbands, ncol]
    //plotvars_2d.push_back(&albedo_dif);
 
    // names of plotted variables
    amrex::Vector<std::string> varnames_2d;
    varnames_2d.push_back("cld");
    varnames_2d.push_back("cldfsnow");
    varnames_2d.push_back("iclwp");
    varnames_2d.push_back("iciwp");
    varnames_2d.push_back("icswp");
    varnames_2d.push_back("dei");
    varnames_2d.push_back("des");
    varnames_2d.push_back("lambdac");
    varnames_2d.push_back("mu");
    varnames_2d.push_back("rei");
    varnames_2d.push_back("rel");
 
    //   SW allsky
    varnames_2d.push_back("sw_fluxes_allsky.flux_up");
    varnames_2d.push_back("sw_fluxes_allsky.flux_dn");
    varnames_2d.push_back("sw_fluxes_allsky.flux_net");
    varnames_2d.push_back("sw_fluxes_allsky.flux_dn_dir");
    //   SW clearsky
    varnames_2d.push_back("sw_fluxes_clrsky.flux_up");
    varnames_2d.push_back("sw_fluxes_clrsky.flux_dn");
    varnames_2d.push_back("sw_fluxes_clrsky.flux_net");
    varnames_2d.push_back("sw_fluxes_clrsky.flux_dn_dir");
 
    //   LW allsky
    varnames_2d.push_back("lw_fluxes_allsky.flux_up");
    varnames_2d.push_back("lw_fluxes_allsky.flux_dn");
    varnames_2d.push_back("lw_fluxes_allsky.flux_net");
    //   LW clearsky
    varnames_2d.push_back("lw_fluxes_clrsky.flux_up");
    varnames_2d.push_back("lw_fluxes_clrsky.flux_dn");
    varnames_2d.push_back("lw_fluxes_clrsky.flux_net");
 
    varnames_2d.push_back("qrs");
    varnames_2d.push_back("qrl");
    varnames_2d.push_back("zi");
    varnames_2d.push_back("clear_rh");
    varnames_2d.push_back("qrsc");
    varnames_2d.push_back("qrlc");
    varnames_2d.push_back("qt");
    varnames_2d.push_back("qi");
    varnames_2d.push_back("qc");
    varnames_2d.push_back("qn");
    varnames_2d.push_back("tmid");
    varnames_2d.push_back("pmid");
    varnames_2d.push_back("pdel");
 
    //varnames_2d.push_back("pint"); // NOTE these have nlev + 1
    //varnames_2d.push_back("tint");
 
    //varnames_2d.push_back("albedo_dir");
    //varnames_2d.push_back("albedo_dif");
 
    // list 3D variables defined over bands (SW and LW)
    //   these are 4D fields split into 3D so they can use the same AMReX plotfile
    amrex::Vector<real3d*> plotvars_3d;
    amrex::Vector<std::string> varnames_3d;
    plotvars_3d.push_back(&cld_tau_bnd_sw);
    plotvars_3d.push_back(&cld_ssa_bnd_sw);
    plotvars_3d.push_back(&cld_asm_bnd_sw);
    plotvars_3d.push_back(&aer_tau_bnd_sw);
    plotvars_3d.push_back(&aer_ssa_bnd_sw);
    plotvars_3d.push_back(&aer_asm_bnd_sw);
    plotvars_3d.push_back(&cld_tau_bnd_lw);
    plotvars_3d.push_back(&aer_tau_bnd_lw);
    //   diagnostic variables:
    plotvars_3d.push_back(&liq_tau_bnd_sw);
    plotvars_3d.push_back(&ice_tau_bnd_sw);
    plotvars_3d.push_back(&snw_tau_bnd_sw);
    plotvars_3d.push_back(&liq_tau_bnd_lw);
    plotvars_3d.push_back(&ice_tau_bnd_lw);
    plotvars_3d.push_back(&snw_tau_bnd_lw);
 
    varnames_3d.push_back("cld_tau_bnd_sw");
    varnames_3d.push_back("cld_ssa_bnd_sw");
    varnames_3d.push_back("cld_asm_bnd_sw");
    varnames_3d.push_back("aer_tau_bnd_sw");
    varnames_3d.push_back("aer_ssa_bnd_sw");
    varnames_3d.push_back("aer_asm_bnd_sw");
    varnames_3d.push_back("cld_tau_bnd_lw");
    varnames_3d.push_back("aer_tau_bnd_lw");
    //   diagnostic variables:
    varnames_3d.push_back("liq_tau_bnd_sw");
    varnames_3d.push_back("ice_tau_bnd_sw");
    varnames_3d.push_back("snw_tau_bnd_sw");
    varnames_3d.push_back("liq_tau_bnd_lw");
    varnames_3d.push_back("ice_tau_bnd_lw");
    varnames_3d.push_back("snw_tau_bnd_lw");
 
    AMREX_ASSERT(varnames_2d.size() == plotvars_2d.size());
    AMREX_ASSERT(varnames_3d.size() == plotvars_3d.size());
 
    int output_size = plotvars_2d.size() + 1; // 2D vars + coszrs
    //  add in total output size of all 3D vars
    for (int i = 0; i < plotvars_3d.size(); i++)
    {
        output_size += plotvars_3d[i]->get_dimensions()(3);
    }
 
    // convert each YAKL array to a MF and combine into a single MF for plotting
    MultiFab fab(m_box, qrad_src->DistributionMap(), output_size, 0);
    // copy 2D vars first
    for (int i = 0; i < plotvars_2d.size(); i++)
    {
        amrex::MultiFab mf;
        yakl_to_mf(*plotvars_2d[i], mf);
        MultiFab::Copy(fab, mf, 0, i, 1, 0);
    }
 
    int dst = plotvars_2d.size(); // output component index
 
    // expand coszrs from 2D to 3D so it can be saved with the same file
    varnames_2d.push_back("coszrs");
    amrex::MultiFab coszrs_mf;
    expand_yakl1d_to_mf(coszrs, coszrs_mf);
    MultiFab::Copy(fab, coszrs_mf, 0, dst, 1, 0);
    dst++;
 
    // copy 3D vars
    for (int i = 0; i < plotvars_3d.size(); i++)
    {
        // split real3ds into real2d for each band for output
        // TODO: find a better way to do this
        const int var_nbands = plotvars_3d[i]->get_dimensions()(3); // either nswbands or nlwbands
        for (int bnd = 1; bnd <= var_nbands; bnd++)
        {
            // add variable name to 2d list
            varnames_2d.push_back(varnames_3d[i] + "_" + std::to_string(bnd));
 
            amrex::MultiFab mf;
            real2d band_var = plotvars_3d[i]->slice<2>(yakl::COLON, yakl::COLON, bnd);
            yakl_to_mf(band_var, mf);
            MultiFab::Copy(fab, mf, 0, dst, 1, 0);
 
            dst++;
        }
    }
    AMREX_ASSERT(dst == fab.nComp());
 
    // this should now match full output size with all 3D variables expanded
    AMREX_ASSERT(varnames_2d.size() == output_size);
 
 
    amrex::WriteSingleLevelPlotfile(plotfilename, fab, varnames_2d, m_geom, time, level_step);
}

yakl_to_mf()

void Radiation::yakl_to_mf	(	const real2d &	data,
		amrex::MultiFab &	mf
	)

{
    // creates a MF from a YAKL real2d mapping from [col, lev] to [x,y,z]
    // by reshaping the yakl array to the geometry of the output qsrc multifab
    mf = amrex::MultiFab(m_box, qrad_src->DistributionMap(), 1, 0);
    if (!data.initialized())
    {
        // yakl array hasn't been created yet, so create an empty MF
        mf.setVal(0.0);
        return;
    }
 
    for (MFIter mfi(mf); mfi.isValid(); ++mfi) {
        auto mf_arr = mf.array(mfi);
        const auto& box3d = mfi.tilebox();
        const int nx = box3d.length(0);
        const int offset = rank_offsets[mfi.LocalIndex()];
        amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            // map [i,j,k] 0-based to [icol, ilev] 1-based
            const int icol = (j-box3d.smallEnd(1))*nx + (i-box3d.smallEnd(0)) + 1 + offset;
            const int ilev = k+1;
            AMREX_ASSERT(icol <= static_cast<int>(data.get_dimensions()(1)));
            AMREX_ASSERT(ilev <= static_cast<int>(data.get_dimensions()(2)));
            mf_arr(i, j, k) = data(icol, ilev);
        });
    }
}

Here is the call graph for this function:

Member Data Documentation

active_gases

const std::vector<std::string> Radiation::active_gases

private

Initial value:

= {
                                                   "H2O", "CO2", "O3", "N2O",
                                                   "CO" , "CH4", "O2", "N2" }

aer_asm_bnd_sw

real3d Radiation::aer_asm_bnd_sw

private

aer_asm_bnd_sw_1d

real1d Radiation::aer_asm_bnd_sw_1d

private

aer_asm_bnd_sw_o_1d

real1d Radiation::aer_asm_bnd_sw_o_1d

private

aer_rad

AerRadProps Radiation::aer_rad

private

aer_ssa_bnd_sw

real3d Radiation::aer_ssa_bnd_sw

private

aer_ssa_bnd_sw_1d

real1d Radiation::aer_ssa_bnd_sw_1d

private

aer_ssa_bnd_sw_o_1d

real1d Radiation::aer_ssa_bnd_sw_o_1d

private

aer_tau_bnd_lw

real3d Radiation::aer_tau_bnd_lw

private

aer_tau_bnd_sw

real3d Radiation::aer_tau_bnd_sw

private

aer_tau_bnd_sw_1d

real1d Radiation::aer_tau_bnd_sw_1d

private

aer_tau_bnd_sw_o_1d

real1d Radiation::aer_tau_bnd_sw_o_1d

private

aernames

std::vector<std::string> Radiation::aernames

private

albedo_dif

real2d Radiation::albedo_dif

private

albedo_dir

real2d Radiation::albedo_dir

private

cld

real2d Radiation::cld

private

cld_asm_bnd_sw

real3d Radiation::cld_asm_bnd_sw

private

cld_asm_bnd_sw_1d

real1d Radiation::cld_asm_bnd_sw_1d

private

cld_asm_bnd_sw_o_1d

real1d Radiation::cld_asm_bnd_sw_o_1d

private

cld_asm_gpt_sw

real3d Radiation::cld_asm_gpt_sw

private

cld_ssa_bnd_sw

real3d Radiation::cld_ssa_bnd_sw

private

cld_ssa_bnd_sw_1d

real1d Radiation::cld_ssa_bnd_sw_1d

private

cld_ssa_bnd_sw_o_1d

real1d Radiation::cld_ssa_bnd_sw_o_1d

private

cld_ssa_gpt_sw

real3d Radiation::cld_ssa_gpt_sw

private

cld_tau_bnd_lw

real3d Radiation::cld_tau_bnd_lw

private

cld_tau_bnd_sw

real3d Radiation::cld_tau_bnd_sw

private

cld_tau_bnd_sw_1d

real1d Radiation::cld_tau_bnd_sw_1d

private

cld_tau_bnd_sw_o_1d

real1d Radiation::cld_tau_bnd_sw_o_1d

private

cld_tau_gpt_lw

real3d Radiation::cld_tau_gpt_lw

private

cld_tau_gpt_sw

real3d Radiation::cld_tau_gpt_sw

private

cldfsnow

real2d Radiation::cldfsnow

private

clear_rh

real2d Radiation::clear_rh

private

coszrs

real1d Radiation::coszrs

private

dei

real2d Radiation::dei

private

des

real2d Radiation::des

private

do_aerosol_rad

bool Radiation::do_aerosol_rad = true

private

do_long_wave_rad

bool Radiation::do_long_wave_rad

private

do_short_wave_rad

bool Radiation::do_short_wave_rad

private

do_snow_optics

bool Radiation::do_snow_optics

private

dohirs

bool Radiation::dohirs = false

private

dt

real Radiation::dt

private

dt_avg

real Radiation::dt_avg = 0.0

private

eccen

constexpr amrex::Real Radiation::eccen = 0.0

staticconstexprprivate

fixed_total_solar_irradiance

real Radiation::fixed_total_solar_irradiance = -1.

private

flns

real1d Radiation::flns

private

flnt

real1d Radiation::flnt

private

fsds

real1d Radiation::fsds

private

fsns

real1d Radiation::fsns

private

fsnt

real1d Radiation::fsnt

private

gas_vmr

real3d Radiation::gas_vmr

private

gasnames

std::vector<std::string> Radiation::gasnames

private

gpoint_bands_lw

int1d Radiation::gpoint_bands_lw

private

gpoint_bands_sw

int1d Radiation::gpoint_bands_sw

private

has_qmoist

bool Radiation::has_qmoist

private

ice_tau_bnd_lw

real3d Radiation::ice_tau_bnd_lw

private

ice_tau_bnd_sw

real3d Radiation::ice_tau_bnd_sw

private

iciwp

real2d Radiation::iciwp

private

iclwp

real2d Radiation::iclwp

private

icswp

real2d Radiation::icswp

private

ihirsfq

int Radiation::ihirsfq = 1

private

is_cmip6_volc

bool Radiation::is_cmip6_volc

private

lambdac

real2d Radiation::lambdac

private

lambm0

constexpr amrex::Real Radiation::lambm0 = -3.2503635878519378e-2

staticconstexprprivate

liq_tau_bnd_lw

real3d Radiation::liq_tau_bnd_lw

private

liq_tau_bnd_sw

real3d Radiation::liq_tau_bnd_sw

private

lw_band_midpoints

real1d Radiation::lw_band_midpoints

private

lw_fluxes_allsky

FluxesByband Radiation::lw_fluxes_allsky

private

lw_fluxes_clrsky

FluxesByband Radiation::lw_fluxes_clrsky

private

m_box

amrex::BoxArray Radiation::m_box

private

m_geom

amrex::Geometry Radiation::m_geom

private

m_lat

amrex::MultiFab* Radiation::m_lat = nullptr

private

m_lon

amrex::MultiFab* Radiation::m_lon = nullptr

private

m_lsm_fluxes

amrex::MultiFab* Radiation::m_lsm_fluxes = nullptr

private

m_lsm_zenith

amrex::MultiFab* Radiation::m_lsm_zenith = nullptr

private

moisture_type

std::string Radiation::moisture_type = "None"

private

mu

real2d Radiation::mu

private

mvelpp

constexpr amrex::Real Radiation::mvelpp = 103.0 * PI / 180.0 + PI

staticconstexprprivate

naer

int Radiation::naer

private

ncol

int Radiation::ncol

private

net_flux

real1d Radiation::net_flux

private

ngas

int Radiation::ngas

private

nlev

int Radiation::nlev

private

nlwbands

int Radiation::nlwbands

private

nlwgpts

int Radiation::nlwgpts

private

nswbands

int Radiation::nswbands

private

nswgpts

int Radiation::nswgpts

private

obliqr

constexpr amrex::Real Radiation::obliqr = 23.0 * PI / 180.0

staticconstexprprivate

optics

Optics Radiation::optics

private

pdel

real2d Radiation::pdel

private

pint

real2d Radiation::pint

private

pmid

real2d Radiation::pmid

private

qc

real2d Radiation::qc

private

qi

real2d Radiation::qi

private

qn

real2d Radiation::qn

private

qrad_src

amrex::MultiFab* Radiation::qrad_src = nullptr

private

qrl

real2d Radiation::qrl

private

qrlc

real2d Radiation::qrlc

private

qrs

real2d Radiation::qrs

private

qrsc

real2d Radiation::qrsc

private

qt

real2d Radiation::qt

private

radiation

Rrtmgp Radiation::radiation

private

rank_offsets

amrex::Vector<int> Radiation::rank_offsets

private

rei

real2d Radiation::rei

private

rel

real2d Radiation::rel

private

rrtmg_to_rrtmgp

int1d Radiation::rrtmg_to_rrtmgp

private

rrtmgp_coefficients_file_lw

std::string Radiation::rrtmgp_coefficients_file_lw

private

rrtmgp_coefficients_file_name_lw

std::string Radiation::rrtmgp_coefficients_file_name_lw = "rrtmgp_coefficients_lw_20181204.nc"

private

rrtmgp_coefficients_file_name_sw

std::string Radiation::rrtmgp_coefficients_file_name_sw = "rrtmgp_coefficients_sw_20181204.nc"

private

rrtmgp_coefficients_file_sw

std::string Radiation::rrtmgp_coefficients_file_sw

private

rrtmgp_data_path

std::string Radiation::rrtmgp_data_path

private

rrtmgp_enable_temperature_warnings

bool Radiation::rrtmgp_enable_temperature_warnings = true

private

snw_tau_bnd_lw

real3d Radiation::snw_tau_bnd_lw

private

snw_tau_bnd_sw

real3d Radiation::snw_tau_bnd_sw

private

spectralflux

bool Radiation::spectralflux = false

private

sw_band_midpoints

real1d Radiation::sw_band_midpoints

private

sw_fluxes_allsky

FluxesByband Radiation::sw_fluxes_allsky

private

sw_fluxes_clrsky

FluxesByband Radiation::sw_fluxes_clrsky

private

tint

real2d Radiation::tint

private

tmid

real2d Radiation::tmid

private

uniform_angle

amrex::Real Radiation::uniform_angle = 78.463

private

use_rad_dt_cosz

bool Radiation::use_rad_dt_cosz = false

private

zhi

int Radiation::zhi

private

zi

real2d Radiation::zi

private

zlo

int Radiation::zlo

private

---

The documentation for this class was generated from the following files:

• Source/Radiation/ERF_Radiation.H
• Source/Radiation/ERF_Radiation.cpp