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

◆ 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 amrex::Real(0.7) micron, or 14286 cm^-one
299  const RealT visible_wavenumber_threshold = amrex::Real(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) = myhalf*(sfc_alb_dir_vis(icol) + sfc_alb_dir_nir(icol));
319  sfc_alb_dif(icol,ibnd) = myhalf*(sfc_alb_dif_vis(icol) + sfc_alb_dif_nir(icol));
320  }
321  });
322 }
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
@ sfc_alb_dir_vis
Definition: ERF_NOAHMP.H:27
@ sfc_alb_dif_nir
Definition: ERF_NOAHMP.H:30
@ sfc_alb_dir_nir
Definition: ERF_NOAHMP.H:28
@ sfc_alb_dif_vis
Definition: ERF_NOAHMP.H:29
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)) + myhalf * sw_bnd_flux_dir(i+1,kbot,10);
341  //sfc_flux_dir_vis(i) = sum(sw_bnd_flux_dir(i+1,kbot,11:14)) + myhalf * sw_bnd_flux_dir(i+1,kbot,10);
342  //sfc_flux_dif_nir(i) = sum(sw_bnd_flux_dif(i+1,kbot,1:9)) + myhalf * sw_bnd_flux_dif(i+1,kbot,10);
343  //sfc_flux_dif_vis(i) = sum(sw_bnd_flux_dif(i+1,kbot,11:14)) + myhalf * 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 amrex::Real(0.7) micron, or 14286 cm^-one
352  const RealT visible_wavenumber_threshold = amrex::Real(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 myhalf
372  // the flux in visible and myhalf in near-infrared fluxes
373  sfc_flux_dir_vis(icol) += myhalf * sw_bnd_flux_dir(icol,kbot,ibnd);
374  sfc_flux_dif_vis(icol) += myhalf * sw_bnd_flux_dif(icol,kbot,ibnd);
375  sfc_flux_dir_nir(icol) += myhalf * sw_bnd_flux_dir(icol,kbot,ibnd);
376  sfc_flux_dif_nir(icol) += myhalf * 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 
)
1057 {
1058  // Subcolumn binary cld mask; if any layers with pressure between pmin and pmax are cloudy
1059  // then 2d subcol mask is 1, otherwise it is 0
1060  real2d_k subcol_mask("subcol_mask", ncol, ngpt);
1061  Kokkos::deep_copy(subcol_mask, zero);
1062  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, ngpt}),
1063  KOKKOS_LAMBDA (int icol, int ilay, int igpt)
1064  {
1065  // NOTE: using plev would need to assume level ordering (top to bottom or bottom to top), but
1066  // using play/pmid does not
1067  if (cld_tau_gpt(icol,ilay,igpt) > 0 && pmid(icol,ilay) >= pmin && pmid(icol,ilay) < pmax) {
1068  subcol_mask(icol,igpt) = one;
1069  }
1070  });
1071  // Compute average over subcols to get cloud area
1072  auto ngpt_inv = one / ngpt;
1073  Kokkos::deep_copy(cld_area, 0);
1074  Kokkos::parallel_for(ncol, KOKKOS_LAMBDA(int icol)
1075  {
1076  // This loop needs to be serial because of the atomic reduction
1077  for (int igpt = 0; igpt < ngpt; ++igpt) {
1078  cld_area(icol) += subcol_mask(icol,igpt) * ngpt_inv;
1079  }
1080  });
1081 }
Kokkos::View< RealT **, layout_t, KokkosDefaultDevice > real2d_k
Definition: ERF_Kokkos.H:19

