ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
rrtmgp Namespace Reference

Functions

optical_props2_t get_cloud_optics_sw (const int ncol, const int nlay, cloud_optics_t &cloud_optics, gas_optics_t &kdist, real2d_k &lwp, real2d_k &iwp, real2d_k &rel, real2d_k &rei)
 
optical_props1_t get_cloud_optics_lw (const int ncol, const int nlay, cloud_optics_t &cloud_optics, gas_optics_t &kdist, real2d_k &lwp, real2d_k &iwp, real2d_k &rel, real2d_k &rei)
 
optical_props2_t get_subsampled_clouds (const int ncol, const int nlay, const int nbnd, const int ngpt, optical_props2_t &cloud_optics, gas_optics_t &kdist, real2d_k &cld, real2d_k &p_lay)
 
optical_props1_t get_subsampled_clouds (const int ncol, const int nlay, const int nbnd, const int ngpt, optical_props1_t &cloud_optics, gas_optics_t &kdist, real2d_k &cld, real2d_k &p_lay)
 
void rrtmgp_initialize (gas_concs_t &gas_concs_k, const std::string &coefficients_file_sw, const std::string &coefficients_file_lw, const std::string &cloud_optics_file_sw, const std::string &cloud_optics_file_lw)
 
void rrtmgp_finalize ()
 
void compute_band_by_band_surface_albedos (const int ncol, const int nswbands, real1d_k &sfc_alb_dir_vis, real1d_k &sfc_alb_dir_nir, real1d_k &sfc_alb_dif_vis, real1d_k &sfc_alb_dif_nir, real2d_k &sfc_alb_dir, real2d_k &sfc_alb_dif)
 
void compute_broadband_surface_fluxes (const int ncol, const int kbot, const int nswbands, real3d_k &sw_bnd_flux_dir, real3d_k &sw_bnd_flux_dif, real1d_k &sfc_flux_dir_vis, real1d_k &sfc_flux_dir_nir, real1d_k &sfc_flux_dif_vis, real1d_k &sfc_flux_dif_nir)
 
void rrtmgp_main (const int ncol, const int nlay, real2d_k &p_lay, real2d_k &t_lay, real2d_k &p_lev, real2d_k &t_lev, gas_concs_t &gas_concs, real2d_k &sfc_alb_dir, real2d_k &sfc_alb_dif, real1d_k &mu0, real1d_k &t_sfc, real2d_k &emis_sfc, real1d_k &lw_src, real2d_k &lwp, real2d_k &iwp, real2d_k &rel, real2d_k &rei, real2d_k &cldfrac, real3d_k &aer_tau_sw, real3d_k &aer_ssa_sw, real3d_k &aer_asm_sw, real3d_k &aer_tau_lw, real3d_k &cld_tau_sw_bnd, real3d_k &cld_tau_lw_bnd, real3d_k &cld_tau_sw_gpt, real3d_k &cld_tau_lw_gpt, real2d_k &sw_flux_up, real2d_k &sw_flux_dn, real2d_k &sw_flux_dn_dir, real2d_k &lw_flux_up, real2d_k &lw_flux_dn, real2d_k &sw_clnclrsky_flux_up, real2d_k &sw_clnclrsky_flux_dn, real2d_k &sw_clnclrsky_flux_dn_dir, real2d_k &sw_clrsky_flux_up, real2d_k &sw_clrsky_flux_dn, real2d_k &sw_clrsky_flux_dn_dir, real2d_k &sw_clnsky_flux_up, real2d_k &sw_clnsky_flux_dn, real2d_k &sw_clnsky_flux_dn_dir, real2d_k &lw_clnclrsky_flux_up, real2d_k &lw_clnclrsky_flux_dn, real2d_k &lw_clrsky_flux_up, real2d_k &lw_clrsky_flux_dn, real2d_k &lw_clnsky_flux_up, real2d_k &lw_clnsky_flux_dn, real3d_k &sw_bnd_flux_up, real3d_k &sw_bnd_flux_dn, real3d_k &sw_bnd_flux_dn_dir, real3d_k &lw_bnd_flux_up, real3d_k &lw_bnd_flux_dn, const RealT tsi_scaling, const bool extra_clnclrsky_diag, const bool extra_clnsky_diag)
 
int3d_k get_subcolumn_mask (const int ncol, const int nlay, const int ngpt, real2d_k &cldf, const int overlap_option, int1d_k &seeds)
 
void rrtmgp_sw (const int ncol, const int nlay, gas_optics_t &k_dist, real2d_k &p_lay, real2d_k &t_lay, real2d_k &p_lev, real2d_k &t_lev, gas_concs_t &gas_concs, real2d_k &sfc_alb_dir, real2d_k &sfc_alb_dif, real1d_k &mu0, optical_props2_t &aerosol, optical_props2_t &clouds, fluxes_t &fluxes, fluxes_broadband_t &clnclrsky_fluxes, fluxes_broadband_t &clrsky_fluxes, fluxes_broadband_t &clnsky_fluxes, const RealT tsi_scaling, const bool extra_clnclrsky_diag, const bool extra_clnsky_diag)
 
void rrtmgp_lw (const int ncol, const int nlay, gas_optics_t &k_dist, real2d_k &p_lay, real2d_k &t_lay, real2d_k &p_lev, real2d_k &t_lev, real1d_k &t_sfc, real2d_k &emis_sfc, real1d_k &lw_src, gas_concs_t &gas_concs, optical_props1_t &aerosol, optical_props1_t &clouds, fluxes_t &fluxes, fluxes_broadband_t &clnclrsky_fluxes, fluxes_broadband_t &clrsky_fluxes, fluxes_broadband_t &clnsky_fluxes, const bool extra_clnclrsky_diag, const bool extra_clnsky_diag)
 
void compute_cloud_area (int ncol, int nlay, int ngpt, const RealT pmin, const RealT pmax, const real2d_k &pmid, const real3d_k &cld_tau_gpt, real1d_k &cld_area)
 
int get_wavelength_index_sw (RealT wavelength)
 
int get_wavelength_index_lw (RealT wavelength)
 
int get_wavelength_index (optical_props_t &kdist, RealT wavelength)
 
void compute_aerocom_cloudtop (int ncol, int nlay, const real2d_k &tmid, const real2d_k &pmid, const real2d_k &p_del, const real2d_k &z_del, const real2d_k &qc, const real2d_k &qi, const real2d_k &rel, const real2d_k &rei, const real2d_k &cldfrac_tot, const real2d_k &nc, real1d_k &T_mid_at_cldtop, real1d_k &p_mid_at_cldtop, real1d_k &cldfrac_ice_at_cldtop, real1d_k &cldfrac_liq_at_cldtop, real1d_k &cldfrac_tot_at_cldtop, real1d_k &cdnc_at_cldtop, real1d_k &eff_radius_qc_at_cldtop, real1d_k &eff_radius_qi_at_cldtop)
 
template<class View1 , class View2 , class View3 , class View4 , class View5 >
void mixing_ratio_to_cloud_mass (View1 const &mixing_ratio, View2 const &cloud_fraction, View3 const &rho, View4 const &dz, View5 const &cloud_mass)
 
template<typename InT , typename OutT , typename T >
void limit_to_bounds_1d (InT const &arr_in, T const lower, T const upper, OutT &arr_out)
 
template<typename InT , typename OutT , typename T >
void limit_to_bounds_2d (InT const &arr_in, T const lower, T const upper, OutT &arr_out)
 
template<class View1 , class View2 , class View3 , class View4 , class View5 >
void compute_heating_rate (View1 const &flux_up, View2 const &flux_dn, View3 const &rho, View4 const &dz, View5 &heating_rate)
 
bool radiation_do (const int irad, const int nstep)
 

Variables

std::unique_ptr< gas_optics_tk_dist_sw_k
 
std::unique_ptr< gas_optics_tk_dist_lw_k
 
std::unique_ptr< cloud_optics_tcloud_optics_sw_k
 
std::unique_ptr< cloud_optics_tcloud_optics_lw_k
 
pool_t kokkos_mem_pool
 
bool initialized = false
 

Function Documentation

◆ compute_aerocom_cloudtop()

