diff --git a/Source/Radiation/Radiation.H b/Source/Radiation/Radiation.H index 8151262bb..8284a97fa 100644 --- a/Source/Radiation/Radiation.H +++ b/Source/Radiation/Radiation.H @@ -2,13 +2,13 @@ * RTE-RRTMGP radiation model interface to ERF * The orginal code is developed by RobertPincus, and the code is open source available at: * https://github.com/earth-system-radiation/rte-rrtmgp - * Please reference to the following paper, + * Please reference to the following paper, * https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019MS001621 * NOTE: we use the C++ version of RTE-RRTMGP, which is the implemention of the original Fortran * code using C++ YAKL for CUDA, HiP and SYCL application by E3SM ECP team, the C++ version * of the rte-rrtmgp code is located at: * https://github.com/E3SM-Project/rte-rrtmgp - * The RTE-RRTMGP uses BSD-3-Clause Open Source License, if you want to make changes, + * The RTE-RRTMGP uses BSD-3-Clause Open Source License, if you want to make changes, * and modifications to the code, please refer to BSD-3-Clause Open Source License. */ #ifndef ERF_RADIATION_H @@ -51,7 +51,7 @@ class Radiation { // run radiation model void run(); - // call back + // call back void on_complete(); void radiation_driver_lw(int ncol, int nlev, @@ -62,23 +62,23 @@ class Radiation { real2d& qrl, real2d& qrlc); void radiation_driver_sw(int ncol, - real3d& gas_vmr, real2d& pmid, real2d& pint, real2d& tmid, + real3d& gas_vmr, real2d& pmid, real2d& pint, real2d& tmid, real2d& albedo_dir, real2d& albedo_dif, real1d& coszrs, real3d& cld_tau_gpt, real3d& cld_ssa_gpt, real3d& cld_asm_gpt, real3d& aer_tau_bnd, real3d& aer_ssa_bnd, real3d& aer_asm_bnd, FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, real2d& qrs, real2d& qrsc); - void set_daynight_indices(real1d& coszrs, - int1d& day_indices, + void set_daynight_indices(real1d& coszrs, + int1d& day_indices, int1d& night_indices); - void get_gas_vmr(const std::vector& gas_names, + void get_gas_vmr(const std::vector& gas_names, real3d& gas_vmr); - void calculate_heating_rate(real2d& flux_up, - real2d& flux_dn, - real2d& pint, + void calculate_heating_rate(real2d& flux_up, + real2d& flux_dn, + real2d& pint, real2d& heating_rate); private: // geometry @@ -91,7 +91,7 @@ class Radiation { int nlev, zlo, zhi; // number of columns in horizontal plane - int ncol; + int ncol; int nlwgpts, nswgpts; int nlwbands, nswbands; @@ -101,11 +101,11 @@ class Radiation { bool do_long_wave_rad; bool do_snow_optics; - // Flag to indicate whether to do aerosol optical calculations. This + // Flag to indicate whether to do aerosol optical calculations. This // zeroes out the aerosol optical properties if False bool do_aerosol_rad = true; - // rrtmgp + // rrtmgp Rrtmgp radiation; // optics radation properties @@ -121,7 +121,7 @@ class Radiation { real1d net_flux; // This should be module data or something specific to aerosol where it is used? - bool is_cmip6_volc; // true if cmip6 style volcanic file is read otherwise false + bool is_cmip6_volc; // true if cmip6 style volcanic file is read otherwise false real dt; // time step(s) - needed for aerosol optics call @@ -131,33 +131,33 @@ class Radiation { real1d flns; //(ncols) ! Srf longwave cooling (up-down) flux real1d flnt; //(ncols) ! Net outgoing lw flux at model top real1d fsds; //(ncols) ! Surface solar down flux - + // radiation data const std::vector active_gases = { - "H2O", "CO2", "O3", "N2O", + "H2O", "CO2", "O3", "N2O", "CO", "CH4", "O2", "N2" }; bool spectralflux = false; // calculate fluxes (up and down) per band. // Flag to indicate whether or not to use the radiation timestep for solar zenith - // angle calculations. If true, use the radiation timestep for all solar zenith + // angle calculations. If true, use the radiation timestep for all solar zenith // angle (cosz) calculations. // TODO: How does this differ if value is .false.? bool use_rad_dt_cosz = false; - // Value for prescribing an invariant solar constant (i.e. total solar + // Value for prescribing an invariant solar constant (i.e. total solar // irradiance at TOA). Used for idealized experiments such as RCE. // Disabled when value is less than 0. real fixed_total_solar_irradiance = -1.; // The RRTMGP warnings are printed when the state variables need to be limted, // such as when the temperature drops too low. This is not normally an issue, - // but in aquaplanet and RCE configurations these situations occur much more + // but in aquaplanet and RCE configurations these situations occur much more // frequently, so this flag was added to be able to disable those messages. bool rrtmgp_enable_temperature_warnings = true; // Output diagnostic brightness temperatures at the top of the - // atmosphere for 7 TOVS/HIRS channels (2,4,6,8,10,11,12) and 4 TOVS/MSU + // atmosphere for 7 TOVS/HIRS channels (2,4,6,8,10,11,12) and 4 TOVS/MSU // channels (1,2,3,4). // TODO: where are these options set? bool dohirs = false; @@ -183,10 +183,10 @@ class Radiation { std::vector aernames; int1d rrtmg_to_rrtmgp; - + // Pointers to heating rates on physics buffer - real2d qrs; // shortwave radiative heating rate - real2d qrl; // longwave radiative heating rate + real2d qrs; // shortwave radiative heating rate + real2d qrl; // longwave radiative heating rate // Pointers to fields on the physics buffer real2d zi; @@ -194,7 +194,7 @@ class Radiation { // 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 - real2d qrsc; + real2d qrsc; real2d qrlc; real2d tmid, pmid, pdel; diff --git a/Source/Radiation/Radiation.cpp b/Source/Radiation/Radiation.cpp index b21f00c7e..f7170f8d1 100644 --- a/Source/Radiation/Radiation.cpp +++ b/Source/Radiation/Radiation.cpp @@ -2,13 +2,13 @@ * RTE-RRTMGP radiation model interface to ERF * The orginal code is developed by RobertPincus, and the code is open source available at: * https://github.com/earth-system-radiation/rte-rrtmgp - * Please reference to the following paper, + * Please reference to the following paper, * https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019MS001621 * NOTE: we use the C++ version of RTE-RRTMGP, which is reimplemented the original Fortran * code using C++ YAKL for CUDA, HiP and SYCL application by E3SM ECP team, the C++ version * of the rte-rrtmgp code is located at: * https://github.com/E3SM-Project/rte-rrtmgp - * The RTE-RRTMGP uses BSD-3-Clause Open Source License, if you want to make changes, + * The RTE-RRTMGP uses BSD-3-Clause Open Source License, if you want to make changes, * and modifications to the code, please refer to BSD-3-Clause Open Source License. */ #include @@ -84,10 +84,10 @@ void Radiation::initialize(const MultiFab& cons_in, MultiFab& qmoist, m_geom = geom; m_box = grids; - + auto dz = m_geom.CellSize(2); auto lowz = m_geom.ProbLo(2); - + dt = dt_advance; do_short_wave_rad = do_sw_rad; @@ -106,12 +106,12 @@ void Radiation::initialize(const MultiFab& cons_in, MultiFab& qmoist, // initialize cloud, aerosol, and radiation optics.initialize(); - radiation.initialize(ngas, active_gases, - rrtmgp_coefficients_file_sw.c_str(), + radiation.initialize(ngas, active_gases, + rrtmgp_coefficients_file_sw.c_str(), rrtmgp_coefficients_file_lw.c_str()); // initialize the radiation data - auto nswbands = radiation.get_nband_sw(); + auto nswbands = radiation.get_nband_sw(); auto nswgpts = radiation.get_ngpt_sw(); auto nlwbands = radiation.get_nband_lw(); auto nlwgpts = radiation.get_ngpt_lw(); @@ -132,10 +132,10 @@ void Radiation::initialize(const MultiFab& cons_in, MultiFab& qmoist, tint = real2d("tint", ncol, nlev+1); albedo_dir = real2d("albedo_dir", nswbands, ncol); - albedo_dif = real2d("albedo_dif", 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 + 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 @@ -153,7 +153,7 @@ void Radiation::initialize(const MultiFab& cons_in, MultiFab& qmoist, << " Fixed solar consant (disabled with -1): " //f10.4/, & << " Enable temperature warnings: "; //,l5/ ) -} +} // run radiation model @@ -164,7 +164,7 @@ void Radiation::run() { // Cosine solar zenith angle for all columns in chunk real1d coszrs("coszrs", ncol); - + // Pointers to fields on the physics buffer real2d cld("cld", ncol, nlev), cldfsnow("cldfsnow", ncol, nlev), iclwp("iclwp", ncol, nlev), iciwp("iciwp", ncol, nlev), @@ -173,25 +173,25 @@ void Radiation::run() { mu("mu", ncol, nlev), rei("rei", ncol, nlev), rel("rel", ncol, nlev); // Cloud, snow, and aerosol optical properties - real3d cld_tau_gpt_sw("cld_tau_gpt_sw", ncol, nlev, nswgpts), - cld_ssa_gpt_sw("cld_ssa_gpt_sw", ncol, nlev, nswgpts), + real3d cld_tau_gpt_sw("cld_tau_gpt_sw", ncol, nlev, nswgpts), + cld_ssa_gpt_sw("cld_ssa_gpt_sw", ncol, nlev, nswgpts), cld_asm_gpt_sw("cld_asm_gpt_sw", ncol, nlev, nswgpts); - real3d cld_tau_bnd_sw("cld_tau_bnd_sw", ncol, nlev, nswbands), - cld_ssa_bnd_sw("cld_ssa_bnd_sw", ncol, nlev, nswbands), + real3d cld_tau_bnd_sw("cld_tau_bnd_sw", ncol, nlev, nswbands), + cld_ssa_bnd_sw("cld_ssa_bnd_sw", ncol, nlev, nswbands), cld_asm_bnd_sw("cld_asm_bnd_sw", ncol, nlev, nswbands); - real3d aer_tau_bnd_sw("aer_tau_bnd_sw", ncol, nlev, nswbands), - aer_ssa_bnd_sw("aer_ssa_bnd_sw", ncol, nlev, nswbands), + real3d aer_tau_bnd_sw("aer_tau_bnd_sw", ncol, nlev, nswbands), + aer_ssa_bnd_sw("aer_ssa_bnd_sw", ncol, nlev, nswbands), aer_asm_bnd_sw("aer_asm_bnd_sw", ncol, nlev, nswbands); - real3d cld_tau_bnd_lw("cld_tau_bnd_lw", ncol, nlev, nlwbands), + real3d cld_tau_bnd_lw("cld_tau_bnd_lw", ncol, nlev, nlwbands), aer_tau_bnd_lw("aer_tau_bnd_lw", ncol, nlev, nlwbands); real3d cld_tau_gpt_lw("cld_tau_gpt_lw", ncol, nlev, nlwgpts); // NOTE: these are diagnostic only - real3d liq_tau_bnd_sw("liq_tau_bnd_sw", ncol, nlev, nswbands), - ice_tau_bnd_sw("ice_tau_bnd_sw", ncol, nlev, nswbands), + real3d liq_tau_bnd_sw("liq_tau_bnd_sw", ncol, nlev, nswbands), + ice_tau_bnd_sw("ice_tau_bnd_sw", ncol, nlev, nswbands), snw_tau_bnd_sw("snw_tau_bnd_sw", ncol, nlev, nswbands); - real3d liq_tau_bnd_lw("liq_tau_bnd_lw", ncol, nlev, nlwbands), - ice_tau_bnd_lw("ice_tau_bnd_lw", ncol, nlev, nlwbands), + real3d liq_tau_bnd_lw("liq_tau_bnd_lw", ncol, nlev, nlwbands), + ice_tau_bnd_lw("ice_tau_bnd_lw", ncol, nlev, nlwbands), snw_tau_bnd_lw("snw_tau_bnd_lw", ncol, nlev, nlwbands); // Gas volume mixing ratios @@ -201,7 +201,7 @@ void Radiation::run() { int nday, nnight; // Number of daylight columns int1d day_indices("day_indices", ncol), night_indices("night_indices", ncol); // Indicies of daylight coumns - // Flag to carry (QRS,QRL)*dp across time steps. + // Flag to carry (QRS,QRL)*dp across time steps. // TODO: what does this mean? bool conserve_energy = true; @@ -221,7 +221,7 @@ void Radiation::run() { internal::initial_fluxes(ncol, nlev, nlwbands, fluxes_allsky); internal::initial_fluxes(ncol, nlev, nlwbands, fluxes_clrsky); - // Get cosine solar zenith angle for current time step. ( still NOT YET implemented here) + // Get cosine solar zenith angle for current time step. ( still NOT YET implemented here) // set_cosine_solar_zenith_angle(state, dt_avg, coszrs(1:ncol)) yakl::memset(coszrs, 0.2); // we set constant value here to avoid numerical overflow. @@ -234,10 +234,10 @@ void Radiation::run() { yakl::memset(cld_ssa_gpt_sw, 0.); yakl::memset(cld_asm_gpt_sw, 0.); - 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, + 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 @@ -264,7 +264,7 @@ void Radiation::run() { 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); - } + } } } @@ -272,10 +272,10 @@ void Radiation::run() { // 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, + 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 @@ -325,7 +325,7 @@ void Radiation::run() { 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 { @@ -335,16 +335,16 @@ void Radiation::run() { } // 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, - fluxes_allsky, fluxes_clrsky, qrs, qrsc + 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, + fluxes_allsky, fluxes_clrsky, qrs, qrsc ); } - } + } else { // Conserve energy if (conserve_energy) { @@ -363,22 +363,22 @@ void Radiation::run() { FluxesByband fluxes_allsky, fluxes_clrsky; internal::initial_fluxes(ncol, nlev, nlwbands, fluxes_allsky); internal::initial_fluxes(ncol, nlev, nlwbands, 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, + 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, + ncol, nlev, nlwgpts, gpoint_bands_lw, + pmid, cld, cldfsnow, cld_tau_bnd_lw, cld_tau_gpt_lw); // Loop over diagnostic calls @@ -413,7 +413,7 @@ void Radiation::run() { } // dolw // Compute heating rate for dtheta/dt - for (auto ilay = 1; ilay < nlev; ++nlev) { + for (auto ilay = 1; ilay < nlev; ++nlev) { for (auto icol = 1; icol < ncol; ++icol) { hr(icol,ilay) = (qrs(icol,ilay) + qrl(icol,ilay)) / Cp_d * (1.e5 / std::pow(pmid(icol,ilay), R_d/Cp_d)); } @@ -428,18 +428,18 @@ void Radiation::run() { } } } -} +} -void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, - real2d& pmid, real2d& pint, real2d& tmid, real2d& albedo_dir, real2d& albedo_dif, - real1d& coszrs, real3d& cld_tau_gpt, real3d& cld_ssa_gpt, real3d& cld_asm_gpt, +void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, + real2d& pmid, real2d& pint, real2d& tmid, real2d& albedo_dir, real2d& albedo_dif, + real1d& coszrs, real3d& cld_tau_gpt, real3d& cld_ssa_gpt, real3d& cld_asm_gpt, real3d& aer_tau_bnd, real3d& aer_ssa_bnd, real3d& aer_asm_bnd, FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, real2d& qrs, real2d& qrsc) { - // Incoming solar radiation, scaled for solar zenith angle + // Incoming solar radiation, scaled for solar zenith angle // and earth-sun distance real2d solar_irradiance_by_gpt("solar_irradiance_by_gpt",ncol,nswgpts); - // Gathered indicies of day and night columns + // Gathered indicies of day and night columns // chunk_column_index = day_indices(daylight_column_index) int1d day_indices("day_indices",ncol), night_indices("night_indices", ncol); // Indicies of daylight coumns @@ -449,14 +449,14 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, 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 gas_vmr_day("gas_vmr_day", ngas, ncol, nlev); real3d gas_vmr_rad("gas_vmr_rad", ngas, ncol, nlev); - real3d cld_tau_gpt_day("cld_tau_gpt_day", ncol, nlev-1, nswgpts); - real3d cld_ssa_gpt_day("cld_ssa_gpt_day", ncol, nlev-1, nswgpts); + real3d cld_tau_gpt_day("cld_tau_gpt_day", ncol, nlev-1, nswgpts); + real3d cld_ssa_gpt_day("cld_ssa_gpt_day", ncol, nlev-1, nswgpts); real3d cld_asm_gpt_day("cld_asm_gpt_day", ncol, nlev-1, nswgpts); real3d aer_tau_bnd_day("aer_tau_bnd_day", ncol, nlev-1, nswbands); - real3d aer_ssa_bnd_day("aer_ssa_bnd_day", ncol, nlev-1, nswbands); + real3d aer_ssa_bnd_day("aer_ssa_bnd_day", ncol, nlev-1, nswbands); real3d aer_asm_bnd_day("aer_asm_bnd_day", ncol, nlev-1, nswbands); real3d cld_tau_gpt_rad("cld_tau_gpt_rad", ncol, nlev+1, nswgpts); @@ -477,7 +477,7 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, // tsi_scaling = get_eccentricity_factor(); } 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 + // 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; } @@ -551,13 +551,13 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, }); // Do shortwave radiative transfer calculations - radiation.run_shortwave_rrtmgp( ngas, nday, nlev, - gas_vmr_rad, pmid_day, tmid_day, pint_day, coszrs_day, albedo_dir_day, albedo_dif_day, + radiation.run_shortwave_rrtmgp( ngas, nday, nlev, + gas_vmr_rad, 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, + 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 @@ -565,8 +565,8 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, 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, + calculate_heating_rate(fluxes_allsky.flux_up, + fluxes_allsky.flux_dn, pint, qrs); calculate_heating_rate(fluxes_clrsky.flux_up, @@ -574,10 +574,10 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, pint, qrsc); } -void Radiation::radiation_driver_lw(int ncol, int nlev, - real3d& gas_vmr, - real2d& pmid, real2d& pint, real2d& tmid, real2d& tint, - real3d& cld_tau_gpt, real3d& aer_tau_bnd, FluxesByband& fluxes_clrsky, +void Radiation::radiation_driver_lw(int ncol, int nlev, + real3d& gas_vmr, + real2d& pmid, real2d& pint, real2d& tmid, real2d& tint, + real3d& cld_tau_gpt, real3d& aer_tau_bnd, FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, real2d& qrl, 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); @@ -611,12 +611,12 @@ void Radiation::radiation_driver_lw(int ncol, int nlev, }); // 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, + 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.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 @@ -647,7 +647,7 @@ void Radiation::set_daynight_indices(real1d& coszrs, int1d& day_indices, int1d& // of bounds, and stopping with an informative error message if we do for some reason. int iday = 0; int inight = 0; - for (auto icol = 0; icol < ncol; ++icol) { + for (auto icol = 0; icol < ncol; ++icol) { if (coszrs(icol) > 0.) { iday += 1; day_indices(iday) = icol; @@ -666,7 +666,7 @@ void Radiation::get_gas_vmr(const std::vector& gas_names, real3d& g // CO and N2, which RRTMG did not have. const std::vector gas_species = {"H2O", "CO2", "O3", "N2O", "CO", "CH4", "O2", "N2"}; - const std::vector mol_weight_gas = {18.01528, 44.0095, 47.9982, 44.0128, + const std::vector mol_weight_gas = {18.01528, 44.0095, 47.9982, 44.0128, 28.0101, 16.04246, 31.998, 28.0134}; // g/mol // Molar weight of air const real mol_weight_air = 28.97; // g/mol @@ -702,7 +702,7 @@ void Radiation::get_gas_vmr(const std::vector& gas_names, real3d& g // first specific humidity (held in the mass_mix_ratio array read // from rad_constituents) is converted to an actual mass mixing ratio. yakl::c::parallel_for(yakl::c::Bounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) { - gas_vmr(igas,icol,ilev) = mmr(icol,ilev) / ( + gas_vmr(igas,icol,ilev) = mmr(icol,ilev) / ( 1. - mmr(icol,ilev))*mol_weight_air / mol_weight_gas[igas]; }); } else { @@ -712,7 +712,7 @@ void Radiation::get_gas_vmr(const std::vector& gas_names, real3d& g // Convert to volume mixing ratio by multiplying by the ratio of // molecular weight of dry air to molecular weight of gas yakl::c::parallel_for(yakl::c::Bounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) { - gas_vmr(igas,icol,ilev) = mmr(icol,ilev) + gas_vmr(igas,icol,ilev) = mmr(icol,ilev) * mol_weight_air / mol_weight_gas[igas]; }); } @@ -730,7 +730,7 @@ void Radiation::get_gas_vmr(const std::vector& gas_names, real3d& g // Why? Something to do with convenience with applying the fluxes to the // heating tendency? void Radiation::calculate_heating_rate(real2d& flux_up, real2d& flux_dn, real2d& pint, real2d& heating_rate) { - for (auto ilev = 0; ilev < nlev; ++ilev) { + for (auto ilev = 0; ilev < nlev; ++ilev) { for (auto icol = 0; icol < ncol; ++icol) { heating_rate(icol,ilev) = (flux_up(icol,ilev+1) - flux_up(icol,ilev)- flux_dn(icol,ilev+1)+flux_dn(icol,ilev)) *CONST_GRAV/(pint(icol,ilev+1)-pint(icol,ilev)); @@ -738,6 +738,6 @@ void Radiation::calculate_heating_rate(real2d& flux_up, real2d& flux_dn, real2d& } } -// call back +// call back void Radiation::on_complete() { }