◆ 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  // NOTE: This calculation defines Fnet = Up - Dn
93  // H = dT/dt = g/Cp dF/dP
94  // Using dF/dz * dz/dP = dF/dz * (1/(-rho*g)) = dF/dP
95  // We have: dT/dt = dF/dz * (1/(-rho*Cp))
96  heating_rate(icol,ilay) = ( flux_up(icol,ilay+1) - flux_up(icol,ilay)
97  - flux_dn(icol,ilay+1) + flux_dn(icol,ilay) )
98  / ( -rho(icol,ilay) * dz(icol,ilay) * Cp_d );
99  });
100 }
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:23
rho
Definition: ERF_InitCustomPert_Bubble.H:106
const Real dz
Definition: ERF_InitCustomPert_Bubble.H:25

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 
)
502 {
503  // Routine will return subcolumn mask with values of 0 indicating no cloud, 1 indicating cloud
504  int3d_k subcolumn_mask("subcolumn_mask", ncol, nlay, ngpt);
505 
506  // Subcolumn generators are a means for producing a variable x(i,j,k), where
507  //
508  // c(i,j,k) = 1 for x(i,j,k) > 1 - cldf(i,j)
509  // c(i,j,k) = 0 for x(i,j,k) <= 1 - cldf(i,j)
510  //
511  // I am going to call this "cldx" to be just slightly less ambiguous
512  real3d_k cldx("cldx", ncol, nlay, ngpt);
513 
514  // Apply overlap assumption to set cldx
515  if (overlap_option == 0) { // Dummy mask, always cloudy
516  Kokkos::deep_copy(cldx, 1);
517  } else { // Default case, maximum-random overlap
518  // Maximum-random overlap:
519  // Uses essentially the algorithm described in eq (14) in Raisanen et al. 2004,
520  // https://rmets.onlinelibrary.wiley.com/doi/epdf/amrex::Real(10.1256)/qj.03.99. Also the same
521  // algorithm used in RRTMG implementation of maximum-random overlap (see
522  // https://github.com/AER-RC/RRTMG_SW/blob/master/src/mcica_subcol_gen_sw.f90)
523  //
524  // First, fill cldx with random numbers.
525  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, ngpt}),
526  KOKKOS_LAMBDA (int icol, int ilay, int igpt)
527  {
528  conv::Random rand(seeds(icol) + ilay*ngpt + igpt);
529  cldx(icol,ilay,igpt) = rand.genFP<RealT>();
530  });
531 
532  // Step down columns and apply algorithm from eq (14)
533  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, ngpt}),
534  KOKKOS_LAMBDA (int icol, int igpt)
535  {
536  for (int ilay = 1; ilay < nlay; ilay++) {
537  // Check cldx in level above and see if it satisfies conditions to create a cloudy subcolumn
538  if (cldx(icol,ilay-1,igpt) > one - cldf(icol,ilay-1)) {
539  // Cloudy subcolumn above, use same random number here so that clouds in these two adjacent
540  // layers are maximimally overlapped
541  cldx(icol,ilay,igpt) = cldx(icol,ilay-1,igpt);
542  } else {
543  // Cloud-less above, use new random number so that clouds are distributed
544  // randomly in this layer. Need to scale new random number to range
545  // [0, one - cldf(ilay-1)] because we have artificially changed the distribution
546  // of random numbers in this layer with the above branch of the conditional,
547  // which would otherwise inflate cloud fraction in this layer.
548  cldx(icol,ilay,igpt) = cldx(icol,ilay ,igpt) * (one - cldf(icol,ilay-1));
549  }
550  }
551  });
552  }
553 
554  // Use cldx array to create subcolumn mask
555  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay, ngpt}),
556  KOKKOS_LAMBDA (int icol, int ilay, int igpt)
557  {
558  if (cldx(icol,ilay,igpt) > one - cldf(icol,ilay)) {
559  subcolumn_mask(icol,ilay,igpt) = 1;
560  } else {
561  subcolumn_mask(icol,ilay,igpt) = 0;
562  }
563  });
564  return subcolumn_mask;
565 }
Kokkos::View< RealT ***, layout_t, KokkosDefaultDevice > real3d_k
Definition: ERF_Kokkos.H:20
Kokkos::View< int ***, layout_t, KokkosDefaultDevice > int3d_k
Definition: ERF_Kokkos.H:23

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, zero);
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:21
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:496
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, zero);
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) = amrex::Real(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 
)
1095 {
1096  // Get wavelength bounds for all wavelength bands
1097  auto band_lims_wvn = kdist.get_band_lims_wavenumber();
1098  real2d_k wavelength_bounds("wavelength_bounds",band_lims_wvn.extent(0), band_lims_wvn.extent(1));
1099  wavelength_bounds = kdist.get_band_lims_wavelength();
1100 
1101  // Find the band index for the specified wavelength
1102  // Note that bands are stored in wavenumber space, units of cm-1, so if we are passed wavelength
1103  // in units of meters, we need a conversion factor of 10^2
1104  int nbnds = kdist.get_nband();
1105  int band_index = -1;
1106  Kokkos::parallel_reduce(nbnds, KOKKOS_LAMBDA(int ibnd, int& band_index_inner)
1107  {
1108  if (wavelength_bounds(0,ibnd) < wavelength_bounds(1,ibnd)) {
1109  if (wavelength_bounds(0,ibnd) <= wavelength * 1e2 && wavelength * 1e2 <= wavelength_bounds(1,ibnd)) {
1110  band_index_inner = ibnd;
1111  }
1112  } else {
1113  if (wavelength_bounds(0,ibnd) >= wavelength * 1e2 && wavelength * 1e2 >= wavelength_bounds(1,ibnd)) {
1114  band_index_inner = ibnd;
1115  }
1116  }
1117  }, Kokkos::Max<int>(band_index));
1118  return band_index;
1119 }
Real wavelength
Definition: ERF_InitCustomPert_MovingTerrain.H:5

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)
int get_wavelength_index(optical_props_t &kdist, RealT wavelength)
Definition: ERF_RRTMGP_Interface.cpp:1093
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)
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 amrex::Real(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(amrex::Real(0.0001), cloud_fraction(icol,ilay)), amrex::Real(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
107 {
108  // If irad == 0, then never do radiation;
109  // Otherwise, we always call radiation at the first step,
110  // and afterwards we do radiation if the timestep is divisible
111  // by irad
112  if (irad == 0) {
113  return false;
114  } else {
115  return ( (nstep == 0) || (nstep % irad == 0) );
116  }
117 }

◆ 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::~Radiation().

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 = 12;
254  const size_t ngpt = std::max(k_dist_sw_k->get_ngpt(),k_dist_lw_k->get_ngpt());
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 * ngpt);
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,
real1d_k sfc_emis,
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 
)
890 {
891  // Problem size
892  int nbnd = k_dist.get_nband();
893 
894  // Associate local pointers for fluxes
895  auto& flux_up = fluxes.flux_up;
896  auto& flux_dn = fluxes.flux_dn;
897  auto& bnd_flux_up = fluxes.bnd_flux_up;
898  auto& bnd_flux_dn = fluxes.bnd_flux_dn;
899  auto& clnclrsky_flux_up = clnclrsky_fluxes.flux_up;
900  auto& clnclrsky_flux_dn = clnclrsky_fluxes.flux_dn;
901  auto& clrsky_flux_up = clrsky_fluxes.flux_up;
902  auto& clrsky_flux_dn = clrsky_fluxes.flux_dn;
903  auto& clnsky_flux_up = clnsky_fluxes.flux_up;
904  auto& clnsky_flux_dn = clnsky_fluxes.flux_dn;
905 
906  // Reset fluxes to zero
907  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, nlay+1}),
908  KOKKOS_LAMBDA (int icol, int ilev)
909  {
910  flux_up(icol, ilev) = 0;
911  flux_dn(icol, ilev) = 0;
912  clrsky_flux_up(icol, ilev) = 0;
913  clrsky_flux_dn(icol, ilev) = 0;
914  });
915  if (extra_clnclrsky_diag) {
916  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, nlay+1}),
917  KOKKOS_LAMBDA (int icol, int ilev)
918  {
919  clnclrsky_flux_up(icol, ilev) = 0;
920  clnclrsky_flux_dn(icol, ilev) = 0;
921  });
922  }
923  if (extra_clnsky_diag) {
924  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, nlay+1}),
925  KOKKOS_LAMBDA (int icol, int ilev)
926  {
927  clnsky_flux_up(icol, ilev) = 0;
928  clnsky_flux_dn(icol, ilev) = 0;
929  });
930  }
931  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0}, {ncol, nlay+1, nbnd}),
932  KOKKOS_LAMBDA (int icol, int ilev, int ibnd)
933  {
934  bnd_flux_up(icol, ilev, ibnd) = 0;
935  bnd_flux_dn(icol, ilev, ibnd) = 0;
936  });
937 
938  // Allocate space for optical properties
939  optical_props1_t optics;
940  optics.alloc_1scl(ncol, nlay, k_dist);
941 
942  optical_props1_t optics_no_aerosols;
943  if (extra_clnsky_diag) {
944  // Allocate space for optical properties (no aerosols)
945  optics_no_aerosols.alloc_1scl(ncol, nlay, k_dist);
946  }
947 
948  bool top_at_1 = false;
949  Kokkos::parallel_reduce(1, KOKKOS_LAMBDA(int, bool& val)
950  {
951  val |= p_lay(0, 0) < p_lay(0, nlay-1);
952  }, Kokkos::LOr<bool>(top_at_1));
953 
954  // Boundary conditions
955  //=====================================================================
956  source_func_t lw_sources;
957  lw_sources.alloc(ncol, nlay, k_dist);
958 
959  /*
960  // Surface LW source
961  // AML NOTE: This is removed in EAMXX, LSM doesn't transfer its lw_src?
962  auto d_lw_src = lw_sources.sfc_source;
963  Kokkos::parallel_for(ncol, KOKKOS_LAMBDA (int icol)
964  {
965  d_lw_src(icol, 0) = lw_src(icol);
966  });
967  */
968 
969  // Surface temperature
970  // AML NOTE: We already populate this when initializing data
971 
972 
973  // Surface emissivity (transposed in RRTMGP)
974  // AML NOTE: This transfer was removed in EAMXX, LSM doesn't transfer its emis_sfc?
975  real2d_k emis_sfc_T("emis_sfc",nbnd,ncol);
976  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {ncol, nbnd}),
977  KOKKOS_LAMBDA (int icol, int ibnd)
978  {
979  emis_sfc_T(ibnd,icol) = sfc_emis(icol);
980  });
981  //Kokkos::deep_copy(emis_sfc_T, amrex::Real(0.98));
982 
983  // Get Gaussian quadrature weights
984  // Weights and angle secants for first order (k=1) Gaussian quadrature.
985  // Values from Table 2, Clough et al, 1992, doi:Real(10.1029)/92JD01419
986  // after Abramowitz & Stegun 1972, page 921
987  int constexpr max_gauss_pts = 4;
988  RealT gauss_Ds_host_raw[max_gauss_pts][max_gauss_pts] = { {amrex::Real(1.66), amrex::Real(1.18350343), amrex::Real(1.09719858), amrex::Real(1.06056257)},
989  {zero , amrex::Real(2.81649655), amrex::Real(1.69338507), amrex::Real(1.38282560)},
990  {zero , zero , amrex::Real(4.70941630), amrex::Real(2.40148179)},
991  {zero , zero , zero , amrex::Real(7.15513024)} };
992  realHost2d_k gauss_Ds_host(&gauss_Ds_host_raw[0][0], max_gauss_pts, max_gauss_pts);
993 
994  RealT gauss_wts_host_raw[max_gauss_pts][max_gauss_pts] = { {myhalf, amrex::Real(0.3180413817), amrex::Real(0.2009319137), amrex::Real(0.1355069134)},
995  {zero , amrex::Real(0.1819586183), amrex::Real(0.2292411064), amrex::Real(0.2034645680)},
996  {zero , zero , amrex::Real(0.0698269799), amrex::Real(0.1298475476)},
997  {zero , zero , zero , amrex::Real(0.0311809710)} };
998  realHost2d_k gauss_wts_host(&gauss_wts_host_raw[0][0],max_gauss_pts,max_gauss_pts);
999 
1000  real2d_k gauss_Ds ("gauss_Ds" ,max_gauss_pts,max_gauss_pts);
1001  real2d_k gauss_wts("gauss_wts",max_gauss_pts,max_gauss_pts);
1002  Kokkos::deep_copy(gauss_Ds, gauss_Ds_host);
1003  Kokkos::deep_copy(gauss_wts, gauss_wts_host);
1004 
1005  // Limit temperatures for gas optics look-up tables
1006  real2d_k t_lay_limited("t_lay_limited", ncol, nlay );
1007  real2d_k t_lev_limited("t_lev_limited", ncol, nlay+1);
1008  limit_to_bounds_2d(t_lay, k_dist.get_temp_min(),
1009  k_dist.get_temp_max(), t_lay_limited);
1010  limit_to_bounds_2d(t_lev, k_dist.get_temp_min(),
1011  k_dist.get_temp_max(), t_lev_limited);
1012 
1013  // Do gas optics
1014  real3d_k col_gas("col_gas", ncol, nlay, k_dist.get_ngas()+1);
1015  k_dist.gas_optics(ncol, nlay, top_at_1, p_lay, p_lev, t_lay_limited,
1016  t_sfc, gas_concs, col_gas, optics, lw_sources, view_t<RealT**>(), t_lev_limited);
1017  if (extra_clnsky_diag) {
1018  k_dist.gas_optics(ncol, nlay, top_at_1, p_lay, p_lev, t_lay_limited,
1019  t_sfc, gas_concs, col_gas, optics_no_aerosols, lw_sources, view_t<RealT**>(), t_lev_limited);
1020  }
1021 
1022  if (extra_clnclrsky_diag) {
1023  // Compute clean-clear-sky fluxes before we add in aerosols and clouds
1024  rte_lw(max_gauss_pts, gauss_Ds, gauss_wts, optics, top_at_1, lw_sources, emis_sfc_T, clnclrsky_fluxes);
1025  }
1026 
1027  // Combine gas and aerosol optics
1028  aerosol.increment(optics);
1029 
1030  // Compute clear-sky fluxes before we add in clouds
1031  rte_lw(max_gauss_pts, gauss_Ds, gauss_wts, optics, top_at_1, lw_sources, emis_sfc_T, clrsky_fluxes);
1032 
1033  // Combine gas and cloud optics
1034  clouds.increment(optics);
1035 
1036  // Compute allsky fluxes
1037  rte_lw(max_gauss_pts, gauss_Ds, gauss_wts, optics, top_at_1, lw_sources, emis_sfc_T, fluxes);
1038 
1039  if (extra_clnsky_diag) {
1040  // First increment clouds in optics_no_aerosols
1041  clouds.increment(optics_no_aerosols);
1042  // Compute clean-sky fluxes
1043  rte_lw(max_gauss_pts, gauss_Ds, gauss_wts, optics_no_aerosols, top_at_1, lw_sources, emis_sfc_T, clnsky_fluxes);
1044  }
1045 }
Kokkos::View< RealT **, layout_t, KokkosHostDevice > realHost2d_k
Definition: ERF_Kokkos.H:17
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
@ sfc_emis
Definition: ERF_NOAHMP.H:26
@ 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,
real1d_k sfc_emis,
real1d_k lw_src,
real2d_k lwp,
real2d_k iwp,
real2d_k rel,
real2d_k rei,
real2d_k cldfrac,
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 
)
404 {
405  // Setup pointers to RRTMGP SW fluxes
406  fluxes_t fluxes_sw;
407  fluxes_sw.flux_up = sw_flux_up;
408  fluxes_sw.flux_dn = sw_flux_dn;
409  fluxes_sw.flux_dn_dir = sw_flux_dn_dir;
410  fluxes_sw.bnd_flux_up = sw_bnd_flux_up;
411  fluxes_sw.bnd_flux_dn = sw_bnd_flux_dn;
412  fluxes_sw.bnd_flux_dn_dir = sw_bnd_flux_dn_dir;
413  // Clean-clear-sky
414  fluxes_broadband_t clnclrsky_fluxes_sw;
415  clnclrsky_fluxes_sw.flux_up = sw_clnclrsky_flux_up;
416  clnclrsky_fluxes_sw.flux_dn = sw_clnclrsky_flux_dn;
417  clnclrsky_fluxes_sw.flux_dn_dir = sw_clnclrsky_flux_dn_dir;
418  // Clear-sky
419  fluxes_broadband_t clrsky_fluxes_sw;
420  clrsky_fluxes_sw.flux_up = sw_clrsky_flux_up;
421  clrsky_fluxes_sw.flux_dn = sw_clrsky_flux_dn;
422  clrsky_fluxes_sw.flux_dn_dir = sw_clrsky_flux_dn_dir;
423  // Clean-sky
424  fluxes_broadband_t clnsky_fluxes_sw;
425  clnsky_fluxes_sw.flux_up = sw_clnsky_flux_up;
426  clnsky_fluxes_sw.flux_dn = sw_clnsky_flux_dn;
427  clnsky_fluxes_sw.flux_dn_dir = sw_clnsky_flux_dn_dir;
428 
429  // Setup pointers to RRTMGP LW fluxes
430  fluxes_t fluxes_lw;
431  fluxes_lw.flux_up = lw_flux_up;
432  fluxes_lw.flux_dn = lw_flux_dn;
433  fluxes_lw.bnd_flux_up = lw_bnd_flux_up;
434  fluxes_lw.bnd_flux_dn = lw_bnd_flux_dn;
435  // Clean-clear-sky
436  fluxes_broadband_t clnclrsky_fluxes_lw;
437  clnclrsky_fluxes_lw.flux_up = lw_clnclrsky_flux_up;
438  clnclrsky_fluxes_lw.flux_dn = lw_clnclrsky_flux_dn;
439  // Clear-sky
440  fluxes_broadband_t clrsky_fluxes_lw;
441  clrsky_fluxes_lw.flux_up = lw_clrsky_flux_up;
442  clrsky_fluxes_lw.flux_dn = lw_clrsky_flux_dn;
443  // Clean-sky
444  fluxes_broadband_t clnsky_fluxes_lw;
445  clnsky_fluxes_lw.flux_up = lw_clnsky_flux_up;
446  clnsky_fluxes_lw.flux_dn = lw_clnsky_flux_dn;
447 
448  auto nswbands = k_dist_sw_k->get_nband();
449  auto nlwbands = k_dist_lw_k->get_nband();
450 
451  // Setup aerosol optical properties
452  optical_props2_t aerosol_sw;
453  optical_props1_t aerosol_lw;
454  aerosol_sw.init(k_dist_sw_k->get_band_lims_wavenumber());
455  aerosol_sw.alloc_2str(ncol, nlay);
456  aerosol_lw.init(k_dist_lw_k->get_band_lims_wavenumber());
457  aerosol_lw.alloc_1scl(ncol, nlay);
458 
459  // Convert cloud physical properties to optical properties for input to RRTMGP
460  optical_props2_t clouds_sw = get_cloud_optics_sw(ncol, nlay,
462  lwp, iwp, rel, rei);
463  optical_props1_t clouds_lw = get_cloud_optics_lw(ncol, nlay,
465  lwp, iwp, rel, rei);
466  // Do subcolumn sampling to map bands -> gpoints based on cloud fraction and overlap assumption;
467  // This implements the Monte Carlo Independent Column Approximation by mapping only a single
468  // subcolumn (cloud state) to each gpoint.
469  auto nswgpts = k_dist_sw_k->get_ngpt();
470  auto clouds_sw_gpt = get_subsampled_clouds(ncol, nlay, nswbands, nswgpts,
471  clouds_sw, *k_dist_sw_k, cldfrac, p_lay);
472  // Longwave
473  auto nlwgpts = k_dist_lw_k->get_ngpt();
474  auto clouds_lw_gpt = get_subsampled_clouds(ncol, nlay, nlwbands, nlwgpts,
475  clouds_lw, *k_dist_lw_k, cldfrac, p_lay);
476 
477  // Do shortwave
478  rrtmgp_sw(ncol, nlay,
479  *k_dist_sw_k, p_lay, t_lay, p_lev, t_lev, gas_concs,
480  sfc_alb_dir, sfc_alb_dif, mu0, aerosol_sw, clouds_sw_gpt,
481  fluxes_sw, clnclrsky_fluxes_sw, clrsky_fluxes_sw, clnsky_fluxes_sw,
482  tsi_scaling, extra_clnclrsky_diag, extra_clnsky_diag);
483 
484  // Do longwave
485  rrtmgp_lw(ncol, nlay,
486  *k_dist_lw_k, p_lay, t_lay, p_lev, t_lev,
487  t_sfc, sfc_emis, lw_src,
488  gas_concs, aerosol_lw, clouds_lw_gpt,
489  fluxes_lw, clnclrsky_fluxes_lw, clrsky_fluxes_lw, clnsky_fluxes_lw,
490  extra_clnclrsky_diag, extra_clnsky_diag);
491 
492 }
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:32
@ lw_flux_dn
Definition: ERF_NOAHMP.H:37
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:569
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, real1d_k &sfc_emis, 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:875

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