void rrtmgp::compute_aerocom_cloudtop ( int  ncol,
int  nlay,
const real2d_k tmid,
const real2d_k pmid,
const real2d_k p_del,
const real2d_k z_del,
const real2d_k qc,
const real2d_k qi,
const real2d_k rel,
const real2d_k rei,
const real2d_k cldfrac_tot,
const real2d_k nc,
real1d_k T_mid_at_cldtop,
real1d_k p_mid_at_cldtop,
real1d_k cldfrac_ice_at_cldtop,
real1d_k cldfrac_liq_at_cldtop,
real1d_k cldfrac_tot_at_cldtop,
real1d_k cdnc_at_cldtop,
real1d_k eff_radius_qc_at_cldtop,
real1d_k eff_radius_qi_at_cldtop 
)
1141 {
1142  /* The goal of this routine is to calculate properties at cloud top
1143  * based on the AeroCom recommendation. See reference for routine
1144  * get_subcolumn_mask above, where equation 14 is used for the
1145  * maximum-random overlap assumption for subcolumn generation. We use
1146  * equation 13, the column counterpart.
1147  */
1148  // Set outputs to zero
1149  Kokkos::deep_copy(T_mid_at_cldtop, 0.0);
1150  Kokkos::deep_copy(p_mid_at_cldtop, 0.0);
1151  Kokkos::deep_copy(cldfrac_ice_at_cldtop, 0.0);
1152  Kokkos::deep_copy(cldfrac_liq_at_cldtop, 0.0);
1153  Kokkos::deep_copy(cldfrac_tot_at_cldtop, 0.0);
1154  Kokkos::deep_copy(cdnc_at_cldtop, 0.0);
1155  Kokkos::deep_copy(eff_radius_qc_at_cldtop, 0.0);
1156  Kokkos::deep_copy(eff_radius_qi_at_cldtop, 0.0);
1157 
1158  // Initialize the 1D "clear fraction" as 1 (totally clear)
1159  real1d_k aerocom_clr("aerocom_clr", ncol);
1160  Kokkos::deep_copy(aerocom_clr, 1.0);
1161 
1162  // TODO: move tunable constant to namelist
1163  constexpr RealT q_threshold = 0.0; // BAD_CONSTANT!
1164 
1165  // TODO: move tunable constant to namelist
1166  constexpr RealT cldfrac_tot_threshold = 0.001; // BAD_CONSTANT!
1167 
1168  // Loop over all columns in parallel
1169  Kokkos::parallel_for(ncol, KOKKOS_LAMBDA(int icol)
1170  {
1171  // Loop over all layers in serial (due to accumulative
1172  // product), starting at 2 (second highest) layer because the
1173  // highest is assumed to have no clouds
1174  for(int ilay = 1; ilay < nlay; ++ilay) {
1175  // Only do the calculation if certain conditions are met
1176  if((qc(icol, ilay) + qi(icol, ilay)) > q_threshold &&
1177  (cldfrac_tot(icol, ilay) > cldfrac_tot_threshold)) {
1178  /* PART I: Probabilistically determining cloud top */
1179  // Populate aerocom_tmp as the clear-sky fraction
1180  // probability of this level, where aerocom_clr is that of
1181  // the previous level
1182  auto aerocom_tmp = aerocom_clr(icol) *
1183  (1.0 - std::max(cldfrac_tot(icol, ilay - 1),
1184  cldfrac_tot(icol, ilay))) /
1185  (1.0 - std::min(cldfrac_tot(icol, ilay - 1),
1186  1.0 - cldfrac_tot_threshold));
1187  // Temporary variable for probability "weights"
1188  auto aerocom_wts = aerocom_clr(icol) - aerocom_tmp;
1189  // Temporary variable for liquid "phase"
1190  auto aerocom_phi = qc(icol, ilay) / (qc(icol, ilay) + qi(icol, ilay));
1191  /* PART II: The inferred properties */
1192  /* In general, converting a 3D property X to a 2D cloud-top
1193  * counterpart x follows: x(i) += X(i,k) * weights * Phase
1194  * but X and Phase are not always needed */
1195  // T_mid_at_cldtop
1196  T_mid_at_cldtop(icol) += tmid(icol, ilay) * aerocom_wts;
1197  // p_mid_at_cldtop
1198  p_mid_at_cldtop(icol) += pmid(icol, ilay) * aerocom_wts;
1199  // cldfrac_ice_at_cldtop
1200  cldfrac_ice_at_cldtop(icol) += (1.0 - aerocom_phi) * aerocom_wts;
1201  // cldfrac_liq_at_cldtop
1202  cldfrac_liq_at_cldtop(icol) += aerocom_phi * aerocom_wts;
1203  // cdnc_at_cldtop
1204  /* We need to convert nc from 1/mass to 1/volume first, and
1205  * from grid-mean to in-cloud, but after that, the
1206  * calculation follows the general logic */
1207  // AML NOTE: p_del/z_del/g should be replaced with RHO for our dycore
1208  auto cdnc = nc(icol, ilay) * p_del(icol, ilay) /
1209  z_del(icol, ilay) / CONST_GRAV /
1210  cldfrac_tot(icol, ilay);
1211  cdnc_at_cldtop(icol) += cdnc * aerocom_phi * aerocom_wts;
1212  // eff_radius_qc_at_cldtop
1213  eff_radius_qc_at_cldtop(icol) += rel(icol, ilay) * aerocom_phi * aerocom_wts;
1214  // eff_radius_qi_at_cldtop
1215  eff_radius_qi_at_cldtop(icol) += rei(icol, ilay) * (1.0 - aerocom_phi) * aerocom_wts;
1216  // Reset aerocom_clr to aerocom_tmp to accumulate
1217  aerocom_clr(icol) = aerocom_tmp;
1218  }
1219  }
1220  // After the serial loop over levels, the cloudy fraction is
1221  // defined as (1 - aerocom_clr). This is true because
1222  // aerocom_clr is the result of accumulative probabilities
1223  // (their products)
1224  cldfrac_tot_at_cldtop(icol) = 1.0 - aerocom_clr(icol);
1225  });
1226 }
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
Kokkos::View< RealT *, KokkosDefaultDevice > real1d_k
Definition: ERF_Kokkos.H:17
amrex::Real RealT
Definition: ERF_Kokkos.H:14
@ nc
Definition: ERF_Morrison.H:44
@ qc
Definition: ERF_SatAdj.H:36

◆ compute_band_by_band_surface_albedos()

void rrtmgp::compute_band_by_band_surface_albedos ( const int  ncol,
const int  nswbands,
real1d_k sfc_alb_dir_vis,
real1d_k sfc_alb_dir_nir,
real1d_k sfc_alb_dif_vis,
real1d_k sfc_alb_dif_nir,
real2d_k sfc_alb_dir,
real2d_k sfc_alb_dif 
)
290 {
291  auto wavenumber_limits = k_dist_sw_k->get_band_lims_wavenumber();
292 
293  // Loop over bands, and determine for each band whether it is broadly in the
294  // visible or infrared part of the spectrum (visible or "not visible")
295  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, nswbands}),
296  KOKKOS_LAMBDA (int icol, int ibnd)
297  {
298  // Threshold between visible and infrared is 0.7 micron, or 14286 cm^-1.
299  const RealT visible_wavenumber_threshold = 14286.0;
300 
301  // Wavenumber is in the visible if it is above the visible wavenumber
302  // threshold, and in the infrared if it is below the threshold
303  const bool is_visible_wave1 = (wavenumber_limits(0, ibnd) > visible_wavenumber_threshold ? true : false);
304  const bool is_visible_wave2 = (wavenumber_limits(1, ibnd) > visible_wavenumber_threshold ? true : false);
305 
306  if (is_visible_wave1 && is_visible_wave2) {
307  // Entire band is in the visible
308  sfc_alb_dir(icol,ibnd) = sfc_alb_dir_vis(icol);
309  sfc_alb_dif(icol,ibnd) = sfc_alb_dif_vis(icol);
310  } else if (!is_visible_wave1 && !is_visible_wave2) {
311  // Entire band is in the longwave (near-infrared)
312  sfc_alb_dir(icol,ibnd) = sfc_alb_dir_nir(icol);
313  sfc_alb_dif(icol,ibnd) = sfc_alb_dif_nir(icol);
314  } else {
315  // Band straddles the visible to near-infrared transition, so we take
316  // the albedo to be the average of the visible and near-infrared
317  // broadband albedos
318  sfc_alb_dir(icol,ibnd) = 0.5*(sfc_alb_dir_vis(icol) + sfc_alb_dir_nir(icol));
319  sfc_alb_dif(icol,ibnd) = 0.5*(sfc_alb_dif_vis(icol) + sfc_alb_dif_nir(icol));
320  }
321  });
322 }
@ sfc_alb_dif_vis
Definition: ERF_NOAHMP.H:29
@ sfc_alb_dir_vis
Definition: ERF_NOAHMP.H:27
@ sfc_alb_dir_nir
Definition: ERF_NOAHMP.H:28
@ sfc_alb_dif_nir
Definition: ERF_NOAHMP.H:30
std::unique_ptr< gas_optics_t > k_dist_sw_k
Definition: ERF_RRTMGP_Interface.cpp:10

Referenced by Radiation::run_impl().

Here is the caller graph for this function:

◆ compute_broadband_surface_fluxes()

void rrtmgp::compute_broadband_surface_fluxes ( const int  ncol,
const int  kbot,
const int  nswbands,
real3d_k sw_bnd_flux_dir,
real3d_k sw_bnd_flux_dif,
real1d_k sfc_flux_dir_vis,
real1d_k sfc_flux_dir_nir,
real1d_k sfc_flux_dif_vis,
real1d_k sfc_flux_dif_nir 
)
335 {
336  // Band 10 straddles the near-IR and visible, so divide contributions from band 10 between both broadband sums
337  // TODO: Hard-coding these band indices is really bad practice. If the bands ever were to change (like when
338  // the RRTMG bands were re-ordered for RRTMGP), we would be using the wrong bands for the IR and UV/VIS. This
339  // should be refactored to grab the correct bands by specifying appropriate wavenumber rather than index.
340  //sfc_flux_dir_nir(i) = sum(sw_bnd_flux_dir(i+1,kbot,1:9)) + 0.5 * sw_bnd_flux_dir(i+1,kbot,10);
341  //sfc_flux_dir_vis(i) = sum(sw_bnd_flux_dir(i+1,kbot,11:14)) + 0.5 * sw_bnd_flux_dir(i+1,kbot,10);
342  //sfc_flux_dif_nir(i) = sum(sw_bnd_flux_dif(i+1,kbot,1:9)) + 0.5 * sw_bnd_flux_dif(i+1,kbot,10);
343  //sfc_flux_dif_vis(i) = sum(sw_bnd_flux_dif(i+1,kbot,11:14)) + 0.5 * sw_bnd_flux_dif(i+1,kbot,10);
344 
345  // Initialize sums over bands
346  Kokkos::deep_copy(sfc_flux_dir_nir, 0);
347  Kokkos::deep_copy(sfc_flux_dir_vis, 0);
348  Kokkos::deep_copy(sfc_flux_dif_nir, 0);
349  Kokkos::deep_copy(sfc_flux_dif_vis, 0);
350 
351  // Threshold between visible and infrared is 0.7 micron, or 14286 cm^-1.
352  const RealT visible_wavenumber_threshold = 14286.0;
353  auto wavenumber_limits = k_dist_sw_k->get_band_lims_wavenumber();
354  Kokkos::parallel_for(ncol, KOKKOS_LAMBDA(const int icol)
355  {
356  for (int ibnd = 0; ibnd < nswbands; ++ibnd) {
357  // Wavenumber is in the visible if it is above the visible wavenumber
358  // threshold, and in the infrared if it is below the threshold
359  const bool is_visible_wave1 = (wavenumber_limits(0, ibnd) > visible_wavenumber_threshold ? true : false);
360  const bool is_visible_wave2 = (wavenumber_limits(1, ibnd) > visible_wavenumber_threshold ? true : false);
361 
362  if (is_visible_wave1 && is_visible_wave2) {
363  // Entire band is in the visible
364  sfc_flux_dir_vis(icol) += sw_bnd_flux_dir(icol,kbot,ibnd);
365  sfc_flux_dif_vis(icol) += sw_bnd_flux_dif(icol,kbot,ibnd);
366  } else if (!is_visible_wave1 && !is_visible_wave2) {
367  // Entire band is in the longwave (near-infrared)
368  sfc_flux_dir_nir(icol) += sw_bnd_flux_dir(icol,kbot,ibnd);
369  sfc_flux_dif_nir(icol) += sw_bnd_flux_dif(icol,kbot,ibnd);
370  } else {
371  // Band straddles the visible to near-infrared transition, so put half
372  // the flux in visible and half in near-infrared fluxes
373  sfc_flux_dir_vis(icol) += 0.5 * sw_bnd_flux_dir(icol,kbot,ibnd);
374  sfc_flux_dif_vis(icol) += 0.5 * sw_bnd_flux_dif(icol,kbot,ibnd);
375  sfc_flux_dir_nir(icol) += 0.5 * sw_bnd_flux_dir(icol,kbot,ibnd);
376  sfc_flux_dif_nir(icol) += 0.5 * sw_bnd_flux_dif(icol,kbot,ibnd);
377  }
378  }
379  });
380 }

Referenced by Radiation::run_impl().

Here is the caller graph for this function:

◆ compute_cloud_area()

void rrtmgp::compute_cloud_area ( int  ncol,
int  nlay,
int  ngpt,
const RealT  pmin,
const RealT  pmax,
const real2d_k pmid,
const real3d_k cld_tau_gpt,
real1d_k cld_area 
)
1065 {
1066  // Subcolumn binary cld mask; if any layers with pressure between pmin and pmax are cloudy
1067  // then 2d subcol mask is 1, otherwise it is 0
1068  real2d_k subcol_mask("subcol_mask", ncol, ngpt);
1069  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, ngpt}),
1070  KOKKOS_LAMBDA (int icol, int ilay, int igpt)
1071  {
1072  // NOTE: using plev would need to assume level ordering (top to bottom or bottom to top), but
1073  // using play/pmid does not
1074  if (cld_tau_gpt(icol,ilay,igpt) > 0 && pmid(icol,ilay) >= pmin && pmid(icol,ilay) < pmax) {
1075  subcol_mask(icol,igpt) = 1.0;
1076  } else {
1077  subcol_mask(icol,igpt) = 0.0;
1078  }
1079  });
1080  // Compute average over subcols to get cloud area
1081  auto ngpt_inv = 1.0 / ngpt;
1082  Kokkos::deep_copy(cld_area, 0);
1083  Kokkos::parallel_for(ncol, KOKKOS_LAMBDA(int icol)
1084  {
1085  // This loop needs to be serial because of the atomic reduction
1086  for (int igpt = 0; igpt < ngpt; ++igpt) {
1087  cld_area(icol) += subcol_mask(icol,igpt) * ngpt_inv;
1088  }
1089  });
1090 }
Kokkos::View< RealT **, layout_t, KokkosDefaultDevice > real2d_k
Definition: ERF_Kokkos.H:18

◆ compute_heating_rate()

template<class View1 , class View2 , class View3 , class View4 , class View5 >
void rrtmgp::compute_heating_rate ( View1 const &  flux_up,
View2 const &  flux_dn,
View3 const &  rho,
View4 const &  dz,
View5 &  heating_rate 
)
86 {
87  const int ncol = (int)flux_up.extent(0);
88  const int nlay = (int)flux_up.extent(1)-1;
89  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, nlay}),
90  KOKKOS_LAMBDA (int icol, int ilay)
91  {
92  heating_rate(icol,ilay) = ( flux_up(icol,ilay+1) - flux_up(icol,ilay)
93  - flux_dn(icol,ilay+1) + flux_dn(icol,ilay) )
94  / ( -rho(icol,ilay) * dz(icol,ilay) * Cp_d );
95  });
96 }
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:12
@ rho
Definition: ERF_Kessler.H:22

Referenced by Radiation::finalize_impl(), and Radiation::run_impl().

Here is the caller graph for this function:

◆ get_cloud_optics_lw()

optical_props1_t rrtmgp::get_cloud_optics_lw ( const int  ncol,
const int  nlay,
cloud_optics_t cloud_optics,
gas_optics_t kdist,
real2d_k lwp,
real2d_k iwp,
real2d_k rel,
real2d_k rei 
)
70 {
71  // Initialize optics
72  optical_props1_t clouds;
73  clouds.init(kdist.get_band_lims_wavenumber());
74  clouds.alloc_1scl(ncol, nlay);
75 
76  // Needed for consistency with all-sky example problem?
77  cloud_optics.set_ice_roughness(2);
78 
79  // Limit effective radii to be within bounds of lookup table
80  real2d_k rel_limited("rel_limited", ncol, nlay);
81  real2d_k rei_limited("rei_limited", ncol, nlay);
82  limit_to_bounds_2d(rel, cloud_optics.radliq_lwr,
83  cloud_optics.radliq_upr, rel_limited);
84  limit_to_bounds_2d(rei, cloud_optics.radice_lwr,
85  cloud_optics.radice_upr, rei_limited);
86 
87  // Calculate cloud optics
88  cloud_optics.cloud_optics(ncol, nlay, lwp, iwp, rel_limited, rei_limited, clouds);
89 
90  // Return optics
91  return clouds;
92 }
OpticalProps1sclK< RealT, layout_t, KokkosDefaultDevice > optical_props1_t
Definition: ERF_RRTMGP_Interface.H:36
void limit_to_bounds_2d(InT const &arr_in, T const lower, T const upper, OutT &arr_out)
Definition: ERF_RRTMGP_Utils.H:55

Referenced by rrtmgp_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_cloud_optics_sw()

optical_props2_t rrtmgp::get_cloud_optics_sw ( const int  ncol,
const int  nlay,
cloud_optics_t cloud_optics,
gas_optics_t kdist,
real2d_k lwp,
real2d_k iwp,
real2d_k rel,
real2d_k rei 
)
36 {
37  // Initialize optics
38  optical_props2_t clouds;
39  clouds.init(kdist.get_band_lims_wavenumber());
40  clouds.alloc_2str(ncol, nlay);
41 
42  // Needed for consistency with all-sky example problem?
43  cloud_optics.set_ice_roughness(2);
44 
45  // Limit effective radii to be within bounds of lookup table
46  real2d_k rel_limited("rel_limited", ncol, nlay);
47  real2d_k rei_limited("rei_limited", ncol, nlay);
48  limit_to_bounds_2d(rel, cloud_optics.radliq_lwr,
49  cloud_optics.radliq_upr, rel_limited);
50  limit_to_bounds_2d(rei, cloud_optics.radice_lwr,
51  cloud_optics.radice_upr, rei_limited);
52 
53  // Calculate cloud optics
54  cloud_optics.cloud_optics(ncol, nlay, lwp, iwp, rel_limited, rei_limited, clouds);
55 
56  // Return optics
57  return clouds;
58 }
OpticalProps2strK< RealT, layout_t, KokkosDefaultDevice > optical_props2_t
Definition: ERF_RRTMGP_Interface.H:37

Referenced by rrtmgp_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_subcolumn_mask()

int3d_k rrtmgp::get_subcolumn_mask ( const int  ncol,
const int  nlay,
const int  ngpt,
real2d_k cldf,
const int  overlap_option,
int1d_k seeds 
)
534 {
535  // Routine will return subcolumn mask with values of 0 indicating no cloud, 1 indicating cloud
536  int3d_k subcolumn_mask("subcolumn_mask", ncol, nlay, ngpt);
537 
538  // Subcolumn generators are a means for producing a variable x(i,j,k), where
539  //
540  // c(i,j,k) = 1 for x(i,j,k) > 1 - cldf(i,j)
541  // c(i,j,k) = 0 for x(i,j,k) <= 1 - cldf(i,j)
542  //
543  // I am going to call this "cldx" to be just slightly less ambiguous
544  real3d_k cldx("cldx", ncol, nlay, ngpt);
545 
546  // Apply overlap assumption to set cldx
547  if (overlap_option == 0) { // Dummy mask, always cloudy
548  Kokkos::deep_copy(cldx, 1);
549  } else { // Default case, maximum-random overlap
550  // Maximum-random overlap:
551  // Uses essentially the algorithm described in eq (14) in Raisanen et al. 2004,
552  // https://rmets.onlinelibrary.wiley.com/doi/epdf/10.1256/qj.03.99. Also the same
553  // algorithm used in RRTMG implementation of maximum-random overlap (see
554  // https://github.com/AER-RC/RRTMG_SW/blob/master/src/mcica_subcol_gen_sw.f90)
555  //
556  // First, fill cldx with random numbers.
557  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, ngpt}),
558  KOKKOS_LAMBDA (int icol, int ilay, int igpt)
559  {
560  conv::Random rand(seeds(icol) + ilay*ngpt + igpt);
561  cldx(icol,ilay,igpt) = rand.genFP<RealT>();
562  });
563 
564  // Step down columns and apply algorithm from eq (14)
565  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, ngpt}),
566  KOKKOS_LAMBDA (int icol, int igpt)
567  {
568  for (int ilay = 1; ilay < nlay; ilay++) {
569  // Check cldx in level above and see if it satisfies conditions to create a cloudy subcolumn
570  if (cldx(icol,ilay-1,igpt) > 1.0 - cldf(icol,ilay-1)) {
571  // Cloudy subcolumn above, use same random number here so that clouds in these two adjacent
572  // layers are maximimally overlapped
573  cldx(icol,ilay,igpt) = cldx(icol,ilay-1,igpt);
574  } else {
575  // Cloud-less above, use new random number so that clouds are distributed
576  // randomly in this layer. Need to scale new random number to range
577  // [0, 1.0 - cldf(ilay-1)] because we have artificially changed the distribution
578  // of random numbers in this layer with the above branch of the conditional,
579  // which would otherwise inflate cloud fraction in this layer.
580  cldx(icol,ilay,igpt) = cldx(icol,ilay ,igpt) * (1.0 - cldf(icol,ilay-1));
581  }
582  }
583  });
584  }
585 
586  // Use cldx array to create subcolumn mask
587  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, ngpt}),
588  KOKKOS_LAMBDA (int icol, int ilay, int igpt)
589  {
590  if (cldx(icol,ilay,igpt) > 1.0 - cldf(icol,ilay)) {
591  subcolumn_mask(icol,ilay,igpt) = 1;
592  } else {
593  subcolumn_mask(icol,ilay,igpt) = 0;
594  }
595  });
596  return subcolumn_mask;
597 }
Kokkos::View< RealT ***, layout_t, KokkosDefaultDevice > real3d_k
Definition: ERF_Kokkos.H:19
Kokkos::View< int ***, layout_t, KokkosDefaultDevice > int3d_k
Definition: ERF_Kokkos.H:22

Referenced by get_subsampled_clouds().

Here is the caller graph for this function:

◆ get_subsampled_clouds() [1/2]

optical_props1_t rrtmgp::get_subsampled_clouds ( const int  ncol,
const int  nlay,
const int  nbnd,
const int  ngpt,
optical_props1_t cloud_optics,
gas_optics_t kdist,
real2d_k cld,
real2d_k p_lay 
)
172 {
173  // Initialized subsampled optics
174  optical_props1_t subsampled_optics;
175  subsampled_optics.init(kdist.get_band_lims_wavenumber(), kdist.get_band_lims_gpoint(), "subsampled_optics");
176  subsampled_optics.alloc_1scl(ncol, nlay);
177 
178  // Check that we do not have clouds with no optical properties; this would get corrected
179  // when we assign optical props, but we want to use a "radiative cloud fraction"
180  // for the subcolumn sampling too because otherwise we can get vertically-contiguous cloud
181  // mask profiles with no actual cloud properties in between, which would just further overestimate
182  // the vertical correlation of cloudy layers. I.e., cloudy layers might look maximally overlapped
183  // even when separated by layers with no cloud properties, when in fact those layers should be
184  // randomly overlapped.
185  real2d_k cldfrac_rad("cldfrac_rad", ncol, nlay);
186  Kokkos::deep_copy(cldfrac_rad, 0.0);
187  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, nbnd}),
188  KOKKOS_LAMBDA (int icol, int ilay, int ibnd)
189  {
190  if (cloud_optics.tau(icol,ilay,ibnd) > 0) {
191  cldfrac_rad(icol,ilay) = cld(icol,ilay);
192  }
193  });
194 
195  // Get subcolumn cloud mask
196  int overlap = 1;
197  // Get unique seeds for each column that are reproducible across different MPI rank layouts;
198  // use decimal part of pressure for this, consistent with the implementation in EAM; use different
199  // seed values for longwave and shortwave
200  int1d_k seeds("seeds", ncol);
201  Kokkos::parallel_for(ncol, KOKKOS_LAMBDA(int icol)
202  {
203  seeds(icol) = 1e9 * (p_lay(icol,nlay-1) - int(p_lay(icol,nlay-1)));
204  });
205  auto cldmask = get_subcolumn_mask(ncol, nlay, ngpt, cldfrac_rad, overlap, seeds);
206 
207  // Assign optical properties to subcolumns (note this implements MCICA)
208  auto gpoint_bands = kdist.get_gpoint_bands();
209  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, ngpt}),
210  KOKKOS_LAMBDA (int icol, int ilay, int igpt)
211  {
212  auto ibnd = gpoint_bands(igpt);
213  if (cldmask(icol,ilay,igpt) == 1) {
214  subsampled_optics.tau(icol,ilay,igpt) = cloud_optics.tau(icol,ilay,ibnd);
215  } else {
216  subsampled_optics.tau(icol,ilay,igpt) = 0;
217  }
218  });
219  return subsampled_optics;
220 }
Kokkos::View< int *, KokkosDefaultDevice > int1d_k
Definition: ERF_Kokkos.H:20
int3d_k get_subcolumn_mask(const int ncol, const int nlay, const int ngpt, real2d_k &cldf, const int overlap_option, int1d_k &seeds)
Definition: ERF_RRTMGP_Interface.cpp:528
Here is the call graph for this function:

◆ get_subsampled_clouds() [2/2]

optical_props2_t rrtmgp::get_subsampled_clouds ( const int  ncol,
const int  nlay,
const int  nbnd,
const int  ngpt,
optical_props2_t cloud_optics,
gas_optics_t kdist,
real2d_k cld,
real2d_k p_lay 
)
104 {
105  // Initialized subsampled optics
106  optical_props2_t subsampled_optics;
107  subsampled_optics.init(kdist.get_band_lims_wavenumber(), kdist.get_band_lims_gpoint(), "subsampled_optics");
108  subsampled_optics.alloc_2str(ncol, nlay);
109 
110  // Check that we do not have clouds with no optical properties; this would get corrected
111  // when we assign optical props, but we want to use a "radiative cloud fraction"
112  // for the subcolumn sampling too because otherwise we can get vertically-contiguous cloud
113  // mask profiles with no actual cloud properties in between, which would just further overestimate
114  // the vertical correlation of cloudy layers. I.e., cloudy layers might look maximally overlapped
115  // even when separated by layers with no cloud properties, when in fact those layers should be
116  // randomly overlapped.
117  real2d_k cldfrac_rad("cldfrac_rad", ncol, nlay);
118  Kokkos::deep_copy(cldfrac_rad, 0.0);
119  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, nbnd}),
120  KOKKOS_LAMBDA (int icol, int ilay, int ibnd)
121  {
122  if (cloud_optics.tau(icol,ilay,ibnd) > 0) {
123  cldfrac_rad(icol,ilay) = cld(icol,ilay);
124  }
125  });
126 
127  // Get subcolumn cloud mask; note that get_subcolumn_mask exposes overlap assumption as an option,
128  // but the only currently supported options are 0 (trivial all-or-nothing cloud) or 1 (max-rand),
129  // so overlap has not been exposed as an option beyond this subcolumn. In the future, we should
130  // support generalized overlap as well, with parameters derived from DPSCREAM simulations with very
131  // high resolution.
132  int overlap = 1;
133 
134  // Get unique seeds for each column that are reproducible across different MPI rank layouts;
135  // use decimal part of pressure for this, consistent with the implementation in EAM
136  int1d_k seeds("seeds", ncol);
137  Kokkos::parallel_for(ncol, KOKKOS_LAMBDA(int icol)
138  {
139  seeds(icol) = 1.0e9 * (p_lay(icol,nlay-1) - int(p_lay(icol,nlay-1)));
140  });
141  auto cldmask = get_subcolumn_mask(ncol, nlay, ngpt, cldfrac_rad, overlap, seeds);
142 
143  // Assign optical properties to subcolumns (note this implements MCICA)
144  auto gpoint_bands = kdist.get_gpoint_bands();
145  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, ngpt}),
146  KOKKOS_LAMBDA (int icol, int ilay, int igpt)
147  {
148  auto ibnd = gpoint_bands(igpt);
149  if (cldmask(icol,ilay,igpt) == 1) {
150  subsampled_optics.tau(icol,ilay,igpt) = cloud_optics.tau(icol,ilay,ibnd);
151  subsampled_optics.ssa(icol,ilay,igpt) = cloud_optics.ssa(icol,ilay,ibnd);
152  subsampled_optics.g (icol,ilay,igpt) = cloud_optics.g (icol,ilay,ibnd);
153  } else {
154  subsampled_optics.tau(icol,ilay,igpt) = 0;
155  subsampled_optics.ssa(icol,ilay,igpt) = 0;
156  subsampled_optics.g (icol,ilay,igpt) = 0;
157  }
158  });
159  return subsampled_optics;
160 }

Referenced by rrtmgp_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_wavelength_index()

int rrtmgp::get_wavelength_index ( optical_props_t kdist,
RealT  wavelength 
)
1104 {
1105  // Get wavelength bounds for all wavelength bands
1106  auto band_lims_wvn = kdist.get_band_lims_wavenumber();
1107  real2d_k wavelength_bounds("wavelength_bounds",band_lims_wvn.extent(0), band_lims_wvn.extent(1));
1108  wavelength_bounds = kdist.get_band_lims_wavelength();
1109 
1110  // Find the band index for the specified wavelength
1111  // Note that bands are stored in wavenumber space, units of cm-1, so if we are passed wavelength
1112  // in units of meters, we need a conversion factor of 10^2
1113  int nbnds = kdist.get_nband();
1114  int band_index = -1;
1115  Kokkos::parallel_reduce(nbnds, KOKKOS_LAMBDA(int ibnd, int& band_index_inner)
1116  {
1117  if (wavelength_bounds(0,ibnd) < wavelength_bounds(1,ibnd)) {
1118  if (wavelength_bounds(0,ibnd) <= wavelength * 1e2 && wavelength * 1e2 <= wavelength_bounds(1,ibnd)) {
1119  band_index_inner = ibnd;
1120  }
1121  } else {
1122  if (wavelength_bounds(0,ibnd) >= wavelength * 1e2 && wavelength * 1e2 >= wavelength_bounds(1,ibnd)) {
1123  band_index_inner = ibnd;
1124  }
1125  }
1126  }, Kokkos::Max<int>(band_index));
1127  return band_index;
1128 }

Referenced by get_wavelength_index_lw(), and get_wavelength_index_sw().

Here is the caller graph for this function:

◆ get_wavelength_index_lw()

int rrtmgp::get_wavelength_index_lw ( RealT  wavelength)
1098 { return get_wavelength_index(*k_dist_lw_k, wavelength); }
int get_wavelength_index(optical_props_t &kdist, RealT wavelength)
Definition: ERF_RRTMGP_Interface.cpp:1102
std::unique_ptr< gas_optics_t > k_dist_lw_k
Definition: ERF_RRTMGP_Interface.cpp:11
Here is the call graph for this function:

◆ get_wavelength_index_sw()

int rrtmgp::get_wavelength_index_sw ( RealT  wavelength)
1094 { return get_wavelength_index(*k_dist_sw_k, wavelength); }
Here is the call graph for this function:

◆ limit_to_bounds_1d()

template<typename InT , typename OutT , typename T >
void rrtmgp::limit_to_bounds_1d ( InT const &  arr_in,
T const  lower,
T const  upper,
OutT &  arr_out 
)
45 {
46  Kokkos::parallel_for(arr_out.size(), KOKKOS_LAMBDA(int i)
47  {
48  arr_out(i) = std::min(std::max(arr_in(i), lower), upper);
49  });
50 }

◆ limit_to_bounds_2d()

template<typename InT , typename OutT , typename T >
void rrtmgp::limit_to_bounds_2d ( InT const &  arr_in,
T const  lower,
T const  upper,
OutT &  arr_out 
)
59 {
60  const int ex0 = (int) arr_out.extent(0);
61  const int ex1 = (int) arr_out.extent(1);
62  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ex0, ex1}),
63  KOKKOS_LAMBDA (int i, int j)
64  {
65  arr_out(i, j) = std::min(std::max(arr_in(i, j), lower), upper);
66  });
67 }

Referenced by get_cloud_optics_lw(), get_cloud_optics_sw(), rrtmgp_lw(), and rrtmgp_sw().

Here is the caller graph for this function:

◆ mixing_ratio_to_cloud_mass()

template<class View1 , class View2 , class View3 , class View4 , class View5 >
void rrtmgp::mixing_ratio_to_cloud_mass ( View1 const &  mixing_ratio,
View2 const &  cloud_fraction,
View3 const &  rho,
View4 const &  dz,
View5 const &  cloud_mass 
)
17 {
18  const int ncol = mixing_ratio.extent(0);
19  const int nlay = mixing_ratio.extent(1);
20  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, nlay}),
21  KOKKOS_LAMBDA (int icol, int ilay)
22  {
23  // Compute in-cloud mixing ratio (mixing ratio of the cloudy part of the layer)
24  // NOTE: these thresholds (from E3SM) seem arbitrary, but included here for consistency
25  // This limits in-cloud mixing ratio to 0.005 kg/kg. According to note in cloud_diagnostics
26  // in EAM, this is consistent with limits in MG2. Is this true for P3?
27  if (cloud_fraction(icol,ilay) > 0) {
28  // Compute layer-integrated cloud mass (per unit area)
29  auto incloud_mixing_ratio = std::min(mixing_ratio(icol,ilay) /
30  std::max(0.0001, cloud_fraction(icol,ilay)), 0.005);
31  cloud_mass(icol,ilay) = incloud_mixing_ratio * rho(icol,ilay) * dz(icol,ilay);
32  } else {
33  cloud_mass(icol,ilay) = 0;
34  }
35  });
36 }

Referenced by Radiation::run_impl().

Here is the caller graph for this function:

◆ radiation_do()

bool rrtmgp::radiation_do ( const int  irad,
const int  nstep 
)
inline
103 {
104  // If irad == 0, then never do radiation;
105  // Otherwise, we always call radiation at the first step,
106  // and afterwards we do radiation if the timestep is divisible
107  // by irad
108  if (irad == 0) {
109  return false;
110  } else {
111  return ( (nstep == 0) || (nstep % irad == 0) );
112  }
113 }

◆ rrtmgp_finalize()

void rrtmgp::rrtmgp_finalize ( )
267 {
268  initialized = false;
269  k_dist_sw_k->finalize();
270  k_dist_lw_k->finalize();
271  cloud_optics_sw_k->finalize();
272  cloud_optics_lw_k->finalize();
273  k_dist_sw_k = nullptr;
274  k_dist_lw_k = nullptr;
275  cloud_optics_sw_k = nullptr;
276  cloud_optics_lw_k = nullptr;
277  pool_t::finalize();
278 }
std::unique_ptr< cloud_optics_t > cloud_optics_sw_k
Definition: ERF_RRTMGP_Interface.cpp:18
std::unique_ptr< cloud_optics_t > cloud_optics_lw_k
Definition: ERF_RRTMGP_Interface.cpp:19
bool initialized
Definition: ERF_RRTMGP_Interface.cpp:24

Referenced by Radiation::finalize_impl().

Here is the caller graph for this function:

◆ rrtmgp_initialize()

void rrtmgp::rrtmgp_initialize ( gas_concs_t gas_concs_k,
const std::string &  coefficients_file_sw,
const std::string &  coefficients_file_lw,
const std::string &  cloud_optics_file_sw,
const std::string &  cloud_optics_file_lw 
)
234 {
235  // Initialize Kokkos
236  if (!Kokkos::is_initialized()) { Kokkos::initialize(); }
237 
238  // Create objects for static ptrs
239  k_dist_sw_k = std::make_unique<gas_optics_t>();
240  k_dist_lw_k = std::make_unique<gas_optics_t>();
241  cloud_optics_sw_k = std::make_unique<cloud_optics_t>();
242  cloud_optics_lw_k = std::make_unique<cloud_optics_t>();
243 
244  // Load and initialize absorption coefficient data
245  load_and_init(*k_dist_sw_k, coefficients_file_sw, gas_concs_k);
246  load_and_init(*k_dist_lw_k, coefficients_file_lw, gas_concs_k);
247 
248  // Load and initialize cloud optical property look-up table information
249  load_cld_lutcoeff(*cloud_optics_sw_k, cloud_optics_file_sw);
250  load_cld_lutcoeff(*cloud_optics_lw_k, cloud_optics_file_lw);
251 
252  // Initialize kokkos rrtmgp pool allocator
253  const size_t nvar = 300;
254  const size_t nbnd = std::max(k_dist_sw_k->get_nband(),k_dist_sw_k->get_nband());
255  const size_t ncol = gas_concs_k.ncol;
256  const size_t nlay = gas_concs_k.nlay;
257  auto my_size_ref = static_cast<unsigned long>(nvar * ncol * nlay * nbnd);
258  pool_t::init(my_size_ref);
259 
260  // We are now initialized!
261  initialized = true;
262 }

Referenced by Radiation::initialize_impl().

Here is the caller graph for this function:

◆ rrtmgp_lw()

void rrtmgp::rrtmgp_lw ( const int  ncol,
const int  nlay,
gas_optics_t k_dist,
real2d_k p_lay,
real2d_k t_lay,
real2d_k p_lev,
real2d_k t_lev,
real1d_k t_sfc,
real2d_k emis_sfc,
real1d_k lw_src,
gas_concs_t gas_concs,
optical_props1_t aerosol,
optical_props1_t clouds,
fluxes_t fluxes,
fluxes_broadband_t clnclrsky_fluxes,
fluxes_broadband_t clrsky_fluxes,
fluxes_broadband_t clnsky_fluxes,
const bool  extra_clnclrsky_diag,
const bool  extra_clnsky_diag 
)
910 {
911  // Problem size
912  int nbnd = k_dist.get_nband();
913 
914  // Associate local pointers for fluxes
915  auto& flux_up = fluxes.flux_up;
916  auto& flux_dn = fluxes.flux_dn;
917  auto& bnd_flux_up = fluxes.bnd_flux_up;
918  auto& bnd_flux_dn = fluxes.bnd_flux_dn;
919  auto& clnclrsky_flux_up = clnclrsky_fluxes.flux_up;
920  auto& clnclrsky_flux_dn = clnclrsky_fluxes.flux_dn;
921  auto& clrsky_flux_up = clrsky_fluxes.flux_up;
922  auto& clrsky_flux_dn = clrsky_fluxes.flux_dn;
923  auto& clnsky_flux_up = clnsky_fluxes.flux_up;
924  auto& clnsky_flux_dn = clnsky_fluxes.flux_dn;
925 
926  // Reset fluxes to zero
927  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, nlay+1}),
928  KOKKOS_LAMBDA (int icol, int ilev)
929  {
930  flux_up(icol, ilev) = 0;
931  flux_dn(icol, ilev) = 0;
932  clnclrsky_flux_up(icol, ilev) = 0;
933  clnclrsky_flux_dn(icol, ilev) = 0;
934  clrsky_flux_up(icol, ilev) = 0;
935  clrsky_flux_dn(icol, ilev) = 0;
936  clnsky_flux_up(icol, ilev) = 0;
937  clnsky_flux_dn(icol, ilev) = 0;
938  });
939  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay+1, nbnd}),
940  KOKKOS_LAMBDA (int icol, int ilev, int ibnd)
941  {
942  bnd_flux_up(icol, ilev, ibnd) = 0;
943  bnd_flux_dn(icol, ilev, ibnd) = 0;
944  });
945 
946  // Allocate space for optical properties
947  optical_props1_t optics;
948  optics.alloc_1scl(ncol, nlay, k_dist);
949 
950  optical_props1_t optics_no_aerosols;
951  if (extra_clnsky_diag) {
952  // Allocate space for optical properties (no aerosols)
953  optics_no_aerosols.alloc_1scl(ncol, nlay, k_dist);
954  }
955 
956  bool top_at_1 = false;
957  Kokkos::parallel_reduce(1, KOKKOS_LAMBDA(int, bool& val)
958  {
959  val |= p_lay(0, 0) < p_lay(0, nlay-1);
960  }, Kokkos::LOr<bool>(top_at_1));
961 
962  // Boundary conditions
963  //=====================================================================
964  source_func_t lw_sources;
965  lw_sources.alloc(ncol, nlay, k_dist);
966 
967  /*
968  // Surface LW source
969  // AML NOTE: This is removed in EAMXX, LSM doesn't transfer its lw_src?
970  auto d_lw_src = lw_sources.sfc_source;
971  Kokkos::parallel_for(ncol, KOKKOS_LAMBDA (int icol)
972  {
973  d_lw_src(icol, 0) = lw_src(icol);
974  });
975  */
976 
977  // Surface temperature
978  // AML NOTE: We already populate this when initializing data
979 
980 
981  // Surface emissivity (transposed in RRTMGP)
982  // AML NOTE: This transfer was removed in EAMXX, LSM doesn't transfer its emis_sfc?
983  real2d_k emis_sfc_T("emis_sfc",nbnd,ncol);
984  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, nbnd}),
985  KOKKOS_LAMBDA (int icol, int ibnd)
986  {
987  emis_sfc_T(ibnd,icol) = emis_sfc(icol,ibnd);
988  });
989  //Kokkos::deep_copy(emis_sfc_T, 0.98);
990 
991  // Get Gaussian quadrature weights
992  // Weights and angle secants for first order (k=1) Gaussian quadrature.
993  // Values from Table 2, Clough et al, 1992, doi:10.1029/92JD01419
994  // after Abramowitz & Stegun 1972, page 921
995  int constexpr max_gauss_pts = 4;
996  RealT gauss_Ds_host_raw[max_gauss_pts][max_gauss_pts] = { {1.66, 1.18350343, 1.09719858, 1.06056257},
997  {0. , 2.81649655, 1.69338507, 1.38282560},
998  {0. , 0. , 4.70941630, 2.40148179},
999  {0. , 0. , 0. , 7.15513024} };
1000  realHost2d_k gauss_Ds_host(&gauss_Ds_host_raw[0][0], max_gauss_pts, max_gauss_pts);
1001 
1002  RealT gauss_wts_host_raw[max_gauss_pts][max_gauss_pts] = { {0.5, 0.3180413817, 0.2009319137, 0.1355069134},
1003  {0. , 0.1819586183, 0.2292411064, 0.2034645680},
1004  {0. , 0. , 0.0698269799, 0.1298475476},
1005  {0. , 0. , 0. , 0.0311809710} };
1006  realHost2d_k gauss_wts_host(&gauss_wts_host_raw[0][0],max_gauss_pts,max_gauss_pts);
1007 
1008  real2d_k gauss_Ds ("gauss_Ds" ,max_gauss_pts,max_gauss_pts);
1009  real2d_k gauss_wts("gauss_wts",max_gauss_pts,max_gauss_pts);
1010  Kokkos::deep_copy(gauss_Ds, gauss_Ds_host);
1011  Kokkos::deep_copy(gauss_wts, gauss_wts_host);
1012 
1013  // Limit temperatures for gas optics look-up tables
1014  real2d_k t_lay_limited("t_lay_limited", ncol, nlay );
1015  real2d_k t_lev_limited("t_lev_limited", ncol, nlay+1);
1016  limit_to_bounds_2d(t_lay, k_dist.get_temp_min(),
1017  k_dist.get_temp_max(), t_lay_limited);
1018  limit_to_bounds_2d(t_lev, k_dist.get_temp_min(),
1019  k_dist.get_temp_max(), t_lev_limited);
1020 
1021  // Do gas optics
1022  real3d_k col_gas("col_gas", ncol, nlay, k_dist.get_ngas()+1);
1023  k_dist.gas_optics(ncol, nlay, top_at_1, p_lay, p_lev, t_lay_limited,
1024  t_sfc, gas_concs, col_gas, optics, lw_sources, view_t<RealT**>(), t_lev_limited);
1025  if (extra_clnsky_diag) {
1026  k_dist.gas_optics(ncol, nlay, top_at_1, p_lay, p_lev, t_lay_limited,
1027  t_sfc, gas_concs, col_gas, optics_no_aerosols, lw_sources, view_t<RealT**>(), t_lev_limited);
1028  }
1029 
1030  if (extra_clnclrsky_diag) {
1031  // Compute clean-clear-sky fluxes before we add in aerosols and clouds
1032  rte_lw(max_gauss_pts, gauss_Ds, gauss_wts, optics, top_at_1, lw_sources, emis_sfc_T, clnclrsky_fluxes);
1033  }
1034 
1035  // Combine gas and aerosol optics
1036  aerosol.increment(optics);
1037 
1038  // Compute clear-sky fluxes before we add in clouds
1039  rte_lw(max_gauss_pts, gauss_Ds, gauss_wts, optics, top_at_1, lw_sources, emis_sfc_T, clrsky_fluxes);
1040 
1041  // Combine gas and cloud optics
1042  clouds.increment(optics);
1043 
1044  // Compute allsky fluxes
1045  rte_lw(max_gauss_pts, gauss_Ds, gauss_wts, optics, top_at_1, lw_sources, emis_sfc_T, fluxes);
1046 
1047  if (extra_clnsky_diag) {
1048  // First increment clouds in optics_no_aerosols
1049  clouds.increment(optics_no_aerosols);
1050  // Compute clean-sky fluxes
1051  rte_lw(max_gauss_pts, gauss_Ds, gauss_wts, optics_no_aerosols, top_at_1, lw_sources, emis_sfc_T, clnsky_fluxes);
1052  }
1053 }
Kokkos::View< RealT **, layout_t, KokkosHostDevice > realHost2d_k
Definition: ERF_Kokkos.H:16
Kokkos::View< T, layout_t, KokkosDefaultDevice > view_t
Definition: ERF_RRTMGP_Interface.H:26
SourceFuncLWK< RealT, layout_t, KokkosDefaultDevice > source_func_t
Definition: ERF_RRTMGP_Interface.H:38
@ t_sfc
Definition: ERF_NOAHMP.H:25

Referenced by rrtmgp_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ rrtmgp_main()

void rrtmgp::rrtmgp_main ( const int  ncol,
const int  nlay,
real2d_k p_lay,
real2d_k t_lay,
real2d_k p_lev,
real2d_k t_lev,
gas_concs_t gas_concs,
real2d_k sfc_alb_dir,
real2d_k sfc_alb_dif,
real1d_k mu0,
real1d_k t_sfc,
real2d_k emis_sfc,
real1d_k lw_src,
real2d_k lwp,
real2d_k iwp,
real2d_k rel,
real2d_k rei,
real2d_k cldfrac,
real3d_k aer_tau_sw,
real3d_k aer_ssa_sw,
real3d_k aer_asm_sw,
real3d_k aer_tau_lw,
real3d_k cld_tau_sw_bnd,
real3d_k cld_tau_lw_bnd,
real3d_k cld_tau_sw_gpt,
real3d_k cld_tau_lw_gpt,
real2d_k sw_flux_up,
real2d_k sw_flux_dn,
real2d_k sw_flux_dn_dir,
real2d_k lw_flux_up,
real2d_k lw_flux_dn,
real2d_k sw_clnclrsky_flux_up,
real2d_k sw_clnclrsky_flux_dn,
real2d_k sw_clnclrsky_flux_dn_dir,
real2d_k sw_clrsky_flux_up,
real2d_k sw_clrsky_flux_dn,
real2d_k sw_clrsky_flux_dn_dir,
real2d_k sw_clnsky_flux_up,
real2d_k sw_clnsky_flux_dn,
real2d_k sw_clnsky_flux_dn_dir,
real2d_k lw_clnclrsky_flux_up,
real2d_k lw_clnclrsky_flux_dn,
real2d_k lw_clrsky_flux_up,
real2d_k lw_clrsky_flux_dn,
real2d_k lw_clnsky_flux_up,
real2d_k lw_clnsky_flux_dn,
real3d_k sw_bnd_flux_up,
real3d_k sw_bnd_flux_dn,
real3d_k sw_bnd_flux_dn_dir,
real3d_k lw_bnd_flux_up,
real3d_k lw_bnd_flux_dn,
const RealT  tsi_scaling,
const bool  extra_clnclrsky_diag,
const bool  extra_clnsky_diag 
)
408 {
409  // Setup pointers to RRTMGP SW fluxes
410  fluxes_t fluxes_sw;
411  fluxes_sw.flux_up = sw_flux_up;
412  fluxes_sw.flux_dn = sw_flux_dn;
413  fluxes_sw.flux_dn_dir = sw_flux_dn_dir;
414  fluxes_sw.bnd_flux_up = sw_bnd_flux_up;
415  fluxes_sw.bnd_flux_dn = sw_bnd_flux_dn;
416  fluxes_sw.bnd_flux_dn_dir = sw_bnd_flux_dn_dir;
417  // Clean-clear-sky
418  fluxes_broadband_t clnclrsky_fluxes_sw;
419  clnclrsky_fluxes_sw.flux_up = sw_clnclrsky_flux_up;
420  clnclrsky_fluxes_sw.flux_dn = sw_clnclrsky_flux_dn;
421  clnclrsky_fluxes_sw.flux_dn_dir = sw_clnclrsky_flux_dn_dir;
422  // Clear-sky
423  fluxes_broadband_t clrsky_fluxes_sw;
424  clrsky_fluxes_sw.flux_up = sw_clrsky_flux_up;
425  clrsky_fluxes_sw.flux_dn = sw_clrsky_flux_dn;
426  clrsky_fluxes_sw.flux_dn_dir = sw_clrsky_flux_dn_dir;
427  // Clean-sky
428  fluxes_broadband_t clnsky_fluxes_sw;
429  clnsky_fluxes_sw.flux_up = sw_clnsky_flux_up;
430  clnsky_fluxes_sw.flux_dn = sw_clnsky_flux_dn;
431  clnsky_fluxes_sw.flux_dn_dir = sw_clnsky_flux_dn_dir;
432 
433  // Setup pointers to RRTMGP LW fluxes
434  fluxes_t fluxes_lw;
435  fluxes_lw.flux_up = lw_flux_up;
436  fluxes_lw.flux_dn = lw_flux_dn;
437  fluxes_lw.bnd_flux_up = lw_bnd_flux_up;
438  fluxes_lw.bnd_flux_dn = lw_bnd_flux_dn;
439  // Clean-clear-sky
440  fluxes_broadband_t clnclrsky_fluxes_lw;
441  clnclrsky_fluxes_lw.flux_up = lw_clnclrsky_flux_up;
442  clnclrsky_fluxes_lw.flux_dn = lw_clnclrsky_flux_dn;
443  // Clear-sky
444  fluxes_broadband_t clrsky_fluxes_lw;
445  clrsky_fluxes_lw.flux_up = lw_clrsky_flux_up;
446  clrsky_fluxes_lw.flux_dn = lw_clrsky_flux_dn;
447  // Clean-sky
448  fluxes_broadband_t clnsky_fluxes_lw;
449  clnsky_fluxes_lw.flux_up = lw_clnsky_flux_up;
450  clnsky_fluxes_lw.flux_dn = lw_clnsky_flux_dn;
451 
452  auto nswbands = k_dist_sw_k->get_nband();
453  auto nlwbands = k_dist_lw_k->get_nband();
454 
455  // Setup aerosol optical properties
456  optical_props2_t aerosol_sw;
457  optical_props1_t aerosol_lw;
458  aerosol_sw.init(k_dist_sw_k->get_band_lims_wavenumber());
459  aerosol_sw.alloc_2str(ncol, nlay);
460  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, nswbands}),
461  KOKKOS_LAMBDA (int icol, int ilay, int ibnd)
462  {
463  aerosol_sw.tau(icol,ilay,ibnd) = aer_tau_sw(icol,ilay,ibnd);
464  aerosol_sw.ssa(icol,ilay,ibnd) = aer_ssa_sw(icol,ilay,ibnd);
465  aerosol_sw.g (icol,ilay,ibnd) = aer_asm_sw(icol,ilay,ibnd);
466  });
467  aerosol_lw.init(k_dist_lw_k->get_band_lims_wavenumber());
468  aerosol_lw.alloc_1scl(ncol, nlay);
469  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, nlwbands}),
470  KOKKOS_LAMBDA (int icol, int ilay, int ibnd)
471  {
472  aerosol_lw.tau(icol,ilay,ibnd) = aer_tau_lw(icol,ilay,ibnd);
473  });
474 
475  // Convert cloud physical properties to optical properties for input to RRTMGP
476  optical_props2_t clouds_sw = get_cloud_optics_sw(ncol, nlay,
478  lwp, iwp, rel, rei);
479  optical_props1_t clouds_lw = get_cloud_optics_lw(ncol, nlay,
481  lwp, iwp, rel, rei);
482  Kokkos::deep_copy(cld_tau_sw_bnd, clouds_sw.tau);
483  Kokkos::deep_copy(cld_tau_lw_bnd, clouds_lw.tau);
484 
485  // Do subcolumn sampling to map bands -> gpoints based on cloud fraction and overlap assumption;
486  // This implements the Monte Carlo Independing Column Approximation by mapping only a single
487  // subcolumn (cloud state) to each gpoint.
488  auto nswgpts = k_dist_sw_k->get_ngpt();
489  auto clouds_sw_gpt = get_subsampled_clouds(ncol, nlay, nswbands, nswgpts,
490  clouds_sw, *k_dist_sw_k, cldfrac, p_lay);
491  // Longwave
492  auto nlwgpts = k_dist_lw_k->get_ngpt();
493  auto clouds_lw_gpt = get_subsampled_clouds(ncol, nlay, nlwbands, nlwgpts,
494  clouds_lw, *k_dist_lw_k, cldfrac, p_lay);
495 
496  // Copy cloud properties to outputs (is this needed, or can we just use pointers?)
497  // Alternatively, just compute and output a subcolumn cloud mask
498  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, nswgpts}),
499  KOKKOS_LAMBDA (int icol, int ilay, int igpt)
500  {
501  cld_tau_sw_gpt(icol,ilay,igpt) = clouds_sw_gpt.tau(icol,ilay,igpt);
502  });
503  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, nlwgpts}),
504  KOKKOS_LAMBDA (int icol, int ilay, int igpt)
505  {
506  cld_tau_lw_gpt(icol,ilay,igpt) = clouds_lw_gpt.tau(icol,ilay,igpt);
507  });
508 
509  // Do shortwave
510  rrtmgp_sw(ncol, nlay,
511  *k_dist_sw_k, p_lay, t_lay, p_lev, t_lev, gas_concs,
512  sfc_alb_dir, sfc_alb_dif, mu0, aerosol_sw, clouds_sw_gpt,
513  fluxes_sw, clnclrsky_fluxes_sw, clrsky_fluxes_sw, clnsky_fluxes_sw,
514  tsi_scaling, extra_clnclrsky_diag, extra_clnsky_diag);
515 
516  // Do longwave
517  rrtmgp_lw(ncol, nlay,
518  *k_dist_lw_k, p_lay, t_lay, p_lev, t_lev,
519  t_sfc, emis_sfc, lw_src,
520  gas_concs, aerosol_lw, clouds_lw_gpt,
521  fluxes_lw, clnclrsky_fluxes_lw, clrsky_fluxes_lw, clnsky_fluxes_lw,
522  extra_clnclrsky_diag, extra_clnsky_diag);
523 
524 }
FluxesBroadbandK< RealT, layout_t, KokkosDefaultDevice > fluxes_broadband_t
Definition: ERF_RRTMGP_Interface.H:34
FluxesBybandK< RealT, layout_t, KokkosDefaultDevice > fluxes_t
Definition: ERF_RRTMGP_Interface.H:33
@ sw_flux_dn
Definition: ERF_NOAHMP.H:31
@ lw_flux_dn
Definition: ERF_NOAHMP.H:32
optical_props1_t get_cloud_optics_lw(const int ncol, const int nlay, cloud_optics_t &cloud_optics, gas_optics_t &kdist, real2d_k &lwp, real2d_k &iwp, real2d_k &rel, real2d_k &rei)
Definition: ERF_RRTMGP_Interface.cpp:62
optical_props1_t get_subsampled_clouds(const int ncol, const int nlay, const int nbnd, const int ngpt, optical_props1_t &cloud_optics, gas_optics_t &kdist, real2d_k &cld, real2d_k &p_lay)
Definition: ERF_RRTMGP_Interface.cpp:164
optical_props2_t get_cloud_optics_sw(const int ncol, const int nlay, cloud_optics_t &cloud_optics, gas_optics_t &kdist, real2d_k &lwp, real2d_k &iwp, real2d_k &rel, real2d_k &rei)
Definition: ERF_RRTMGP_Interface.cpp:28
void rrtmgp_sw(const int ncol, const int nlay, gas_optics_t &k_dist, real2d_k &p_lay, real2d_k &t_lay, real2d_k &p_lev, real2d_k &t_lev, gas_concs_t &gas_concs, real2d_k &sfc_alb_dir, real2d_k &sfc_alb_dif, real1d_k &mu0, optical_props2_t &aerosol, optical_props2_t &clouds, fluxes_t &fluxes, fluxes_broadband_t &clnclrsky_fluxes, fluxes_broadband_t &clrsky_fluxes, fluxes_broadband_t &clnsky_fluxes, const RealT tsi_scaling, const bool extra_clnclrsky_diag, const bool extra_clnsky_diag)
Definition: ERF_RRTMGP_Interface.cpp:601
void rrtmgp_lw(const int ncol, const int nlay, gas_optics_t &k_dist, real2d_k &p_lay, real2d_k &t_lay, real2d_k &p_lev, real2d_k &t_lev, real1d_k &t_sfc, real2d_k &emis_sfc, real1d_k &lw_src, gas_concs_t &gas_concs, optical_props1_t &aerosol, optical_props1_t &clouds, fluxes_t &fluxes, fluxes_broadband_t &clnclrsky_fluxes, fluxes_broadband_t &clrsky_fluxes, fluxes_broadband_t &clnsky_fluxes, const bool extra_clnclrsky_diag, const bool extra_clnsky_diag)
Definition: ERF_RRTMGP_Interface.cpp:895

Referenced by Radiation::run_impl().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ rrtmgp_sw()

void rrtmgp::rrtmgp_sw ( const int  ncol,
const int  nlay,
gas_optics_t k_dist,
real2d_k p_lay,
real2d_k t_lay,
real2d_k p_lev,
real2d_k t_lev,
gas_concs_t gas_concs,
real2d_k sfc_alb_dir,
real2d_k sfc_alb_dif,
real1d_k mu0,
optical_props2_t aerosol,
optical_props2_t clouds,
fluxes_t fluxes,
fluxes_broadband_t clnclrsky_fluxes,
fluxes_broadband_t clrsky_fluxes,
fluxes_broadband_t clnsky_fluxes,
const RealT  tsi_scaling,
const bool  extra_clnclrsky_diag,
const bool  extra_clnsky_diag 
)
618 {
619  // Get problem sizes
620  int nbnd = k_dist.get_nband();
621  int ngpt = k_dist.get_ngpt();
622  int ngas = gas_concs.get_num_gases();
623 
624  // Associate local pointers for fluxes
625  auto& flux_up = fluxes.flux_up;
626  auto& flux_dn = fluxes.flux_dn;
627  auto& flux_dn_dir = fluxes.flux_dn_dir;
628  auto& bnd_flux_up = fluxes.bnd_flux_up;
629  auto& bnd_flux_dn = fluxes.bnd_flux_dn;
630  auto& bnd_flux_dn_dir = fluxes.bnd_flux_dn_dir;
631  auto& clnclrsky_flux_up = clnclrsky_fluxes.flux_up;
632  auto& clnclrsky_flux_dn = clnclrsky_fluxes.flux_dn;
633  auto& clnclrsky_flux_dn_dir = clnclrsky_fluxes.flux_dn_dir;
634  auto& clrsky_flux_up = clrsky_fluxes.flux_up;
635  auto& clrsky_flux_dn = clrsky_fluxes.flux_dn;
636  auto& clrsky_flux_dn_dir = clrsky_fluxes.flux_dn_dir;
637  auto& clnsky_flux_up = clnsky_fluxes.flux_up;
638  auto& clnsky_flux_dn = clnsky_fluxes.flux_dn;
639  auto& clnsky_flux_dn_dir = clnsky_fluxes.flux_dn_dir;
640 
641  // Reset fluxes to zero
642  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, nlay+1}),
643  KOKKOS_LAMBDA (int icol, int ilev)
644  {
645  flux_up (icol,ilev) = 0;
646  flux_dn (icol,ilev) = 0;
647  flux_dn_dir(icol,ilev) = 0;
648  clnclrsky_flux_up (icol,ilev) = 0;
649  clnclrsky_flux_dn (icol,ilev) = 0;
650  clnclrsky_flux_dn_dir(icol,ilev) = 0;
651  clrsky_flux_up (icol,ilev) = 0;
652  clrsky_flux_dn (icol,ilev) = 0;
653  clrsky_flux_dn_dir(icol,ilev) = 0;
654  clnsky_flux_up (icol,ilev) = 0;
655  clnsky_flux_dn (icol,ilev) = 0;
656  clnsky_flux_dn_dir(icol,ilev) = 0;
657  });
658  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay+1, nbnd}),
659  KOKKOS_LAMBDA (int icol, int ilev, int ibnd)
660  {
661  bnd_flux_up (icol,ilev,ibnd) = 0;
662  bnd_flux_dn (icol,ilev,ibnd) = 0;
663  bnd_flux_dn_dir(icol,ilev,ibnd) = 0;
664  });
665 
666  // Get daytime indices
667  int1d_k dayIndices("dayIndices", ncol);
668  Kokkos::deep_copy(dayIndices, -1);
669 
670  // Serialized for now.
671  int nday = 0;
672  Kokkos::parallel_reduce(1, KOKKOS_LAMBDA(int, int& nday_inner)
673  {
674  for (int icol = 0; icol < ncol; ++icol) {
675  if (mu0(icol) > 0) {
676  dayIndices(nday_inner++) = icol;
677  }
678  }
679  }, Kokkos::Sum<int>(nday));
680 
681  // Copy data back to the device
682  if (nday == 0) {
683  // No daytime columns in this chunk, skip the rest of this routine
684  return;
685  }
686 
687  // Subset mu0
688  real1d_k mu0_day("mu0_day", nday);
689  Kokkos::parallel_for(nday, KOKKOS_LAMBDA(int iday)
690  {
691  mu0_day(iday) = mu0(dayIndices(iday));
692  });
693 
694  // subset state variables
695  real2d_k p_lay_day("p_lay_day", nday, nlay);
696  real2d_k t_lay_day("t_lay_day", nday, nlay);
697  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {nday, nlay}),
698  KOKKOS_LAMBDA (int iday, int ilay)
699  {
700  p_lay_day(iday,ilay) = p_lay(dayIndices(iday),ilay);
701  t_lay_day(iday,ilay) = t_lay(dayIndices(iday),ilay);
702  });
703  real2d_k p_lev_day("p_lev_day", nday, nlay+1);
704  real2d_k t_lev_day("t_lev_day", nday, nlay+1);
705  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {nday, nlay+1}),
706  KOKKOS_LAMBDA (int iday, int ilay)
707  {
708  p_lev_day(iday,ilay) = p_lev(dayIndices(iday),ilay);
709  t_lev_day(iday,ilay) = t_lev(dayIndices(iday),ilay);
710  });
711 
712  // Subset gases
713  auto gas_names = gas_concs.get_gas_names();
714  gas_concs_t gas_concs_day;
715  gas_concs_day.init(gas_names, nday, nlay);
716  for (int igas = 0; igas < ngas; igas++) {
717  real2d_k vmr_day("vmr_day", nday, nlay);
718  real2d_k vmr("vmr" , ncol, nlay);
719  gas_concs.get_vmr(gas_names[igas], vmr);
720  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {nday, nlay}),
721  KOKKOS_LAMBDA (int iday, int ilay)
722  {
723  vmr_day(iday,ilay) = vmr(dayIndices(iday),ilay);
724  });
725  gas_concs_day.set_vmr(gas_names[igas], vmr_day);
726  }
727 
728  // Subset aerosol optics
729  optical_props2_t aerosol_day;
730  aerosol_day.init(k_dist.get_band_lims_wavenumber());
731  aerosol_day.alloc_2str(nday, nlay);
732  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {nday, nlay, nbnd}),
733  KOKKOS_LAMBDA (int iday, int ilay, int ibnd)
734  {
735  aerosol_day.tau(iday,ilay,ibnd) = aerosol.tau(dayIndices(iday),ilay,ibnd);
736  aerosol_day.ssa(iday,ilay,ibnd) = aerosol.ssa(dayIndices(iday),ilay,ibnd);
737  aerosol_day.g (iday,ilay,ibnd) = aerosol.g (dayIndices(iday),ilay,ibnd);
738  });
739 
740  // Subset cloud optics
741  // TODO: nbnd -> ngpt once we pass sub-sampled cloud state
742  optical_props2_t clouds_day;
743  clouds_day.init(k_dist.get_band_lims_wavenumber(), k_dist.get_band_lims_gpoint());
744  clouds_day.alloc_2str(nday, nlay);
745  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {nday, nlay, ngpt}),
746  KOKKOS_LAMBDA (int iday, int ilay, int igpt)
747  {
748  clouds_day.tau(iday,ilay,igpt) = clouds.tau(dayIndices(iday),ilay,igpt);
749  clouds_day.ssa(iday,ilay,igpt) = clouds.ssa(dayIndices(iday),ilay,igpt);
750  clouds_day.g (iday,ilay,igpt) = clouds.g (dayIndices(iday),ilay,igpt);
751  });
752 
753  // RRTMGP assumes surface albedos have a screwy dimension ordering
754  // for some strange reason, so we need to transpose these; also do
755  // daytime subsetting in the same kernel
756  real2d_k sfc_alb_dir_T("sfc_alb_dir", nbnd, nday);
757  real2d_k sfc_alb_dif_T("sfc_alb_dif", nbnd, nday);
758  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {nbnd, nday}),
759  KOKKOS_LAMBDA (int ibnd, int icol)
760  {
761  sfc_alb_dir_T(ibnd,icol) = sfc_alb_dir(dayIndices(icol),ibnd);
762  sfc_alb_dif_T(ibnd,icol) = sfc_alb_dif(dayIndices(icol),ibnd);
763  });
764 
765  // Temporaries we need for daytime-only fluxes
766  real2d_k flux_up_day("flux_up_day", nday, nlay+1);
767  real2d_k flux_dn_day("flux_dn_day", nday, nlay+1);
768  real2d_k flux_dn_dir_day("flux_dn_dir_day", nday, nlay+1);
769  real3d_k bnd_flux_up_day("bnd_flux_up_day", nday, nlay+1, nbnd);
770  real3d_k bnd_flux_dn_day("bnd_flux_dn_day", nday, nlay+1, nbnd);
771  real3d_k bnd_flux_dn_dir_day("bnd_flux_dn_dir_day", nday, nlay+1, nbnd);
772  fluxes_t fluxes_day;
773  fluxes_day.flux_up = flux_up_day;
774  fluxes_day.flux_dn = flux_dn_day;
775  fluxes_day.flux_dn_dir = flux_dn_dir_day;
776  fluxes_day.bnd_flux_up = bnd_flux_up_day;
777  fluxes_day.bnd_flux_dn = bnd_flux_dn_day;
778  fluxes_day.bnd_flux_dn_dir = bnd_flux_dn_dir_day;
779 
780  // Allocate space for optical properties
781  optical_props2_t optics;
782  optics.alloc_2str(nday, nlay, k_dist);
783 
784  optical_props2_t optics_no_aerosols;
785  if (extra_clnsky_diag) {
786  // Allocate space for optical properties (no aerosols)
787  optics_no_aerosols.alloc_2str(nday, nlay, k_dist);
788  }
789 
790  // Limit temperatures for gas optics look-up tables
791  real2d_k t_lay_limited("t_lay_limited", nday, nlay);
792  limit_to_bounds_2d(t_lay_day, k_dist_sw_k->get_temp_min(),
793  k_dist_sw_k->get_temp_max(), t_lay_limited);
794 
795  // Do gas optics
796  real2d_k toa_flux("toa_flux", nday, ngpt);
797  real3d_k col_gas("col_gas", ncol, nlay, k_dist.get_ngas()+1);
798  bool top_at_1 = false;
799  Kokkos::parallel_reduce(1, KOKKOS_LAMBDA(int, bool& val)
800  {
801  val |= p_lay(0, 0) < p_lay(0, nlay-1);
802  }, Kokkos::LOr<bool>(top_at_1));
803 
804  k_dist.gas_optics(nday, nlay, top_at_1, p_lay_day, p_lev_day,
805  t_lay_limited, gas_concs_day, col_gas, optics, toa_flux);
806  if (extra_clnsky_diag) {
807  k_dist.gas_optics(nday, nlay, top_at_1, p_lay_day, p_lev_day,
808  t_lay_limited, gas_concs_day, col_gas, optics_no_aerosols, toa_flux);
809  }
810 
811  // Apply tsi_scaling
812  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {nday, ngpt}),
813  KOKKOS_LAMBDA (int iday, int igpt)
814  {
815  toa_flux(iday,igpt) = tsi_scaling * toa_flux(iday,igpt);
816  });
817 
818  if (extra_clnclrsky_diag) {
819  // Compute clear-clean-sky (just gas) fluxes on daytime columns
820  rte_sw(optics, top_at_1, mu0_day, toa_flux, sfc_alb_dir_T, sfc_alb_dif_T, fluxes_day);
821  // Expand daytime fluxes to all columns
822  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {nday, nlay+1}),
823  KOKKOS_LAMBDA (int iday, int ilev)
824  {
825  int icol = dayIndices(iday);
826  clnclrsky_flux_up (icol,ilev) = flux_up_day (iday,ilev);
827  clnclrsky_flux_dn (icol,ilev) = flux_dn_day (iday,ilev);
828  clnclrsky_flux_dn_dir(icol,ilev) = flux_dn_dir_day(iday,ilev);
829  });
830  }
831 
832  // Combine gas and aerosol optics
833  aerosol_day.delta_scale();
834  aerosol_day.increment(optics);
835 
836  // Compute clearsky (gas + aerosol) fluxes on daytime columns
837  rte_sw(optics, top_at_1, mu0_day, toa_flux, sfc_alb_dir_T, sfc_alb_dif_T, fluxes_day);
838 
839  // Expand daytime fluxes to all columns
840  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {nday, nlay+1}),
841  KOKKOS_LAMBDA (int iday, int ilev)
842  {
843  int icol = dayIndices(iday);
844  clrsky_flux_up (icol,ilev) = flux_up_day (iday,ilev);
845  clrsky_flux_dn (icol,ilev) = flux_dn_day (iday,ilev);
846  clrsky_flux_dn_dir(icol,ilev) = flux_dn_dir_day(iday,ilev);
847  });
848 
849  // Now merge in cloud optics and do allsky calculations
850 
851  // Combine gas and cloud optics
852  clouds_day.delta_scale();
853  clouds_day.increment(optics);
854 
855  // Compute fluxes on daytime columns
856  rte_sw(optics, top_at_1, mu0_day, toa_flux, sfc_alb_dir_T, sfc_alb_dif_T, fluxes_day);
857 
858  // Expand daytime fluxes to all columns
859  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {nday, nlay+1}),
860  KOKKOS_LAMBDA (int iday, int ilev)
861  {
862  int icol = dayIndices(iday);
863  flux_up (icol,ilev) = flux_up_day (iday,ilev);
864  flux_dn (icol,ilev) = flux_dn_day (iday,ilev);
865  flux_dn_dir(icol,ilev) = flux_dn_dir_day(iday,ilev);
866  });
867  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {nday, nlay+1, nbnd}),
868  KOKKOS_LAMBDA (int iday, int ilev, int ibnd)
869  {
870  int icol = dayIndices(iday);
871  bnd_flux_up (icol,ilev,ibnd) = bnd_flux_up_day (iday,ilev,ibnd);
872  bnd_flux_dn (icol,ilev,ibnd) = bnd_flux_dn_day (iday,ilev,ibnd);
873  bnd_flux_dn_dir(icol,ilev,ibnd) = bnd_flux_dn_dir_day(iday,ilev,ibnd);
874  });
875 
876  if (extra_clnsky_diag) {
877  // First increment clouds in optics_no_aerosols
878  clouds_day.increment(optics_no_aerosols);
879  // Compute cleansky (gas + clouds) fluxes on daytime columns
880  rte_sw(optics_no_aerosols, top_at_1, mu0_day, toa_flux, sfc_alb_dir_T, sfc_alb_dif_T, fluxes_day);
881  // Expand daytime fluxes to all columns
882  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {nday, nlay+1}),
883  KOKKOS_LAMBDA (int iday, int ilev)
884  {
885  int icol = dayIndices(iday);
886  clnsky_flux_up (icol,ilev) = flux_up_day (iday,ilev);
887  clnsky_flux_dn (icol,ilev) = flux_dn_day (iday,ilev);
888  clnsky_flux_dn_dir(icol,ilev) = flux_dn_dir_day(iday,ilev);
889  });
890  }
891 }
GasConcsK< RealT, layout_t, KokkosDefaultDevice > gas_concs_t
Definition: ERF_RRTMGP_Interface.H:32

Referenced by rrtmgp_main().

Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ cloud_optics_lw_k

std::unique_ptr< cloud_optics_t > rrtmgp::cloud_optics_lw_k

◆ cloud_optics_sw_k

std::unique_ptr< cloud_optics_t > rrtmgp::cloud_optics_sw_k

◆ initialized

bool rrtmgp::initialized = false

◆ k_dist_lw_k

std::unique_ptr< gas_optics_t > rrtmgp::k_dist_lw_k

◆ k_dist_sw_k

◆ kokkos_mem_pool

pool_t rrtmgp::kokkos_mem_pool