From d26f1c80deb3b282ebfc87246cb0294761f04347 Mon Sep 17 00:00:00 2001 From: Menno Date: Wed, 4 Mar 2026 14:50:42 +0100 Subject: [PATCH 01/11] recategorize .ini switches of the non-raytracing executables --- src_test/test_rte_rrtmgp.cpp | 32 ++++++++++++++++++-------------- src_test/test_rte_rrtmgp.cu | 31 +++++++++++++++++-------------- 2 files changed, 35 insertions(+), 28 deletions(-) diff --git a/src_test/test_rte_rrtmgp.cpp b/src_test/test_rte_rrtmgp.cpp index 2f8a1aa..dfb7c5b 100644 --- a/src_test/test_rte_rrtmgp.cpp +++ b/src_test/test_rte_rrtmgp.cpp @@ -137,14 +137,18 @@ void solve_radiation(int argc, char** argv) const auto settings = toml::parse(case_name + ".ini"); ////// FLOW CONTROL SWITCHES ////// - const bool switch_shortwave = get_ini_value(settings, "switches", "shortwave", true); - const bool switch_longwave = get_ini_value(settings, "switches", "longwave", true); - const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); - const bool switch_liquid_cloud_optics = get_ini_value(settings, "switches", "liquid_cloud_optics", false); - const bool switch_ice_cloud_optics = get_ini_value(settings, "switches", "ice_cloud_optics", false); - const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol_optics", false); - const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta_cloud", true); - const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta_aerosol", false); + const bool switch_liquid_cloud_optics = get_ini_value(settings, "clouds", "liquid_cloud_optics", false); + const bool switch_ice_cloud_optics = get_ini_value(settings, "clouds", "ice_cloud_optics", false); + const bool switch_delta_cloud = get_ini_value(settings, "clouds", "delta_cloud", false); + + const bool switch_shortwave = get_ini_value(settings, "shortwave", "shortwave", true); + const bool switch_sw_fluxes = get_ini_value(settings, "shortwave", "plane_parallel", true); + + const bool switch_longwave = get_ini_value(settings, "longwave", "longwave", true); + const bool switch_lw_fluxes = get_ini_value(settings, "longwave", "plane_parallel", true); + + const bool switch_aerosol_optics = get_ini_value(settings, "aerosols", "aerosol_optics", false); + const bool switch_delta_aerosol = get_ini_value(settings, "aerosols", "delta_aerosol", false); ////// READ THE ATMOSPHERIC DATA ////// @@ -289,7 +293,7 @@ void solve_radiation(int argc, char** argv) Array lw_flux_dn; Array lw_flux_net; - if (switch_fluxes) + if (switch_lw_fluxes) { lw_flux_up .set_dims({n_col, n_lev}); lw_flux_dn .set_dims({n_col, n_lev}); @@ -302,7 +306,7 @@ void solve_radiation(int argc, char** argv) auto time_start = std::chrono::high_resolution_clock::now(); rad_lw.solve( - switch_fluxes, + switch_lw_fluxes, switch_cloud_optics, gas_concs, p_lay, p_lev, @@ -328,7 +332,7 @@ void solve_radiation(int argc, char** argv) auto nc_lw_band_lims_wvn = output_nc.add_variable("lw_band_lims_wvn", {"band_lw", "pair"}); nc_lw_band_lims_wvn.insert(rad_lw.get_band_lims_wavenumber().v(), {0, 0}); - if (switch_fluxes) + if (switch_lw_fluxes) { auto nc_lw_flux_up = output_nc.add_variable("lw_flux_up" , {"lev", "y", "x"}); auto nc_lw_flux_dn = output_nc.add_variable("lw_flux_dn" , {"lev", "y", "x"}); @@ -383,7 +387,7 @@ void solve_radiation(int argc, char** argv) Array sw_flux_dn_dir; Array sw_flux_net; - if (switch_fluxes) + if (switch_sw_fluxes) { sw_flux_up .set_dims({n_col, n_lev}); sw_flux_dn .set_dims({n_col, n_lev}); @@ -397,7 +401,7 @@ void solve_radiation(int argc, char** argv) auto time_start = std::chrono::high_resolution_clock::now(); rad_sw.solve( - switch_fluxes, + switch_sw_fluxes, switch_cloud_optics, switch_aerosol_optics, switch_delta_cloud, @@ -430,7 +434,7 @@ void solve_radiation(int argc, char** argv) auto nc_sw_band_lims_wvn = output_nc.add_variable("sw_band_lims_wvn", {"band_sw", "pair"}); nc_sw_band_lims_wvn.insert(rad_sw.get_band_lims_wavenumber().v(), {0, 0}); - if (switch_fluxes) + if (switch_sw_fluxes) { auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); auto nc_sw_flux_dn = output_nc.add_variable("sw_flux_dn" , {"lev", "y", "x"}); diff --git a/src_test/test_rte_rrtmgp.cu b/src_test/test_rte_rrtmgp.cu index 0bb86ca..57c1aca 100644 --- a/src_test/test_rte_rrtmgp.cu +++ b/src_test/test_rte_rrtmgp.cu @@ -163,15 +163,18 @@ void solve_radiation(int argc, char** argv) const auto settings = toml::parse(case_name + ".ini"); ////// FLOW CONTROL SWITCHES ////// - const bool switch_shortwave = get_ini_value(settings, "switches", "shortwave", true); - const bool switch_longwave = get_ini_value(settings, "switches", "longwave", true); - const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); - const bool switch_liquid_cloud_optics = get_ini_value(settings, "switches", "liquid_cloud_optics", false); - const bool switch_ice_cloud_optics = get_ini_value(settings, "switches", "ice_cloud_optics", false); - const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol_optics", false); - const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta_cloud", true); - const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta_aerosol", false); + const bool switch_liquid_cloud_optics = get_ini_value(settings, "clouds", "liquid_cloud_optics", false); + const bool switch_ice_cloud_optics = get_ini_value(settings, "clouds", "ice_cloud_optics", false); + const bool switch_delta_cloud = get_ini_value(settings, "clouds", "delta_cloud", false); + const bool switch_shortwave = get_ini_value(settings, "shortwave", "shortwave", true); + const bool switch_sw_fluxes = get_ini_value(settings, "shortwave", "plane_parallel", true); + + const bool switch_longwave = get_ini_value(settings, "longwave", "longwave", true); + const bool switch_lw_fluxes = get_ini_value(settings, "longwave", "plane_parallel", true); + + const bool switch_aerosol_optics = get_ini_value(settings, "aerosols", "aerosol_optics", false); + const bool switch_delta_aerosol = get_ini_value(settings, "aerosols", "delta_aerosol", false); ////// READ THE ATMOSPHERIC DATA ////// Status::print_message("Reading atmospheric input data from NetCDF."); @@ -336,7 +339,7 @@ void solve_radiation(int argc, char** argv) Array_gpu lw_flux_dn; Array_gpu lw_flux_net; - if (switch_fluxes) + if (switch_lw_fluxes) { lw_flux_up .set_dims({n_col, n_lev}); lw_flux_dn .set_dims({n_col, n_lev}); @@ -370,7 +373,7 @@ void solve_radiation(int argc, char** argv) cudaEventRecord(start, 0); rad_lw.solve_gpu( - switch_fluxes, + switch_lw_fluxes, switch_cloud_optics, gas_concs_gpu, p_lay_gpu, p_lev_gpu, @@ -407,7 +410,7 @@ void solve_radiation(int argc, char** argv) auto nc_lw_band_lims_wvn = output_nc.add_variable("lw_band_lims_wvn", {"band_lw", "pair"}); nc_lw_band_lims_wvn.insert(rad_lw.get_band_lims_wavenumber_gpu().v(), {0, 0}); - if (switch_fluxes) + if (switch_lw_fluxes) { auto nc_lw_flux_up = output_nc.add_variable("lw_flux_up" , {"lev", "y", "x"}); auto nc_lw_flux_dn = output_nc.add_variable("lw_flux_dn" , {"lev", "y", "x"}); @@ -464,7 +467,7 @@ void solve_radiation(int argc, char** argv) Array_gpu sw_flux_dn_dir; Array_gpu sw_flux_net; - if (switch_fluxes) + if (switch_sw_fluxes) { sw_flux_up .set_dims({n_col, n_lev}); sw_flux_dn .set_dims({n_col, n_lev}); @@ -503,7 +506,7 @@ void solve_radiation(int argc, char** argv) cudaEventRecord(start, 0); rad_sw.solve_gpu( - switch_fluxes, + switch_sw_fluxes, switch_cloud_optics, switch_aerosol_optics, switch_delta_cloud, @@ -548,7 +551,7 @@ void solve_radiation(int argc, char** argv) auto nc_sw_band_lims_wvn = output_nc.add_variable("sw_band_lims_wvn", {"band_sw", "pair"}); nc_sw_band_lims_wvn.insert(rad_sw.get_band_lims_wavenumber_gpu().v(), {0, 0}); - if (switch_fluxes) + if (switch_sw_fluxes) { auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); auto nc_sw_flux_dn = output_nc.add_variable("sw_flux_dn" , {"lev", "y", "x"}); From e402bfc6e25df146d04e4ff37bc0cc155c6f6f1d Mon Sep 17 00:00:00 2001 From: Menno Date: Wed, 4 Mar 2026 16:47:29 +0100 Subject: [PATCH 02/11] large makeover of .ini inputs --- include_test/radiation_solver_rt.h | 7 +- src_test/radiation_solver_rt.cu | 296 +++++++++++----------- src_test/test_rte_rrtmgp_rt.cu | 391 ++++++++++++++--------------- 3 files changed, 338 insertions(+), 356 deletions(-) diff --git a/include_test/radiation_solver_rt.h b/include_test/radiation_solver_rt.h index 9df1af2..2271cda 100644 --- a/include_test/radiation_solver_rt.h +++ b/include_test/radiation_solver_rt.h @@ -46,13 +46,12 @@ class Radiation_solver_longwave #ifdef __CUDACC__ void solve_gpu( - const bool switch_fluxes, const bool switch_raytracing, + const bool switch_plane_parallel, const bool switch_cloud_optics, const bool switch_aerosol_optics, const bool switch_delta_cloud, const bool switch_delta_aerosol, - const bool switch_single_gpt, const bool switch_lw_scattering, const bool switch_independent_column, const int single_gpt, @@ -122,14 +121,12 @@ class Radiation_solver_shortwave #ifdef __CUDACC__ void solve_gpu( - const bool switch_fluxes, - const bool switch_twostream, const bool switch_raytracing, + const bool switch_plane_parallel, const bool switch_independent_column, const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, - const bool switch_single_gpt, const bool switch_delta_cloud, const bool switch_delta_aerosol, const bool switch_attenuate_tica, diff --git a/src_test/radiation_solver_rt.cu b/src_test/radiation_solver_rt.cu index 7278181..d4540af 100644 --- a/src_test/radiation_solver_rt.cu +++ b/src_test/radiation_solver_rt.cu @@ -479,13 +479,12 @@ Radiation_solver_longwave::Radiation_solver_longwave( void Radiation_solver_longwave::solve_gpu( - const bool switch_fluxes, + const bool switch_plane_parallel, const bool switch_raytracing, const bool switch_cloud_optics, const bool switch_aerosol_optics, const bool switch_delta_cloud, const bool switch_delta_aerosol, - const bool switch_single_gpt, const bool switch_lw_scattering, const bool switch_independent_column, const int single_gpt, @@ -533,15 +532,13 @@ void Radiation_solver_longwave::solve_gpu( Gas_optics_rrtmgp_rt::get_col_dry(col_dry, gas_concs.get_vmr("h2o"), p_lev); } - if (switch_fluxes) + if (switch_plane_parallel) { Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, n_col, lw_flux_up.ptr()); Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, n_col, lw_flux_dn.ptr()); Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, n_col, lw_flux_net.ptr()); if (switch_raytracing) { } - - } const Array& band_limits_gpt(this->kdist_gpu->get_band_lims_gpoint()); @@ -692,7 +689,7 @@ void Radiation_solver_longwave::solve_gpu( } // Store the optical properties, if desired. - if (switch_single_gpt && igpt == single_gpt) + if (igpt == single_gpt) { lay_source = (*sources).get_lay_source(); lev_source = (*sources).get_lev_source(); @@ -713,13 +710,15 @@ void Radiation_solver_longwave::solve_gpu( } - if (switch_fluxes) - { - constexpr int n_ang = 1; + constexpr int n_ang = 1; + + const bool use_raytracer = switch_raytracing && (lowest_gas_mean_free_path / grid_d_xy_min) > min_mfp_grid_ratio; - std::unique_ptr fluxes = - std::make_unique(grid_cells.x, grid_cells.y, grid_cells.z, n_lev); + std::unique_ptr fluxes = + std::make_unique(grid_cells.x, grid_cells.y, grid_cells.z, n_lev); + if (switch_plane_parallel || !use_raytracer) + { rte_lw.rte_lw( optical_props, top_at_1, @@ -733,78 +732,80 @@ void Radiation_solver_longwave::solve_gpu( (*fluxes).net_flux(); // Copy the data to the output. - Gpt_combine_kernels_cuda_rt::add_from_gpoint( + if (switch_plane_parallel) + { + Gpt_combine_kernels_cuda_rt::add_from_gpoint( n_col, n_lev, lw_flux_up.ptr(), lw_flux_dn.ptr(), lw_flux_net.ptr(), (*fluxes).get_flux_up().ptr(), (*fluxes).get_flux_dn().ptr(), (*fluxes).get_flux_net().ptr()); + } - if (switch_single_gpt && igpt == single_gpt) + if (igpt == single_gpt) { lw_gpt_flux_up = (*fluxes).get_flux_up(); lw_gpt_flux_dn = (*fluxes).get_flux_dn(); lw_gpt_flux_net = (*fluxes).get_flux_net(); } + } - if (switch_raytracing) + if (switch_raytracing) + { + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_tod_dn().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_tod_up().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_sfc_dif().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_sfc_up().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, grid_cells.z, (*fluxes).get_flux_abs_dif().ptr()); + + if (use_raytracer) { - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_tod_dn().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_tod_up().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_sfc_dif().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_sfc_up().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, grid_cells.z, (*fluxes).get_flux_abs_dif().ptr()); + ++monte_carlo_gpoints; - if ( (lowest_gas_mean_free_path / grid_d_xy_min) > min_mfp_grid_ratio) - { - ++monte_carlo_gpoints; - - raytracer_lw.trace_rays( - igpt, - switch_independent_column, - ray_count, - grid_cells, - grid_d, - kn_grid, - dynamic_cast(*optical_props).get_tau(), - dynamic_cast(*optical_props).get_ssa(), - dynamic_cast(*cloud_optical_props).get_tau(), - dynamic_cast(*cloud_optical_props).get_ssa(), - dynamic_cast(*cloud_optical_props).get_g(), - dynamic_cast(*aerosol_optical_props).get_tau(), - dynamic_cast(*aerosol_optical_props).get_ssa(), - dynamic_cast(*aerosol_optical_props).get_g(), - (*sources).get_lay_source(), - (*sources).get_sfc_source(), - emis_sfc, - (*fluxes).get_flux_dn()({1, grid_cells.z}), - (*fluxes).get_flux_tod_dn(), - (*fluxes).get_flux_tod_up(), - (*fluxes).get_flux_sfc_dif(), - (*fluxes).get_flux_sfc_up(), - (*fluxes).get_flux_abs_dif()); - } - else - { - convert_1d_to_rt_output( - n_col, n_lay, grid_cells.z, grid_d.z, - (*fluxes).get_flux_up(), - (*fluxes).get_flux_dn(), - (*fluxes).get_flux_net(), + raytracer_lw.trace_rays( + igpt, + switch_independent_column, + ray_count, + grid_cells, + grid_d, + kn_grid, + dynamic_cast(*optical_props).get_tau(), + dynamic_cast(*optical_props).get_ssa(), + dynamic_cast(*cloud_optical_props).get_tau(), + dynamic_cast(*cloud_optical_props).get_ssa(), + dynamic_cast(*cloud_optical_props).get_g(), + dynamic_cast(*aerosol_optical_props).get_tau(), + dynamic_cast(*aerosol_optical_props).get_ssa(), + dynamic_cast(*aerosol_optical_props).get_g(), + (*sources).get_lay_source(), + (*sources).get_sfc_source(), + emis_sfc, + (*fluxes).get_flux_dn()({1, grid_cells.z}), (*fluxes).get_flux_tod_dn(), (*fluxes).get_flux_tod_up(), (*fluxes).get_flux_sfc_dif(), (*fluxes).get_flux_sfc_up(), (*fluxes).get_flux_abs_dif()); + } + else + { + convert_1d_to_rt_output( + n_col, n_lay, grid_cells.z, grid_d.z, + (*fluxes).get_flux_up(), + (*fluxes).get_flux_dn(), + (*fluxes).get_flux_net(), + (*fluxes).get_flux_tod_dn(), + (*fluxes).get_flux_tod_up(), + (*fluxes).get_flux_sfc_dif(), + (*fluxes).get_flux_sfc_up(), + (*fluxes).get_flux_abs_dif()); - } - - Gpt_combine_kernels_cuda_rt::add_from_gpoint( - grid_cells.x, grid_cells.y, - rt_flux_tod_dn.ptr(), rt_flux_tod_up.ptr(), rt_flux_sfc_dn.ptr(), rt_flux_sfc_up.ptr(), - (*fluxes).get_flux_tod_dn().ptr(), (*fluxes).get_flux_tod_up().ptr(), (*fluxes).get_flux_sfc_dif().ptr(), (*fluxes).get_flux_sfc_up().ptr()); + } - Gpt_combine_kernels_cuda_rt::add_from_gpoint( - n_col, grid_cells.z, rt_flux_abs.ptr(), (*fluxes).get_flux_abs_dif().ptr()); + Gpt_combine_kernels_cuda_rt::add_from_gpoint( + grid_cells.x, grid_cells.y, + rt_flux_tod_dn.ptr(), rt_flux_tod_up.ptr(), rt_flux_sfc_dn.ptr(), rt_flux_sfc_up.ptr(), + (*fluxes).get_flux_tod_dn().ptr(), (*fluxes).get_flux_tod_up().ptr(), (*fluxes).get_flux_sfc_dif().ptr(), (*fluxes).get_flux_sfc_up().ptr()); - } + Gpt_combine_kernels_cuda_rt::add_from_gpoint( + n_col, grid_cells.z, rt_flux_abs.ptr(), (*fluxes).get_flux_abs_dif().ptr()); } Tools_gpu::free_gpu(max_tau_gas_g); @@ -846,14 +847,12 @@ void Radiation_solver_shortwave::load_mie_tables( } void Radiation_solver_shortwave::solve_gpu( - const bool switch_fluxes, - const bool switch_twostream, + const bool switch_plane_parallel, const bool switch_raytracing, const bool switch_independent_column, const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, - const bool switch_single_gpt, const bool switch_delta_cloud, const bool switch_delta_aerosol, const bool switch_attenuate_tica, @@ -914,24 +913,22 @@ void Radiation_solver_shortwave::solve_gpu( Array_gpu mie_cdfs_sub; Array_gpu mie_angs_sub; - if (switch_fluxes) + if (switch_plane_parallel) { - if (switch_twostream) - { - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_up.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn_dir.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_net.ptr()); - } - if (switch_raytracing) - { - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_tod_up.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_dir.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_dif.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_up.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.z, grid_cells.y, grid_cells.x, rt_flux_abs_dir.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.z, grid_cells.y, grid_cells.x, rt_flux_abs_dif.ptr()); - } + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_up.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn_dir.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_net.ptr()); + } + + if (switch_raytracing) + { + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_tod_up.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_dir.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_dif.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_up.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.z, grid_cells.y, grid_cells.x, rt_flux_abs_dir.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.z, grid_cells.y, grid_cells.x, rt_flux_abs_dif.ptr()); } const Array& band_limits_gpt(this->kdist_gpu->get_band_lims_gpoint()); @@ -1064,7 +1061,7 @@ void Radiation_solver_shortwave::solve_gpu( } // Store the optical properties, if desired - if (switch_single_gpt && igpt == single_gpt) + if (igpt == single_gpt) { tot_tau_out = optical_props->get_tau(); tot_ssa_out = optical_props->get_ssa(); @@ -1076,77 +1073,75 @@ void Radiation_solver_shortwave::solve_gpu( aer_asy_out = aerosol_optical_props->get_g(); } - if (switch_fluxes) + std::unique_ptr fluxes = + std::make_unique(grid_cells.x, grid_cells.y, grid_cells.z, n_lev); + + if (switch_plane_parallel) { - std::unique_ptr fluxes = - std::make_unique(grid_cells.x, grid_cells.y, grid_cells.z, n_lev); - if (switch_twostream) - { + rte_sw.rte_sw( + optical_props, + top_at_1, + mu0, + toa_src, + sfc_alb_dir.subset({{ {band, band}, {1, n_col}}}), + sfc_alb_dif.subset({{ {band, band}, {1, n_col}}}), + Array_gpu(), // Add an empty array, no inc_flux. + (*fluxes).get_flux_up(), + (*fluxes).get_flux_dn(), + (*fluxes).get_flux_dn_dir()); + } - rte_sw.rte_sw( - optical_props, - top_at_1, - mu0, - toa_src, - sfc_alb_dir.subset({{ {band, band}, {1, n_col}}}), - sfc_alb_dif.subset({{ {band, band}, {1, n_col}}}), - Array_gpu(), // Add an empty array, no inc_flux. - (*fluxes).get_flux_up(), - (*fluxes).get_flux_dn(), - (*fluxes).get_flux_dn_dir()); - } + if (switch_raytracing) + { + Float zenith_angle = std::acos(mu0({1})); + Float azimuth_angle = azi({1}); // sun approximately from south - if (switch_raytracing) + if (switch_cloud_mie) { - Float zenith_angle = std::acos(mu0({1})); - Float azimuth_angle = azi({1}); // sun approximately from south - - if (switch_cloud_mie) - { - mie_cdfs_sub = mie_cdfs.subset({{ {1, n_mie}, {band, band} }}); - mie_angs_sub = mie_angs.subset({{ {1, n_mie}, {1, n_re}, {band, band} }}); - } - - raytracer.trace_rays( - igpt, - switch_independent_column, - ray_count, - grid_cells, - grid_d, - kn_grid, - mie_cdfs_sub, - mie_angs_sub, - dynamic_cast(*optical_props).get_tau(), - dynamic_cast(*optical_props).get_ssa(), - dynamic_cast(*cloud_optical_props).get_tau(), - dynamic_cast(*cloud_optical_props).get_ssa(), - dynamic_cast(*cloud_optical_props).get_g(), - dynamic_cast(*aerosol_optical_props).get_tau(), - dynamic_cast(*aerosol_optical_props).get_ssa(), - dynamic_cast(*aerosol_optical_props).get_g(), - rel, - sfc_alb_dir.subset({{ {band, band}, {1, n_col}}}), - zenith_angle, - azimuth_angle, - toa_src({1}) * mu0({1}), - Float(0.), - (*fluxes).get_flux_tod_dn(), - (*fluxes).get_flux_tod_up(), - (*fluxes).get_flux_sfc_dir(), - (*fluxes).get_flux_sfc_dif(), - (*fluxes).get_flux_sfc_up(), - (*fluxes).get_flux_abs_dir(), - (*fluxes).get_flux_abs_dif()); + mie_cdfs_sub = mie_cdfs.subset({{ {1, n_mie}, {band, band} }}); + mie_angs_sub = mie_angs.subset({{ {1, n_mie}, {1, n_re}, {band, band} }}); } - (*fluxes).net_flux(); - if (switch_twostream) - { - Gpt_combine_kernels_cuda_rt::add_from_gpoint( - n_col, n_lev, sw_flux_up.ptr(), sw_flux_dn.ptr(), sw_flux_dn_dir.ptr(), sw_flux_net.ptr(), - (*fluxes).get_flux_up().ptr(), (*fluxes).get_flux_dn().ptr(), (*fluxes).get_flux_dn_dir().ptr(), (*fluxes).get_flux_net().ptr()); - } + raytracer.trace_rays( + igpt, + switch_independent_column, + ray_count, + grid_cells, + grid_d, + kn_grid, + mie_cdfs_sub, + mie_angs_sub, + dynamic_cast(*optical_props).get_tau(), + dynamic_cast(*optical_props).get_ssa(), + dynamic_cast(*cloud_optical_props).get_tau(), + dynamic_cast(*cloud_optical_props).get_ssa(), + dynamic_cast(*cloud_optical_props).get_g(), + dynamic_cast(*aerosol_optical_props).get_tau(), + dynamic_cast(*aerosol_optical_props).get_ssa(), + dynamic_cast(*aerosol_optical_props).get_g(), + rel, + sfc_alb_dir.subset({{ {band, band}, {1, n_col}}}), + zenith_angle, + azimuth_angle, + toa_src({1}) * mu0({1}), + Float(0.), + (*fluxes).get_flux_tod_dn(), + (*fluxes).get_flux_tod_up(), + (*fluxes).get_flux_sfc_dir(), + (*fluxes).get_flux_sfc_dif(), + (*fluxes).get_flux_sfc_up(), + (*fluxes).get_flux_abs_dir(), + (*fluxes).get_flux_abs_dif()); + } + + (*fluxes).net_flux(); + if (switch_plane_parallel) + { + Gpt_combine_kernels_cuda_rt::add_from_gpoint( + n_col, n_lev, sw_flux_up.ptr(), sw_flux_dn.ptr(), sw_flux_dn_dir.ptr(), sw_flux_net.ptr(), + (*fluxes).get_flux_up().ptr(), (*fluxes).get_flux_dn().ptr(), (*fluxes).get_flux_dn_dir().ptr(), (*fluxes).get_flux_net().ptr()); + } if (switch_raytracing) { @@ -1159,14 +1154,13 @@ void Radiation_solver_shortwave::solve_gpu( (*fluxes).get_flux_abs_dir().ptr(), (*fluxes).get_flux_abs_dif().ptr()); } - if (switch_single_gpt && igpt == single_gpt) + if (igpt == single_gpt) { sw_gpt_flux_up = (*fluxes).get_flux_up(); sw_gpt_flux_dn = (*fluxes).get_flux_dn(); sw_gpt_flux_dn_dir = (*fluxes).get_flux_dn_dir(); sw_gpt_flux_net = (*fluxes).get_flux_net(); } - } previous_band = band; } } diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index 4dd5911..00ba0f5 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -166,49 +166,75 @@ void solve_radiation(int argc, char** argv) const auto settings = toml::parse(case_name + ".ini"); - ////// FLOW CONTROL SWITCHES ////// - const bool switch_shortwave = get_ini_value(settings, "switches", "shortwave", true); - const bool switch_longwave = get_ini_value(settings, "switches", "longwave", true); - const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); - bool switch_sw_twostream = get_ini_value(settings, "switches", "sw_two_stream", false); - bool switch_sw_raytracing = get_ini_value(settings, "switches", "sw_raytracing", true); - bool switch_lw_raytracing = get_ini_value(settings, "switches", "lw_raytracing", true); - bool switch_independent_column = get_ini_value(settings, "switches", "independent_column", false); - bool switch_liquid_cloud_optics = get_ini_value(settings, "switches", "liquid_cloud_optics", false); - bool switch_ice_cloud_optics = get_ini_value(settings, "switches", "ice_cloud_optics", false); - const bool switch_lw_scattering = get_ini_value(settings, "switches", "lw_scattering", false); - const bool switch_min_mfp_grid_ratio= get_ini_value(settings, "switches", "min_mfp_grid_ratio", true); - const bool switch_cloud_mie = get_ini_value(settings, "switches", "cloud_mie", false); - const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol_optics", false); - const bool switch_single_gpt = get_ini_value(settings, "switches", "single_gpt", false); - const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta_cloud", false); - const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta_aerosol", false); - const bool switch_tica = get_ini_value(settings, "switches", "tica", false); - - const bool switch_bw_raytracing = get_ini_value(settings, "switches", "bw_raytracing", true); - const bool switch_lu_albedo = get_ini_value(settings, "switches", "lu_albedo", false); - const bool switch_image = get_ini_value(settings, "switches", "image", true); - const bool switch_broadband = get_ini_value(settings, "switches", "broadband", false); - const bool switch_cloud_cam = get_ini_value(settings, "switches", "cloud_cam", false); - - const int photons_per_pixel = get_ini_value(settings, "ints", "bw_raytracing", 1); - - if (!switch_shortwave) - switch_sw_raytracing = false; - - if (!switch_longwave) - switch_lw_raytracing = false; - - if (switch_shortwave && !switch_sw_twostream && !switch_sw_raytracing) + ////// READ INI FILE ////// + const bool switch_liquid_cloud_optics = get_ini_value(settings, "clouds", "liquid_cloud_optics", false); + const bool switch_ice_cloud_optics = get_ini_value(settings, "clouds", "ice_cloud_optics", false); + const bool switch_delta_cloud = get_ini_value(settings, "clouds", "delta_cloud", false); + const bool switch_tilted_columns = get_ini_value(settings, "clouds", "tilted_columns", false); + const bool switch_cloud_mie = get_ini_value(settings, "clouds", "cloud_mie", false); // use mie scattering. Curerntly only for shortwave and only for liquid cloud droplets (requires ice_cloud_optics = false) + + const bool switch_aerosol_optics = get_ini_value(settings, "aerosols", "aerosol_optics", false); + const bool switch_delta_aerosol = get_ini_value(settings, "aerosols", "delta_aerosol", false); + + const bool switch_shortwave = get_ini_value(settings, "shortwave", "shortwave", true); + const bool switch_sw_raytracing = get_ini_value(settings, "shortwave", "raytracing", true); // compute and output ray tracer fluxes + const bool switch_sw_plane_parallel = get_ini_value(settings, "shortwave", "plane_parallel", true); // compute and output plane parallel 1D fluxes (two-stream) + bool switch_sw_independent_column = get_ini_value(settings, "shortwave", "independent_column", false); // solve ray tracer in independent column mode + const int sw_single_gpt = get_ini_value(settings, "shortwave", "single_gpt", -1); // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) + + const bool switch_longwave = get_ini_value(settings, "longwave", "longwave", true); + const bool switch_lw_raytracing = get_ini_value(settings, "longwave", "raytracing", true); // compute and output ray tracer fluxes + const bool switch_lw_plane_parallel = get_ini_value(settings, "longwave", "plane_parallel", true); // compute and output plane parallel 1D fluxes (two-stream or no-scattering solution) + const bool switch_lw_scattering = get_ini_value(settings, "longwave", "scattering", false); // enable scattering in longwave solver (only possible in ray tracer executable) + bool switch_lw_independent_column = get_ini_value(settings, "longwave", "independent_column", false); // solve ray tracer in independent column mode + const int lw_single_gpt = get_ini_value(settings, "longwave", "single_gpt", -1); // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) + const Float min_mfp_grid_ratio = get_ini_value(settings, "longwave", "min_mfp_grid_ratio", Float(1.0)); // minimum ratio between the lowest gasous mean fee path and the horizontal grid spacing at which ray tracer is still used. Set to 0. to use ray tracer for all g-points + + const bool switch_bw_raytracing = get_ini_value(settings, "backward", "bw_raytracing", true); + const bool switch_lu_albedo = get_ini_value(settings, "backward", "lu_albedo", false); // read surface type from land_use_map variabele and compute spectral albedo and reflection type (lambertian/specular) accordingly + const bool switch_image = get_ini_value(settings, "backward", "image", true); // solve visible bands and convert to XYZ tristimulus values + const bool switch_broadband = get_ini_value(settings, "backward", "broadband", false); // solve broadband radiances + const bool switch_cloud_cam = get_ini_value(settings, "backward", "cloud_cam", false); // output additional cloud statistics for each camera pixel + + const Float input_sza = get_ini_value(settings, "solar_angles", "sza", -1.0); // if >0, overwrite zenith angle from input netcdf file + const Float input_azi = get_ini_value(settings, "solar_angles", "azi", -1.0); // if >0, overwrite azimuth angle from input netcdf file + + Camera camera; + if (switch_bw_raytracing) { - std::string error = "With shortwave enable, need to run either two-stream solver or ray tracer "; + camera.fov = get_ini_value(settings, "camera", "cam_field-of-view", 80.); + camera.cam_type = get_ini_value(settings, "camera", "cam_type", 0); + camera.position = {get_ini_value(settings, "camera", "cam_px", 0.), + get_ini_value(settings, "camera", "cam_py", 0.), + get_ini_value(settings, "camera", "cam_pz", 100.)}; + camera.nx = get_ini_value(settings, "camera", "cam_nx", 0); + camera.ny = get_ini_value(settings, "camera", "cam_ny", 0); + camera.setup_rotation_matrix(get_ini_value(settings, "camera", "cam_yaw", 0.), + get_ini_value(settings, "camera", "cam_pitch", 0.), + get_ini_value(settings, "camera", "cam_roll", 0.)); + camera.setup_normal_camera(camera); + + camera.npix = Int(camera.nx * camera.ny); + } + + + if (switch_shortwave && !(switch_sw_plane_parallel || switch_sw_raytracing)) + { + std::string error = "With shortwave enabled, need to run either shortwave plane parallel solver or ray tracer "; throw std::runtime_error(error); } + if (switch_longwave && !(switch_lw_plane_parallel || switch_lw_raytracing)) + { + std::string error = "With longwave enabled, need to run either longwave plane parallel solver or ray tracer "; + throw std::runtime_error(error); + } + + // read samples counts (if applicable) Int sw_photons_per_pixel; if (switch_sw_raytracing) { - sw_photons_per_pixel = get_ini_value(settings, "ints", "sw_raytracing", Int(256)); + sw_photons_per_pixel = get_ini_value(settings, "shortwave", "samples", Int(256)); if (Float(int(std::log2(Float(sw_photons_per_pixel)))) != std::log2(Float(sw_photons_per_pixel))) { std::string error = "number of photons per pixel should be a power of 2 "; @@ -216,12 +242,15 @@ void solve_radiation(int argc, char** argv) } } + Int bw_photons_per_pixel; + if (switch_bw_raytracing) + bw_photons_per_pixel = get_ini_value(settings, "backward", "samples", Int(1)); Int lw_photon_power; Int lw_photon_count; if (switch_lw_raytracing) { - lw_photon_power = get_ini_value(settings, "ints", "lw_raytracing", Int(22)); + lw_photon_power = get_ini_value(settings, "longwave", "samples", Int(22)); lw_photon_count = 1 << lw_photon_power; } @@ -230,8 +259,11 @@ void solve_radiation(int argc, char** argv) Status::print_warning("No longwave radiation implemented in the backward ray tracer"); } - if (switch_tica) - switch_independent_column = true; + if (switch_tilted_columns) + { + switch_sw_independent_column = true; + switch_lw_independent_column = true; + } if (switch_cloud_cam && !(switch_liquid_cloud_optics || switch_ice_cloud_optics)) throw std::runtime_error("Enable cloud optics (liquid and/or ice) when cloud-cam is switch on!"); @@ -242,10 +274,6 @@ void solve_radiation(int argc, char** argv) throw std::runtime_error(error); } - int single_gpt = get_ini_value(settings, "ints", "single-gpt", 1); - - const Float min_mfp_grid_ratio = switch_min_mfp_grid_ratio ? get_ini_value(settings, "floats", "min_mfp_grid_ratio", Float(1.)) : Float(0.); - if (switch_sw_raytracing) Status::print_message("Shortwave: using "+ std::to_string(sw_photons_per_pixel) + " rays per pixel per g-point"); @@ -253,9 +281,6 @@ void solve_radiation(int argc, char** argv) Status::print_message("Longwave: using 2**"+std::to_string(lw_photon_power) + " ("+std::to_string(lw_photon_count) + ") rays per g-point"); - const Float input_sza = get_ini_value(settings, "floats", "sza", -1.0); - const Float input_azi = get_ini_value(settings, "floats", "azi", -1.0); - ////// READ THE ATMOSPHERIC DATA ////// Status::print_message("Reading atmospheric input data from NetCDF."); @@ -287,26 +312,6 @@ void solve_radiation(int argc, char** argv) input_nc.get_variable("ngrid_y"), input_nc.get_variable("ngrid_z")}; - Camera camera; - if (switch_bw_raytracing) - { - camera.fov = get_ini_value(settings, "bw_camera", "field-of-view", 80.); - camera.cam_type = get_ini_value(settings, "bw_camera", "cam-type", 0); - camera.position = {get_ini_value(settings, "bw_camera", "px", 0.), - get_ini_value(settings, "bw_camera", "py", 0.), - get_ini_value(settings, "bw_camera", "pz", 100.)}; - camera.nx = get_ini_value(settings, "bw_camera", "nx", 0); - camera.ny = get_ini_value(settings, "bw_camera", "ny", 0); - camera.setup_rotation_matrix(get_ini_value(settings, "bw_camera", "yaw", 0.), - get_ini_value(settings, "bw_camera", "pitch", 0.), - get_ini_value(settings, "bw_camera", "roll", 0.)); - camera.setup_normal_camera(camera); - - camera.npix = Int(camera.nx * camera.ny); - - const Float input_sza = get_ini_value(settings, "floats", "sza", -1.0); - const Float input_azi = get_ini_value(settings, "floats", "azi", -1.0); - } // Read the atmospheric fields. Array p_lay(input_nc.get_variable("p_lay", {n_lay, n_col_y, n_col_x}), {n_col, n_lay}); @@ -314,14 +319,14 @@ void solve_radiation(int argc, char** argv) Array p_lev(input_nc.get_variable("p_lev", {n_lev, n_col_y, n_col_x}), {n_col, n_lev}); Array t_lev(input_nc.get_variable("t_lev", {n_lev, n_col_y, n_col_x}), {n_col, n_lev}); - if (input_nc.variable_exists("col_dry") && switch_tica) + if (input_nc.variable_exists("col_dry") && switch_tilted_columns) { std::string error = "col_dry is not supported in tica mode"; throw std::runtime_error(error); } - if (input_nc.variable_exists("tsi") && switch_tica) + if (input_nc.variable_exists("tsi") && switch_tilted_columns) { std::string error = "tsi is overwritten in tica mode"; throw std::runtime_error(error); @@ -435,7 +440,7 @@ void solve_radiation(int argc, char** argv) azi.fill(input_azi / Float(180.0) * M_PI); - if (switch_tica) + if (switch_tilted_columns) { tica_sza = acos(mu0.v()[0]); tica_azi = azi.v()[0]; @@ -628,7 +633,7 @@ void solve_radiation(int argc, char** argv) Array_gpu lev_source; Array_gpu sfc_source; - if (switch_single_gpt) + if (lw_single_gpt > 0) { lw_tau_tot .set_dims({n_col, n_lay}); lw_tau_cld .set_dims({n_col, n_lay}); @@ -650,7 +655,7 @@ void solve_radiation(int argc, char** argv) Array_gpu lw_flux_dn; Array_gpu lw_flux_net; - if (switch_fluxes) + if (switch_lw_plane_parallel) { lw_flux_up .set_dims({n_col, n_lev}); lw_flux_dn .set_dims({n_col, n_lev}); @@ -661,7 +666,7 @@ void solve_radiation(int argc, char** argv) Array_gpu lw_gpt_flux_dn; Array_gpu lw_gpt_flux_net; - if (switch_single_gpt) + if (lw_single_gpt > 0) { lw_gpt_flux_up .set_dims({n_col, n_lev}); lw_gpt_flux_dn .set_dims({n_col, n_lev}); @@ -718,16 +723,15 @@ void solve_radiation(int argc, char** argv) cudaEventRecord(start, 0); rad_lw.solve_gpu( - switch_fluxes, switch_lw_raytracing, + switch_lw_plane_parallel, switch_cloud_optics, switch_aerosol_optics, switch_delta_cloud, switch_delta_aerosol, - switch_single_gpt, switch_lw_scattering, - switch_independent_column, - single_gpt, + switch_lw_independent_column, + lw_single_gpt, min_mfp_grid_ratio, lw_photon_count, grid_cells, @@ -801,7 +805,7 @@ void solve_radiation(int argc, char** argv) auto nc_lw_band_lims_wvn = output_nc.add_variable("lw_band_lims_wvn", {"band_lw", "pair"}); nc_lw_band_lims_wvn.insert(rad_lw.get_band_lims_wavenumber_gpu().v(), {0, 0}); - if (switch_single_gpt) + if (lw_single_gpt > 0) { auto nc_lw_band_lims_gpt = output_nc.add_variable("lw_band_lims_gpt", {"band_lw", "pair"}); nc_lw_band_lims_gpt.insert(rad_lw.get_band_lims_gpoint_gpu().v(), {0, 0}); @@ -844,7 +848,7 @@ void solve_radiation(int argc, char** argv) nc_sfc_source.insert(sfc_source_cpu.v(), {0, 0}); } - if (switch_fluxes) + if (switch_lw_plane_parallel) { auto nc_lw_flux_up = output_nc.add_variable("lw_flux_up" , {"lev", "y", "x"}); auto nc_lw_flux_dn = output_nc.add_variable("lw_flux_dn" , {"lev", "y", "x"}); @@ -854,7 +858,7 @@ void solve_radiation(int argc, char** argv) nc_lw_flux_dn .insert(lw_flux_dn_cpu .v(), {0, 0, 0}); nc_lw_flux_net.insert(lw_flux_net_cpu.v(), {0, 0, 0}); - if (switch_single_gpt) + if (lw_single_gpt > 0) { auto nc_lw_gpt_flux_up = output_nc.add_variable("lw_gpt_flux_up" , {"lev", "y", "x"}); auto nc_lw_gpt_flux_dn = output_nc.add_variable("lw_gpt_flux_dn" , {"lev", "y", "x"}); @@ -864,21 +868,21 @@ void solve_radiation(int argc, char** argv) nc_lw_gpt_flux_dn .insert(lw_gpt_flux_dn_cpu.v(), {0, 0, 0}); nc_lw_gpt_flux_net.insert(lw_gpt_flux_net_cpu.v(), {0, 0, 0}); } + } - if (switch_lw_raytracing) - { - auto rt_flux_tod_up = output_nc.add_variable("rt_lw_flux_tod_up" , { "y", "x"}); - auto rt_flux_tod_dn = output_nc.add_variable("rt_lw_flux_tod_dn" , { "y", "x"}); - auto rt_flux_sfc_up = output_nc.add_variable("rt_lw_flux_sfc_up" , { "y", "x"}); - auto rt_flux_sfc_dn = output_nc.add_variable("rt_lw_flux_sfc_dn" , { "y", "x"}); - auto rt_flux_abs = output_nc.add_variable("rt_lw_flux_abs" , {"z", "y", "x"}); - - rt_flux_tod_up.insert(rt_flux_tod_up_cpu.v(), {0, 0}); - rt_flux_tod_dn.insert(rt_flux_tod_dn_cpu.v(), {0, 0}); - rt_flux_sfc_up.insert(rt_flux_sfc_up_cpu.v(), {0, 0}); - rt_flux_sfc_dn.insert(rt_flux_sfc_dn_cpu.v(), {0, 0}); - rt_flux_abs .insert(rt_flux_abs_cpu.v(), {0, 0, 0}); - } + if (switch_lw_raytracing) + { + auto rt_flux_tod_up = output_nc.add_variable("rt_lw_flux_tod_up" , { "y", "x"}); + auto rt_flux_tod_dn = output_nc.add_variable("rt_lw_flux_tod_dn" , { "y", "x"}); + auto rt_flux_sfc_up = output_nc.add_variable("rt_lw_flux_sfc_up" , { "y", "x"}); + auto rt_flux_sfc_dn = output_nc.add_variable("rt_lw_flux_sfc_dn" , { "y", "x"}); + auto rt_flux_abs = output_nc.add_variable("rt_lw_flux_abs" , {"z", "y", "x"}); + + rt_flux_tod_up.insert(rt_flux_tod_up_cpu.v(), {0, 0}); + rt_flux_tod_dn.insert(rt_flux_tod_dn_cpu.v(), {0, 0}); + rt_flux_sfc_up.insert(rt_flux_sfc_up_cpu.v(), {0, 0}); + rt_flux_sfc_dn.insert(rt_flux_sfc_dn_cpu.v(), {0, 0}); + rt_flux_abs .insert(rt_flux_abs_cpu.v(), {0, 0, 0}); } } @@ -925,7 +929,7 @@ void solve_radiation(int argc, char** argv) tsi_scaling({icol}) = Float(1.); } - if (switch_tica) + if (switch_tilted_columns) { for (int icol=1; icol<=n_col; ++icol) tsi_scaling({icol}) = std::cos(tica_sza); @@ -941,7 +945,7 @@ void solve_radiation(int argc, char** argv) Array_gpu sw_aer_ssa; Array_gpu sw_aer_asy; - if (switch_single_gpt) + if (sw_single_gpt > 0) { sw_tot_tau .set_dims({n_col, n_lay}); sw_tot_ssa .set_dims({n_col, n_lay}); @@ -966,26 +970,22 @@ void solve_radiation(int argc, char** argv) Array_gpu rt_flux_abs_dif; - if (switch_fluxes) + if(switch_sw_plane_parallel) { - if(switch_sw_twostream) - { - sw_flux_up .set_dims({n_col, n_lev}); - sw_flux_dn .set_dims({n_col, n_lev}); - sw_flux_dn_dir.set_dims({n_col, n_lev}); - sw_flux_net .set_dims({n_col, n_lev}); - } - - if (switch_sw_raytracing) - { - rt_flux_tod_up .set_dims({n_col_x, n_col_y}); - rt_flux_sfc_dir.set_dims({n_col_x, n_col_y}); - rt_flux_sfc_dif.set_dims({n_col_x, n_col_y}); - rt_flux_sfc_up .set_dims({n_col_x, n_col_y}); - rt_flux_abs_dir.set_dims({n_col_x, n_col_y, n_z}); - rt_flux_abs_dif.set_dims({n_col_x, n_col_y, n_z}); - } + sw_flux_up .set_dims({n_col, n_lev}); + sw_flux_dn .set_dims({n_col, n_lev}); + sw_flux_dn_dir.set_dims({n_col, n_lev}); + sw_flux_net .set_dims({n_col, n_lev}); + } + if (switch_sw_raytracing) + { + rt_flux_tod_up .set_dims({n_col_x, n_col_y}); + rt_flux_sfc_dir.set_dims({n_col_x, n_col_y}); + rt_flux_sfc_dif.set_dims({n_col_x, n_col_y}); + rt_flux_sfc_up .set_dims({n_col_x, n_col_y}); + rt_flux_abs_dir.set_dims({n_col_x, n_col_y, n_z}); + rt_flux_abs_dif.set_dims({n_col_x, n_col_y, n_z}); } Array_gpu sw_gpt_flux_up; @@ -993,7 +993,7 @@ void solve_radiation(int argc, char** argv) Array_gpu sw_gpt_flux_dn_dir; Array_gpu sw_gpt_flux_net; - if (switch_single_gpt) + if (sw_single_gpt > 0) { sw_gpt_flux_up .set_dims({n_col, n_lev}); sw_gpt_flux_dn .set_dims({n_col, n_lev}); @@ -1036,18 +1036,16 @@ void solve_radiation(int argc, char** argv) cudaEventRecord(start, 0); rad_sw.solve_gpu( - switch_fluxes, - switch_sw_twostream, switch_sw_raytracing, - switch_independent_column, + switch_sw_plane_parallel, + switch_sw_independent_column, switch_cloud_optics, switch_cloud_mie, switch_aerosol_optics, - switch_single_gpt, switch_delta_cloud, switch_delta_aerosol, - switch_tica, - single_gpt, + switch_tilted_columns, + sw_single_gpt, sw_photons_per_pixel, grid_cells, grid_d, @@ -1120,7 +1118,7 @@ void solve_radiation(int argc, char** argv) auto nc_sw_band_lims_wvn = output_nc.add_variable("sw_band_lims_wvn", {"band_sw", "pair"}); nc_sw_band_lims_wvn.insert(rad_sw.get_band_lims_wavenumber_gpu().v(), {0, 0}); - if (switch_single_gpt) + if (sw_single_gpt > 0) { auto nc_sw_band_lims_gpt = output_nc.add_variable("sw_band_lims_gpt", {"band_sw", "pair"}); nc_sw_band_lims_gpt.insert(rad_sw.get_band_lims_gpoint_gpu().v(), {0, 0}); @@ -1143,14 +1141,14 @@ void solve_radiation(int argc, char** argv) nc_aer_ssa.insert(sw_aer_ssa_cpu.v(), {0, 0, 0}); nc_aer_asy.insert(sw_aer_asy_cpu.v(), {0, 0, 0}); - nc_tot_tau.add_attribute("long_name","Total optical depth at g-point "+std::to_string(single_gpt)); - nc_tot_ssa.add_attribute("long_name","Total single scattering albedo at g-point "+std::to_string(single_gpt)); - nc_cld_tau.add_attribute("long_name","Cloud optical depth at g-point "+std::to_string(single_gpt)); - nc_cld_ssa.add_attribute("long_name","Cloud single scattering albedo at g-point "+std::to_string(single_gpt)); - nc_cld_asy.add_attribute("long_name","Cloud asymmetry parameter at g-point "+std::to_string(single_gpt)); - nc_aer_tau.add_attribute("long_name","Aerosol optical depth at g-point "+std::to_string(single_gpt)); - nc_aer_ssa.add_attribute("long_name","Aerosol single scattering albedo at g-point "+std::to_string(single_gpt)); - nc_aer_asy.add_attribute("long_name","Aerosol asymmetry parameter at g-point "+std::to_string(single_gpt)); + nc_tot_tau.add_attribute("long_name","Total optical depth at g-point "+std::to_string(sw_single_gpt)); + nc_tot_ssa.add_attribute("long_name","Total single scattering albedo at g-point "+std::to_string(sw_single_gpt)); + nc_cld_tau.add_attribute("long_name","Cloud optical depth at g-point "+std::to_string(sw_single_gpt)); + nc_cld_ssa.add_attribute("long_name","Cloud single scattering albedo at g-point "+std::to_string(sw_single_gpt)); + nc_cld_asy.add_attribute("long_name","Cloud asymmetry parameter at g-point "+std::to_string(sw_single_gpt)); + nc_aer_tau.add_attribute("long_name","Aerosol optical depth at g-point "+std::to_string(sw_single_gpt)); + nc_aer_ssa.add_attribute("long_name","Aerosol single scattering albedo at g-point "+std::to_string(sw_single_gpt)); + nc_aer_asy.add_attribute("long_name","Aerosol asymmetry parameter at g-point "+std::to_string(sw_single_gpt)); nc_tot_tau.add_attribute("units", "-"); nc_tot_ssa.add_attribute("units", "-"); @@ -1163,99 +1161,92 @@ void solve_radiation(int argc, char** argv) } - if (switch_fluxes) + if (switch_sw_plane_parallel) { - if (switch_sw_twostream) - { - auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); - auto nc_sw_flux_dn = output_nc.add_variable("sw_flux_dn" , {"lev", "y", "x"}); - auto nc_sw_flux_dn_dir = output_nc.add_variable("sw_flux_dn_dir", {"lev", "y", "x"}); - auto nc_sw_flux_net = output_nc.add_variable("sw_flux_net" , {"lev", "y", "x"}); + auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); + auto nc_sw_flux_dn = output_nc.add_variable("sw_flux_dn" , {"lev", "y", "x"}); + auto nc_sw_flux_dn_dir = output_nc.add_variable("sw_flux_dn_dir", {"lev", "y", "x"}); + auto nc_sw_flux_net = output_nc.add_variable("sw_flux_net" , {"lev", "y", "x"}); - nc_sw_flux_up .insert(sw_flux_up_cpu .v(), {0, 0, 0}); - nc_sw_flux_dn .insert(sw_flux_dn_cpu .v(), {0, 0, 0}); - nc_sw_flux_dn_dir.insert(sw_flux_dn_dir_cpu.v(), {0, 0, 0}); - nc_sw_flux_net .insert(sw_flux_net_cpu .v(), {0, 0, 0}); + nc_sw_flux_up .insert(sw_flux_up_cpu .v(), {0, 0, 0}); + nc_sw_flux_dn .insert(sw_flux_dn_cpu .v(), {0, 0, 0}); + nc_sw_flux_dn_dir.insert(sw_flux_dn_dir_cpu.v(), {0, 0, 0}); + nc_sw_flux_net .insert(sw_flux_net_cpu .v(), {0, 0, 0}); - nc_sw_flux_up.add_attribute("long_name","Upwelling shortwave fluxes (TwoStream solver)"); - nc_sw_flux_up.add_attribute("units","W m-2"); + nc_sw_flux_up.add_attribute("long_name","Upwelling shortwave fluxes (TwoStream solver)"); + nc_sw_flux_up.add_attribute("units","W m-2"); - nc_sw_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes (TwoStream solver)"); - nc_sw_flux_dn.add_attribute("units","W m-2"); + nc_sw_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes (TwoStream solver)"); + nc_sw_flux_dn.add_attribute("units","W m-2"); - nc_sw_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes (TwoStream solver)"); - nc_sw_flux_dn_dir.add_attribute("units","W m-2"); + nc_sw_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes (TwoStream solver)"); + nc_sw_flux_dn_dir.add_attribute("units","W m-2"); - nc_sw_flux_net.add_attribute("long_name","Net shortwave fluxes (TwoStream solver)"); - nc_sw_flux_net.add_attribute("units","W m-2"); + nc_sw_flux_net.add_attribute("long_name","Net shortwave fluxes (TwoStream solver)"); + nc_sw_flux_net.add_attribute("units","W m-2"); - } - - if (switch_sw_raytracing) + if (sw_single_gpt > 0) { - auto nc_rt_flux_tod_up = output_nc.add_variable("rt_flux_tod_up", {"y","x"}); - auto nc_rt_flux_sfc_dir = output_nc.add_variable("rt_flux_sfc_dir", {"y","x"}); - auto nc_rt_flux_sfc_dif = output_nc.add_variable("rt_flux_sfc_dif", {"y","x"}); - auto nc_rt_flux_sfc_up = output_nc.add_variable("rt_flux_sfc_up", {"y","x"}); - auto nc_rt_flux_abs_dir = output_nc.add_variable("rt_flux_abs_dir", {"z","y","x"}); - auto nc_rt_flux_abs_dif = output_nc.add_variable("rt_flux_abs_dif", {"z","y","x"}); + auto nc_sw_gpt_flux_up = output_nc.add_variable("sw_gpt_flux_up" , {"lev", "y", "x"}); + auto nc_sw_gpt_flux_dn = output_nc.add_variable("sw_gpt_flux_dn" , {"lev", "y", "x"}); + auto nc_sw_gpt_flux_dn_dir = output_nc.add_variable("sw_gpt_flux_dn_dir", {"lev", "y", "x"}); + auto nc_sw_gpt_flux_net = output_nc.add_variable("sw_gpt_flux_net" , {"lev", "y", "x"}); - nc_rt_flux_tod_up .insert(rt_flux_tod_up_cpu .v(), {0,0}); - nc_rt_flux_sfc_dir.insert(rt_flux_sfc_dir_cpu.v(), {0,0}); - nc_rt_flux_sfc_dif.insert(rt_flux_sfc_dif_cpu.v(), {0,0}); - nc_rt_flux_sfc_up .insert(rt_flux_sfc_up_cpu .v(), {0,0}); + nc_sw_gpt_flux_up .insert(sw_gpt_flux_up_cpu .v(), {0, 0, 0}); + nc_sw_gpt_flux_dn .insert(sw_gpt_flux_dn_cpu .v(), {0, 0, 0}); + nc_sw_gpt_flux_dn_dir.insert(sw_gpt_flux_dn_dir_cpu.v(), {0, 0, 0}); + nc_sw_gpt_flux_net .insert(sw_gpt_flux_net_cpu .v(), {0, 0, 0}); - nc_rt_flux_abs_dir.insert(rt_flux_abs_dir_cpu.v(), {0,0,0}); - nc_rt_flux_abs_dif.insert(rt_flux_abs_dif_cpu.v(), {0,0,0}); + nc_sw_gpt_flux_up.add_attribute("long_name","Upwelling shortwave fluxes for g-point "+std::to_string(sw_single_gpt)+" (TwoStream solver)"); + nc_sw_gpt_flux_up.add_attribute("units","W m-2"); - nc_rt_flux_tod_up.add_attribute("long_name","Upwelling shortwave top-of-domain fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_tod_up.add_attribute("units","W m-2"); + nc_sw_gpt_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes for g-point "+std::to_string(sw_single_gpt)+" (TwoStream solver)"); + nc_sw_gpt_flux_dn.add_attribute("units","W m-2"); - nc_rt_flux_sfc_dir.add_attribute("long_name","Downwelling direct shortwave surface fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_sfc_dir.add_attribute("units","W m-2"); + nc_sw_gpt_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes for g-point "+std::to_string(sw_single_gpt)+" (TwoStream solver)"); + nc_sw_gpt_flux_dn_dir.add_attribute("units","W m-2"); - nc_rt_flux_sfc_dif.add_attribute("long_name","Downwelling diffuse shortwave surface fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_sfc_dif.add_attribute("units","W m-2"); + nc_sw_gpt_flux_net.add_attribute("long_name","Net shortwave fluxes for g-point "+std::to_string(sw_single_gpt)+" (TwoStream solver)"); + nc_sw_gpt_flux_net.add_attribute("units","W m-2"); - nc_rt_flux_sfc_up.add_attribute("long_name","Upwelling shortwave surface fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_sfc_up.add_attribute("units","W m-2"); + } - nc_rt_flux_abs_dir.add_attribute("long_name","Absorbed direct shortwave fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_abs_dir.add_attribute("units","W m-3"); + } - nc_rt_flux_abs_dif.add_attribute("long_name","Absorbed diffuse shortwave fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_abs_dif.add_attribute("units","W m-3"); + if (switch_sw_raytracing) + { + auto nc_rt_flux_tod_up = output_nc.add_variable("rt_flux_tod_up", {"y","x"}); + auto nc_rt_flux_sfc_dir = output_nc.add_variable("rt_flux_sfc_dir", {"y","x"}); + auto nc_rt_flux_sfc_dif = output_nc.add_variable("rt_flux_sfc_dif", {"y","x"}); + auto nc_rt_flux_sfc_up = output_nc.add_variable("rt_flux_sfc_up", {"y","x"}); + auto nc_rt_flux_abs_dir = output_nc.add_variable("rt_flux_abs_dir", {"z","y","x"}); + auto nc_rt_flux_abs_dif = output_nc.add_variable("rt_flux_abs_dif", {"z","y","x"}); - } + nc_rt_flux_tod_up .insert(rt_flux_tod_up_cpu .v(), {0,0}); + nc_rt_flux_sfc_dir.insert(rt_flux_sfc_dir_cpu.v(), {0,0}); + nc_rt_flux_sfc_dif.insert(rt_flux_sfc_dif_cpu.v(), {0,0}); + nc_rt_flux_sfc_up .insert(rt_flux_sfc_up_cpu .v(), {0,0}); + nc_rt_flux_abs_dir.insert(rt_flux_abs_dir_cpu.v(), {0,0,0}); + nc_rt_flux_abs_dif.insert(rt_flux_abs_dif_cpu.v(), {0,0,0}); - if (switch_single_gpt) - { - if (switch_sw_twostream) - { - auto nc_sw_gpt_flux_up = output_nc.add_variable("sw_gpt_flux_up" , {"lev", "y", "x"}); - auto nc_sw_gpt_flux_dn = output_nc.add_variable("sw_gpt_flux_dn" , {"lev", "y", "x"}); - auto nc_sw_gpt_flux_dn_dir = output_nc.add_variable("sw_gpt_flux_dn_dir", {"lev", "y", "x"}); - auto nc_sw_gpt_flux_net = output_nc.add_variable("sw_gpt_flux_net" , {"lev", "y", "x"}); + nc_rt_flux_tod_up.add_attribute("long_name","Upwelling shortwave top-of-domain fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_tod_up.add_attribute("units","W m-2"); - nc_sw_gpt_flux_up .insert(sw_gpt_flux_up_cpu .v(), {0, 0, 0}); - nc_sw_gpt_flux_dn .insert(sw_gpt_flux_dn_cpu .v(), {0, 0, 0}); - nc_sw_gpt_flux_dn_dir.insert(sw_gpt_flux_dn_dir_cpu.v(), {0, 0, 0}); - nc_sw_gpt_flux_net .insert(sw_gpt_flux_net_cpu .v(), {0, 0, 0}); + nc_rt_flux_sfc_dir.add_attribute("long_name","Downwelling direct shortwave surface fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_sfc_dir.add_attribute("units","W m-2"); - nc_sw_gpt_flux_up.add_attribute("long_name","Upwelling shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_up.add_attribute("units","W m-2"); + nc_rt_flux_sfc_dif.add_attribute("long_name","Downwelling diffuse shortwave surface fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_sfc_dif.add_attribute("units","W m-2"); - nc_sw_gpt_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_dn.add_attribute("units","W m-2"); + nc_rt_flux_sfc_up.add_attribute("long_name","Upwelling shortwave surface fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_sfc_up.add_attribute("units","W m-2"); - nc_sw_gpt_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_dn_dir.add_attribute("units","W m-2"); + nc_rt_flux_abs_dir.add_attribute("long_name","Absorbed direct shortwave fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_abs_dir.add_attribute("units","W m-3"); - nc_sw_gpt_flux_net.add_attribute("long_name","Net shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_net.add_attribute("units","W m-2"); - } - } + nc_rt_flux_abs_dif.add_attribute("long_name","Absorbed diffuse shortwave fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_abs_dif.add_attribute("units","W m-3"); } } @@ -1374,7 +1365,7 @@ void solve_radiation(int argc, char** argv) grid_cells, grid_d, kn_grid, - photons_per_pixel, + bw_photons_per_pixel, gas_concs_gpu, p_lay_gpu, p_lev_gpu, t_lay_gpu, t_lev_gpu, @@ -1448,7 +1439,7 @@ void solve_radiation(int argc, char** argv) grid_cells, grid_d, kn_grid, - photons_per_pixel, + bw_photons_per_pixel, gas_concs_gpu, p_lay_gpu, p_lev_gpu, t_lay_gpu, t_lev_gpu, From 6fcb8d5cfe083df53e5ef716246ff4c7f4a2d821 Mon Sep 17 00:00:00 2001 From: Menno Date: Wed, 4 Mar 2026 16:59:39 +0100 Subject: [PATCH 03/11] change the way single_gpt works, now only solve and store the requested g-point --- include_test/radiation_solver_rt.h | 3 -- src_test/radiation_solver_rt.cu | 30 ++++-------- src_test/test_rte_rrtmgp_rt.cu | 73 ------------------------------ 3 files changed, 8 insertions(+), 98 deletions(-) diff --git a/include_test/radiation_solver_rt.h b/include_test/radiation_solver_rt.h index 2271cda..d26029c 100644 --- a/include_test/radiation_solver_rt.h +++ b/include_test/radiation_solver_rt.h @@ -74,7 +74,6 @@ class Radiation_solver_longwave Array_gpu& aer_tau_out,Array_gpu& aer_ssa_out, Array_gpu& aer_asy_out, Array_gpu& lay_source, Array_gpu& lev_source, Array_gpu& sfc_source, Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net, - Array_gpu& lw_gpt_flux_up, Array_gpu& lw_gpt_flux_dn, Array_gpu& lw_gpt_flux_net, Array_gpu& rt_flux_tod_up, Array_gpu& rt_flux_tod_dn, Array_gpu& rt_flux_sfc_up, Array_gpu& rt_flux_sfc_dn, Array_gpu& rt_flux_abs); @@ -151,8 +150,6 @@ class Radiation_solver_shortwave Array_gpu& aer_tau_out, Array_gpu& aer_ssa_out, Array_gpu& aer_asy_out, Array_gpu& sw_flux_up, Array_gpu& sw_flux_dn, Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net, - Array_gpu& sw_gpt_flux_up, Array_gpu& sw_gpt_flux_dn, - Array_gpu& sw_gpt_flux_dn_dir, Array_gpu& sw_gpt_flux_net, Array_gpu& rt_flux_tod_up, Array_gpu& rt_flux_sfc_dir, Array_gpu& rt_flux_sfc_dif, diff --git a/src_test/radiation_solver_rt.cu b/src_test/radiation_solver_rt.cu index d4540af..fb49e76 100644 --- a/src_test/radiation_solver_rt.cu +++ b/src_test/radiation_solver_rt.cu @@ -507,7 +507,6 @@ void Radiation_solver_longwave::solve_gpu( Array_gpu& aer_tau_out, Array_gpu& aer_ssa_out, Array_gpu& aer_asy_out, Array_gpu& lay_source, Array_gpu& lev_source, Array_gpu& sfc_source, Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net, - Array_gpu& lw_gpt_flux_up, Array_gpu& lw_gpt_flux_dn, Array_gpu& lw_gpt_flux_net, Array_gpu& rt_flux_tod_up, Array_gpu& rt_flux_tod_dn, Array_gpu& rt_flux_sfc_up, Array_gpu& rt_flux_sfc_dn, Array_gpu& rt_flux_abs) { @@ -548,6 +547,8 @@ void Radiation_solver_longwave::solve_gpu( for (int igpt=1; igpt<=n_gpt; ++igpt) { + if ((single_gpt > 0) && (igpt != single_gpt)) continue; + int band = 0; for (int ibnd=1; ibnd<=n_bnd; ++ibnd) { @@ -688,8 +689,8 @@ void Radiation_solver_longwave::solve_gpu( Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_g().ptr()); } - // Store the optical properties, if desired. - if (igpt == single_gpt) + // Store the optical properties, only in single_gpt mode + if (single_gpt > 0) { lay_source = (*sources).get_lay_source(); lev_source = (*sources).get_lev_source(); @@ -738,13 +739,6 @@ void Radiation_solver_longwave::solve_gpu( n_col, n_lev, lw_flux_up.ptr(), lw_flux_dn.ptr(), lw_flux_net.ptr(), (*fluxes).get_flux_up().ptr(), (*fluxes).get_flux_dn().ptr(), (*fluxes).get_flux_net().ptr()); } - - if (igpt == single_gpt) - { - lw_gpt_flux_up = (*fluxes).get_flux_up(); - lw_gpt_flux_dn = (*fluxes).get_flux_dn(); - lw_gpt_flux_net = (*fluxes).get_flux_net(); - } } if (switch_raytracing) @@ -877,8 +871,6 @@ void Radiation_solver_shortwave::solve_gpu( Array_gpu& aer_tau_out, Array_gpu& aer_ssa_out, Array_gpu& aer_asy_out, Array_gpu& sw_flux_up, Array_gpu& sw_flux_dn, Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net, - Array_gpu& sw_gpt_flux_up, Array_gpu& sw_gpt_flux_dn, - Array_gpu& sw_gpt_flux_dn_dir, Array_gpu& sw_gpt_flux_net, Array_gpu& rt_flux_tod_up, Array_gpu& rt_flux_sfc_dir, Array_gpu& rt_flux_sfc_dif, @@ -937,6 +929,8 @@ void Radiation_solver_shortwave::solve_gpu( for (int igpt=1; igpt<=n_gpt; ++igpt) { + if ((single_gpt > 0) && (igpt != single_gpt)) continue; + int band = 0; for (int ibnd=1; ibnd<=n_bnd; ++ibnd) { @@ -1060,8 +1054,8 @@ void Radiation_solver_shortwave::solve_gpu( Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_g().ptr()); } - // Store the optical properties, if desired - if (igpt == single_gpt) + // Store the optical properties, only in single_gpt mode + if (single_gpt > 0) { tot_tau_out = optical_props->get_tau(); tot_ssa_out = optical_props->get_ssa(); @@ -1078,7 +1072,6 @@ void Radiation_solver_shortwave::solve_gpu( if (switch_plane_parallel) { - rte_sw.rte_sw( optical_props, top_at_1, @@ -1154,13 +1147,6 @@ void Radiation_solver_shortwave::solve_gpu( (*fluxes).get_flux_abs_dir().ptr(), (*fluxes).get_flux_abs_dif().ptr()); } - if (igpt == single_gpt) - { - sw_gpt_flux_up = (*fluxes).get_flux_up(); - sw_gpt_flux_dn = (*fluxes).get_flux_dn(); - sw_gpt_flux_dn_dir = (*fluxes).get_flux_dn_dir(); - sw_gpt_flux_net = (*fluxes).get_flux_net(); - } previous_band = band; } } diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index 00ba0f5..bb6124d 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -662,17 +662,6 @@ void solve_radiation(int argc, char** argv) lw_flux_net.set_dims({n_col, n_lev}); } - Array_gpu lw_gpt_flux_up; - Array_gpu lw_gpt_flux_dn; - Array_gpu lw_gpt_flux_net; - - if (lw_single_gpt > 0) - { - lw_gpt_flux_up .set_dims({n_col, n_lev}); - lw_gpt_flux_dn .set_dims({n_col, n_lev}); - lw_gpt_flux_net.set_dims({n_col, n_lev}); - } - Array_gpu rt_flux_tod_up; Array_gpu rt_flux_tod_dn; Array_gpu rt_flux_sfc_up; @@ -750,7 +739,6 @@ void solve_radiation(int argc, char** argv) lw_tau_aer, lw_ssa_aer, lw_asy_aer, lay_source, lev_source, sfc_source, lw_flux_up, lw_flux_dn, lw_flux_net, - lw_gpt_flux_up, lw_gpt_flux_dn, lw_gpt_flux_net, rt_flux_tod_up, rt_flux_tod_dn, rt_flux_sfc_up, rt_flux_sfc_dn, rt_flux_abs); @@ -792,9 +780,6 @@ void solve_radiation(int argc, char** argv) Array lw_flux_up_cpu(lw_flux_up); Array lw_flux_dn_cpu(lw_flux_dn); Array lw_flux_net_cpu(lw_flux_net); - Array lw_gpt_flux_up_cpu(lw_gpt_flux_up); - Array lw_gpt_flux_dn_cpu(lw_gpt_flux_dn); - Array lw_gpt_flux_net_cpu(lw_gpt_flux_net); Array rt_flux_tod_up_cpu(rt_flux_tod_up); Array rt_flux_tod_dn_cpu(rt_flux_tod_dn); @@ -857,17 +842,6 @@ void solve_radiation(int argc, char** argv) nc_lw_flux_up .insert(lw_flux_up_cpu .v(), {0, 0, 0}); nc_lw_flux_dn .insert(lw_flux_dn_cpu .v(), {0, 0, 0}); nc_lw_flux_net.insert(lw_flux_net_cpu.v(), {0, 0, 0}); - - if (lw_single_gpt > 0) - { - auto nc_lw_gpt_flux_up = output_nc.add_variable("lw_gpt_flux_up" , {"lev", "y", "x"}); - auto nc_lw_gpt_flux_dn = output_nc.add_variable("lw_gpt_flux_dn" , {"lev", "y", "x"}); - auto nc_lw_gpt_flux_net = output_nc.add_variable("lw_gpt_flux_net", {"lev", "y", "x"}); - - nc_lw_gpt_flux_up .insert(lw_gpt_flux_up_cpu.v(), {0, 0, 0}); - nc_lw_gpt_flux_dn .insert(lw_gpt_flux_dn_cpu.v(), {0, 0, 0}); - nc_lw_gpt_flux_net.insert(lw_gpt_flux_net_cpu.v(), {0, 0, 0}); - } } if (switch_lw_raytracing) @@ -988,20 +962,6 @@ void solve_radiation(int argc, char** argv) rt_flux_abs_dif.set_dims({n_col_x, n_col_y, n_z}); } - Array_gpu sw_gpt_flux_up; - Array_gpu sw_gpt_flux_dn; - Array_gpu sw_gpt_flux_dn_dir; - Array_gpu sw_gpt_flux_net; - - if (sw_single_gpt > 0) - { - sw_gpt_flux_up .set_dims({n_col, n_lev}); - sw_gpt_flux_dn .set_dims({n_col, n_lev}); - sw_gpt_flux_dn_dir.set_dims({n_col, n_lev}); - sw_gpt_flux_net .set_dims({n_col, n_lev}); - } - - // Solve the radiation. Status::print_message("Solving the shortwave radiation."); @@ -1065,8 +1025,6 @@ void solve_radiation(int argc, char** argv) sw_aer_tau, sw_aer_ssa, sw_aer_asy, sw_flux_up, sw_flux_dn, sw_flux_dn_dir, sw_flux_net, - sw_gpt_flux_up, sw_gpt_flux_dn, - sw_gpt_flux_dn_dir, sw_gpt_flux_net, rt_flux_tod_up, rt_flux_sfc_dir, rt_flux_sfc_dif, @@ -1103,10 +1061,6 @@ void solve_radiation(int argc, char** argv) Array sw_flux_dn_cpu(sw_flux_dn); Array sw_flux_dn_dir_cpu(sw_flux_dn_dir); Array sw_flux_net_cpu(sw_flux_net); - Array sw_gpt_flux_up_cpu(sw_gpt_flux_up); - Array sw_gpt_flux_dn_cpu(sw_gpt_flux_dn); - Array sw_gpt_flux_dn_dir_cpu(sw_gpt_flux_dn_dir); - Array sw_gpt_flux_net_cpu(sw_gpt_flux_net); Array rt_flux_tod_up_cpu(rt_flux_tod_up); Array rt_flux_sfc_dir_cpu(rt_flux_sfc_dir); @@ -1184,33 +1138,6 @@ void solve_radiation(int argc, char** argv) nc_sw_flux_net.add_attribute("long_name","Net shortwave fluxes (TwoStream solver)"); nc_sw_flux_net.add_attribute("units","W m-2"); - - if (sw_single_gpt > 0) - { - auto nc_sw_gpt_flux_up = output_nc.add_variable("sw_gpt_flux_up" , {"lev", "y", "x"}); - auto nc_sw_gpt_flux_dn = output_nc.add_variable("sw_gpt_flux_dn" , {"lev", "y", "x"}); - auto nc_sw_gpt_flux_dn_dir = output_nc.add_variable("sw_gpt_flux_dn_dir", {"lev", "y", "x"}); - auto nc_sw_gpt_flux_net = output_nc.add_variable("sw_gpt_flux_net" , {"lev", "y", "x"}); - - nc_sw_gpt_flux_up .insert(sw_gpt_flux_up_cpu .v(), {0, 0, 0}); - nc_sw_gpt_flux_dn .insert(sw_gpt_flux_dn_cpu .v(), {0, 0, 0}); - nc_sw_gpt_flux_dn_dir.insert(sw_gpt_flux_dn_dir_cpu.v(), {0, 0, 0}); - nc_sw_gpt_flux_net .insert(sw_gpt_flux_net_cpu .v(), {0, 0, 0}); - - nc_sw_gpt_flux_up.add_attribute("long_name","Upwelling shortwave fluxes for g-point "+std::to_string(sw_single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_up.add_attribute("units","W m-2"); - - nc_sw_gpt_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes for g-point "+std::to_string(sw_single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_dn.add_attribute("units","W m-2"); - - nc_sw_gpt_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes for g-point "+std::to_string(sw_single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_dn_dir.add_attribute("units","W m-2"); - - nc_sw_gpt_flux_net.add_attribute("long_name","Net shortwave fluxes for g-point "+std::to_string(sw_single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_net.add_attribute("units","W m-2"); - - } - } if (switch_sw_raytracing) From 7f1bb0d17ed7cf3b6de5913a8611a6a87793d112 Mon Sep 17 00:00:00 2001 From: Menno Date: Wed, 4 Mar 2026 17:11:43 +0100 Subject: [PATCH 04/11] add rt_ to independent_column flag --- src_test/test_rte_rrtmgp_rt.cu | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index bb6124d..d9a811b 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -179,14 +179,14 @@ void solve_radiation(int argc, char** argv) const bool switch_shortwave = get_ini_value(settings, "shortwave", "shortwave", true); const bool switch_sw_raytracing = get_ini_value(settings, "shortwave", "raytracing", true); // compute and output ray tracer fluxes const bool switch_sw_plane_parallel = get_ini_value(settings, "shortwave", "plane_parallel", true); // compute and output plane parallel 1D fluxes (two-stream) - bool switch_sw_independent_column = get_ini_value(settings, "shortwave", "independent_column", false); // solve ray tracer in independent column mode + bool switch_sw_independent_column = get_ini_value(settings, "shortwave", "rt_independent_column", false); // solve ray tracer in independent column mode const int sw_single_gpt = get_ini_value(settings, "shortwave", "single_gpt", -1); // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) const bool switch_longwave = get_ini_value(settings, "longwave", "longwave", true); const bool switch_lw_raytracing = get_ini_value(settings, "longwave", "raytracing", true); // compute and output ray tracer fluxes const bool switch_lw_plane_parallel = get_ini_value(settings, "longwave", "plane_parallel", true); // compute and output plane parallel 1D fluxes (two-stream or no-scattering solution) const bool switch_lw_scattering = get_ini_value(settings, "longwave", "scattering", false); // enable scattering in longwave solver (only possible in ray tracer executable) - bool switch_lw_independent_column = get_ini_value(settings, "longwave", "independent_column", false); // solve ray tracer in independent column mode + bool switch_lw_independent_column = get_ini_value(settings, "longwave", "rt_independent_column", false); // solve ray tracer in independent column mode const int lw_single_gpt = get_ini_value(settings, "longwave", "single_gpt", -1); // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) const Float min_mfp_grid_ratio = get_ini_value(settings, "longwave", "min_mfp_grid_ratio", Float(1.0)); // minimum ratio between the lowest gasous mean fee path and the horizontal grid spacing at which ray tracer is still used. Set to 0. to use ray tracer for all g-points From 8a549ef50f71094c3bd042fa389cc719923b3a76 Mon Sep 17 00:00:00 2001 From: Menno Date: Wed, 4 Mar 2026 17:13:02 +0100 Subject: [PATCH 05/11] change default single_gpt to 0, gpt counters start at 1 --- src_test/test_rte_rrtmgp_rt.cu | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index d9a811b..6760dca 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -180,14 +180,14 @@ void solve_radiation(int argc, char** argv) const bool switch_sw_raytracing = get_ini_value(settings, "shortwave", "raytracing", true); // compute and output ray tracer fluxes const bool switch_sw_plane_parallel = get_ini_value(settings, "shortwave", "plane_parallel", true); // compute and output plane parallel 1D fluxes (two-stream) bool switch_sw_independent_column = get_ini_value(settings, "shortwave", "rt_independent_column", false); // solve ray tracer in independent column mode - const int sw_single_gpt = get_ini_value(settings, "shortwave", "single_gpt", -1); // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) + const int sw_single_gpt = get_ini_value(settings, "shortwave", "single_gpt", 0); // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) const bool switch_longwave = get_ini_value(settings, "longwave", "longwave", true); const bool switch_lw_raytracing = get_ini_value(settings, "longwave", "raytracing", true); // compute and output ray tracer fluxes const bool switch_lw_plane_parallel = get_ini_value(settings, "longwave", "plane_parallel", true); // compute and output plane parallel 1D fluxes (two-stream or no-scattering solution) const bool switch_lw_scattering = get_ini_value(settings, "longwave", "scattering", false); // enable scattering in longwave solver (only possible in ray tracer executable) bool switch_lw_independent_column = get_ini_value(settings, "longwave", "rt_independent_column", false); // solve ray tracer in independent column mode - const int lw_single_gpt = get_ini_value(settings, "longwave", "single_gpt", -1); // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) + const int lw_single_gpt = get_ini_value(settings, "longwave", "single_gpt", 0); // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) const Float min_mfp_grid_ratio = get_ini_value(settings, "longwave", "min_mfp_grid_ratio", Float(1.0)); // minimum ratio between the lowest gasous mean fee path and the horizontal grid spacing at which ray tracer is still used. Set to 0. to use ray tracer for all g-points const bool switch_bw_raytracing = get_ini_value(settings, "backward", "bw_raytracing", true); From 352e1ca1a462d6bba98d9531c05f4fd4e20a712c Mon Sep 17 00:00:00 2001 From: Menno Date: Thu, 5 Mar 2026 11:43:26 +0100 Subject: [PATCH 06/11] fix some swapped switches --- src_test/radiation_solver_rt.cu | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src_test/radiation_solver_rt.cu b/src_test/radiation_solver_rt.cu index fb49e76..21b8b5c 100644 --- a/src_test/radiation_solver_rt.cu +++ b/src_test/radiation_solver_rt.cu @@ -479,8 +479,8 @@ Radiation_solver_longwave::Radiation_solver_longwave( void Radiation_solver_longwave::solve_gpu( - const bool switch_plane_parallel, const bool switch_raytracing, + const bool switch_plane_parallel, const bool switch_cloud_optics, const bool switch_aerosol_optics, const bool switch_delta_cloud, @@ -536,8 +536,6 @@ void Radiation_solver_longwave::solve_gpu( Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, n_col, lw_flux_up.ptr()); Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, n_col, lw_flux_dn.ptr()); Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, n_col, lw_flux_net.ptr()); - - if (switch_raytracing) { } } const Array& band_limits_gpt(this->kdist_gpu->get_band_lims_gpoint()); @@ -604,7 +602,6 @@ void Radiation_solver_longwave::solve_gpu( gas_optics_subset(col_s, n_col_residual); } - // Find maximum gasous optical depth to and compute lowest mean free path on the clearsky atmosphere // Allocate temporary storage @@ -841,8 +838,8 @@ void Radiation_solver_shortwave::load_mie_tables( } void Radiation_solver_shortwave::solve_gpu( - const bool switch_plane_parallel, const bool switch_raytracing, + const bool switch_plane_parallel, const bool switch_independent_column, const bool switch_cloud_optics, const bool switch_cloud_mie, From 3c81f349dddb15fb2acb54885ca2d4c7eb471fab Mon Sep 17 00:00:00 2001 From: Menno Date: Thu, 5 Mar 2026 12:13:28 +0100 Subject: [PATCH 07/11] add explanatory comments above get_ini_value lines and add comments to camera reads --- src_test/radiation_solver_rt.cu | 18 +++++----- src_test/test_rte_rrtmgp_rt.cu | 64 +++++++++++++++++++++++---------- 2 files changed, 54 insertions(+), 28 deletions(-) diff --git a/src_test/radiation_solver_rt.cu b/src_test/radiation_solver_rt.cu index 21b8b5c..fa9e9ec 100644 --- a/src_test/radiation_solver_rt.cu +++ b/src_test/radiation_solver_rt.cu @@ -1133,16 +1133,16 @@ void Radiation_solver_shortwave::solve_gpu( (*fluxes).get_flux_up().ptr(), (*fluxes).get_flux_dn().ptr(), (*fluxes).get_flux_dn_dir().ptr(), (*fluxes).get_flux_net().ptr()); } - if (switch_raytracing) - { - Gpt_combine_kernels_cuda_rt::add_from_gpoint( - grid_cells.x, grid_cells.y, rt_flux_tod_up.ptr(), rt_flux_sfc_dir.ptr(), rt_flux_sfc_dif.ptr(), rt_flux_sfc_up.ptr(), - (*fluxes).get_flux_tod_up().ptr(), (*fluxes).get_flux_sfc_dir().ptr(), (*fluxes).get_flux_sfc_dif().ptr(), (*fluxes).get_flux_sfc_up().ptr()); + if (switch_raytracing) + { + Gpt_combine_kernels_cuda_rt::add_from_gpoint( + grid_cells.x, grid_cells.y, rt_flux_tod_up.ptr(), rt_flux_sfc_dir.ptr(), rt_flux_sfc_dif.ptr(), rt_flux_sfc_up.ptr(), + (*fluxes).get_flux_tod_up().ptr(), (*fluxes).get_flux_sfc_dir().ptr(), (*fluxes).get_flux_sfc_dif().ptr(), (*fluxes).get_flux_sfc_up().ptr()); - Gpt_combine_kernels_cuda_rt::add_from_gpoint( - n_col, grid_cells.z, rt_flux_abs_dir.ptr(), rt_flux_abs_dif.ptr(), - (*fluxes).get_flux_abs_dir().ptr(), (*fluxes).get_flux_abs_dif().ptr()); - } + Gpt_combine_kernels_cuda_rt::add_from_gpoint( + n_col, grid_cells.z, rt_flux_abs_dir.ptr(), rt_flux_abs_dif.ptr(), + (*fluxes).get_flux_abs_dir().ptr(), (*fluxes).get_flux_abs_dif().ptr()); + } previous_band = band; } diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index 6760dca..f995e82 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -171,44 +171,70 @@ void solve_radiation(int argc, char** argv) const bool switch_ice_cloud_optics = get_ini_value(settings, "clouds", "ice_cloud_optics", false); const bool switch_delta_cloud = get_ini_value(settings, "clouds", "delta_cloud", false); const bool switch_tilted_columns = get_ini_value(settings, "clouds", "tilted_columns", false); - const bool switch_cloud_mie = get_ini_value(settings, "clouds", "cloud_mie", false); // use mie scattering. Curerntly only for shortwave and only for liquid cloud droplets (requires ice_cloud_optics = false) + + // use mie scattering. Currently only for shortwave and only for liquid cloud droplets (requires ice_cloud_optics = false) + const bool switch_cloud_mie = get_ini_value(settings, "clouds", "cloud_mie", false); const bool switch_aerosol_optics = get_ini_value(settings, "aerosols", "aerosol_optics", false); const bool switch_delta_aerosol = get_ini_value(settings, "aerosols", "delta_aerosol", false); const bool switch_shortwave = get_ini_value(settings, "shortwave", "shortwave", true); - const bool switch_sw_raytracing = get_ini_value(settings, "shortwave", "raytracing", true); // compute and output ray tracer fluxes - const bool switch_sw_plane_parallel = get_ini_value(settings, "shortwave", "plane_parallel", true); // compute and output plane parallel 1D fluxes (two-stream) - bool switch_sw_independent_column = get_ini_value(settings, "shortwave", "rt_independent_column", false); // solve ray tracer in independent column mode - const int sw_single_gpt = get_ini_value(settings, "shortwave", "single_gpt", 0); // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) + // compute and output ray tracer fluxes + const bool switch_sw_raytracing = get_ini_value(settings, "shortwave", "raytracing", true); + // compute and output plane parallel 1D fluxes (two-stream) + const bool switch_sw_plane_parallel = get_ini_value(settings, "shortwave", "plane_parallel", true); + // solve ray tracer in independent column mode + bool switch_sw_independent_column = get_ini_value(settings, "shortwave", "rt_independent_column", false); + // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) + const int sw_single_gpt = get_ini_value(settings, "shortwave", "single_gpt", 0); const bool switch_longwave = get_ini_value(settings, "longwave", "longwave", true); - const bool switch_lw_raytracing = get_ini_value(settings, "longwave", "raytracing", true); // compute and output ray tracer fluxes - const bool switch_lw_plane_parallel = get_ini_value(settings, "longwave", "plane_parallel", true); // compute and output plane parallel 1D fluxes (two-stream or no-scattering solution) - const bool switch_lw_scattering = get_ini_value(settings, "longwave", "scattering", false); // enable scattering in longwave solver (only possible in ray tracer executable) - bool switch_lw_independent_column = get_ini_value(settings, "longwave", "rt_independent_column", false); // solve ray tracer in independent column mode - const int lw_single_gpt = get_ini_value(settings, "longwave", "single_gpt", 0); // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) - const Float min_mfp_grid_ratio = get_ini_value(settings, "longwave", "min_mfp_grid_ratio", Float(1.0)); // minimum ratio between the lowest gasous mean fee path and the horizontal grid spacing at which ray tracer is still used. Set to 0. to use ray tracer for all g-points + // compute and output ray tracer fluxes + const bool switch_lw_raytracing = get_ini_value(settings, "longwave", "raytracing", true); + // compute and output plane parallel 1D fluxes (two-stream or no-scattering solution) + const bool switch_lw_plane_parallel = get_ini_value(settings, "longwave", "plane_parallel", true); + // enable scattering in longwave solver (only possible in ray tracer executable) + const bool switch_lw_scattering = get_ini_value(settings, "longwave", "scattering", false); + // solve ray tracer in independent column mode + bool switch_lw_independent_column = get_ini_value(settings, "longwave", "rt_independent_column", false); + // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) + const int lw_single_gpt = get_ini_value(settings, "longwave", "single_gpt", 0); + // minimum ratio between the lowest gaseous mean free path and the horizontal grid spacing at which ray tracer is still used. Set to 0. to use ray tracer for all g-points + const Float min_mfp_grid_ratio = get_ini_value(settings, "longwave", "min_mfp_grid_ratio", Float(1.0)); const bool switch_bw_raytracing = get_ini_value(settings, "backward", "bw_raytracing", true); - const bool switch_lu_albedo = get_ini_value(settings, "backward", "lu_albedo", false); // read surface type from land_use_map variabele and compute spectral albedo and reflection type (lambertian/specular) accordingly - const bool switch_image = get_ini_value(settings, "backward", "image", true); // solve visible bands and convert to XYZ tristimulus values - const bool switch_broadband = get_ini_value(settings, "backward", "broadband", false); // solve broadband radiances - const bool switch_cloud_cam = get_ini_value(settings, "backward", "cloud_cam", false); // output additional cloud statistics for each camera pixel - - const Float input_sza = get_ini_value(settings, "solar_angles", "sza", -1.0); // if >0, overwrite zenith angle from input netcdf file - const Float input_azi = get_ini_value(settings, "solar_angles", "azi", -1.0); // if >0, overwrite azimuth angle from input netcdf file + // read surface type from land_use_map variable and compute spectral albedo and reflection type (lambertian/specular) accordingly + const bool switch_lu_albedo = get_ini_value(settings, "backward", "lu_albedo", false); + // solve visible bands and convert to XYZ tristimulus values + const bool switch_image = get_ini_value(settings, "backward", "image", true); + // solve broadband radiances + const bool switch_broadband = get_ini_value(settings, "backward", "broadband", false); + // output additional cloud statistics for each camera pixel + const bool switch_cloud_cam = get_ini_value(settings, "backward", "cloud_cam", false); + + // if >0, overwrite zenith angle from input netcdf file + const Float input_sza = get_ini_value(settings, "solar_angles", "sza", -1.0); + // if >0, overwrite azimuth angle from input netcdf file + const Float input_azi = get_ini_value(settings, "solar_angles", "azi", -1.0); Camera camera; if (switch_bw_raytracing) { - camera.fov = get_ini_value(settings, "camera", "cam_field-of-view", 80.); + camera.fov = get_ini_value(settings, "camera", "cam_field_of_view", 80.); + + // camera type: (0) fish eye camera, (1) rectangular camera, (2) top-of-atmosphere upwelling radiances camera.cam_type = get_ini_value(settings, "camera", "cam_type", 0); + + // x,y,z positions of virtual camera camera.position = {get_ini_value(settings, "camera", "cam_px", 0.), get_ini_value(settings, "camera", "cam_py", 0.), get_ini_value(settings, "camera", "cam_pz", 100.)}; + + // width, height (pixels) of virtual camera or number of zenith and azimuth angles of fish camera camera.nx = get_ini_value(settings, "camera", "cam_nx", 0); camera.ny = get_ini_value(settings, "camera", "cam_ny", 0); + + // yaw, pitch and roll angles (degrees) of the virtual camera camera.setup_rotation_matrix(get_ini_value(settings, "camera", "cam_yaw", 0.), get_ini_value(settings, "camera", "cam_pitch", 0.), get_ini_value(settings, "camera", "cam_roll", 0.)); From ddac2d9996b4e33a05ce35daf7d9ca353664680e Mon Sep 17 00:00:00 2001 From: Menno Date: Thu, 5 Mar 2026 12:20:00 +0100 Subject: [PATCH 08/11] printing number of backward photons to screen --- src_test/test_rte_rrtmgp_rt.cu | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index f995e82..8d597dd 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -246,13 +246,13 @@ void solve_radiation(int argc, char** argv) if (switch_shortwave && !(switch_sw_plane_parallel || switch_sw_raytracing)) { - std::string error = "With shortwave enabled, need to run either shortwave plane parallel solver or ray tracer "; + std::string error = "With shortwave enabled, need to run shortwave plane parallel solver and/or ray tracer "; throw std::runtime_error(error); } if (switch_longwave && !(switch_lw_plane_parallel || switch_lw_raytracing)) { - std::string error = "With longwave enabled, need to run either longwave plane parallel solver or ray tracer "; + std::string error = "With longwave enabled, need to run longwave plane parallel solver and/or ray tracer "; throw std::runtime_error(error); } @@ -306,6 +306,8 @@ void solve_radiation(int argc, char** argv) if (switch_lw_raytracing) Status::print_message("Longwave: using 2**"+std::to_string(lw_photon_power) + " ("+std::to_string(lw_photon_count) + ") rays per g-point"); + if (switch_bw_raytracing) + Status::print_message("Backward: using "+ std::to_string(bw_photons_per_pixel) + " rays per pixel per g-point"); ////// READ THE ATMOSPHERIC DATA ////// Status::print_message("Reading atmospheric input data from NetCDF."); From 037523d46833add9c43a1795ce90d03f9b903c4e Mon Sep 17 00:00:00 2001 From: Menno Date: Thu, 5 Mar 2026 12:38:11 +0100 Subject: [PATCH 09/11] update and simplify test case input files. FUll range of options will be provided in example case --- allsky/allsky.ini | 49 ++++++++++++++++++++++------------------------- rcemip/rcemip.ini | 49 ++++++++++++++++++++++------------------------- rfmip/rfmip.ini | 49 ++++++++++++++++++++++------------------------- 3 files changed, 69 insertions(+), 78 deletions(-) diff --git a/allsky/allsky.ini b/allsky/allsky.ini index 8e1fea5..8990a87 100644 --- a/allsky/allsky.ini +++ b/allsky/allsky.ini @@ -1,32 +1,29 @@ -[switches] -shortwave = true -longwave = true -fluxes = true +[clouds] liquid_cloud_optics = true ice_cloud_optics = true +delta_cloud = false +tilted_columns = false +cloud_mie = false + +[aerosols] aerosol_optics = false -delta_cloud = true delta_aerosol = false -# gpu_only (test_rte_rrtmgp_gpu) -timings = false +[shortwave] +shortwave = true + +# compute ray tracer fluxes using "samples" photons per pixel per spectral quadrature point +raytracing = false + +# compute plane parallel (two-stream) fluxes +plane_parallel = true + +[longwave] +longwave = true + +# compute ray tracer fluxes using 2**samples photons per spectral quadrature point +raytracing = false + +# compute plane parallel fluxes (two-stream or no-scat solution depending on "scattering") +plane_parallel = true -# raytracer (test_rte_rrtmgp_rt_gpu) -sw_two_stream = false -sw_raytracing = true -lw_raytracing = true -independent_column = false -lw_scattering = false -cloud_mie = false -single_gpt = false -profiling = false -min_mfp_grid_ratio = true -tica = false - -[ints] -sw_raytracing = 256 -lw_raytracing = 22 -single_gpt = 1 - -[floats] -min_mfp_grid_ratio = 1.0 diff --git a/rcemip/rcemip.ini b/rcemip/rcemip.ini index c81c2e4..0229057 100644 --- a/rcemip/rcemip.ini +++ b/rcemip/rcemip.ini @@ -1,32 +1,29 @@ -[switches] -shortwave = true -longwave = true -fluxes = true +[clouds] liquid_cloud_optics = false ice_cloud_optics = false +delta_cloud = false +tilted_columns = false +cloud_mie = false + +[aerosols] aerosol_optics = false -delta_cloud = true delta_aerosol = false -# gpu_only (test_rte_rrtmgp_gpu) -timings = false +[shortwave] +shortwave = true + +# compute ray tracer fluxes using "samples" photons per pixel per spectral quadrature point +raytracing = false + +# compute plane parallel (two-stream) fluxes +plane_parallel = true + +[longwave] +longwave = true + +# compute ray tracer fluxes using 2**samples photons per spectral quadrature point +raytracing = false + +# compute plane parallel fluxes (two-stream or no-scat solution depending on "scattering") +plane_parallel = true -# raytracer (test_rte_rrtmgp_rt_gpu) -sw_two_stream = false -sw_raytracing = true -lw_raytracing = true -independent_column = false -lw_scattering = false -cloud_mie = false -single_gpt = false -profiling = false -min_mfp_grid_ratio = true -tica = false - -[ints] -sw_raytracing = 256 -lw_raytracing = 22 -single_gpt = 1 - -[floats] -min_mfp_grid_ratio = 1.0 diff --git a/rfmip/rfmip.ini b/rfmip/rfmip.ini index c81c2e4..0229057 100644 --- a/rfmip/rfmip.ini +++ b/rfmip/rfmip.ini @@ -1,32 +1,29 @@ -[switches] -shortwave = true -longwave = true -fluxes = true +[clouds] liquid_cloud_optics = false ice_cloud_optics = false +delta_cloud = false +tilted_columns = false +cloud_mie = false + +[aerosols] aerosol_optics = false -delta_cloud = true delta_aerosol = false -# gpu_only (test_rte_rrtmgp_gpu) -timings = false +[shortwave] +shortwave = true + +# compute ray tracer fluxes using "samples" photons per pixel per spectral quadrature point +raytracing = false + +# compute plane parallel (two-stream) fluxes +plane_parallel = true + +[longwave] +longwave = true + +# compute ray tracer fluxes using 2**samples photons per spectral quadrature point +raytracing = false + +# compute plane parallel fluxes (two-stream or no-scat solution depending on "scattering") +plane_parallel = true -# raytracer (test_rte_rrtmgp_rt_gpu) -sw_two_stream = false -sw_raytracing = true -lw_raytracing = true -independent_column = false -lw_scattering = false -cloud_mie = false -single_gpt = false -profiling = false -min_mfp_grid_ratio = true -tica = false - -[ints] -sw_raytracing = 256 -lw_raytracing = 22 -single_gpt = 1 - -[floats] -min_mfp_grid_ratio = 1.0 From 9cf5c590387b064abe0b5f841fc464f56878cfd0 Mon Sep 17 00:00:00 2001 From: Menno Date: Thu, 5 Mar 2026 12:39:33 +0100 Subject: [PATCH 10/11] update test case input settings --- les_cloudfield/test.ini | 125 +++++++++++++++++++++++----------------- 1 file changed, 71 insertions(+), 54 deletions(-) diff --git a/les_cloudfield/test.ini b/les_cloudfield/test.ini index f35ae36..40d58df 100644 --- a/les_cloudfield/test.ini +++ b/les_cloudfield/test.ini @@ -1,77 +1,94 @@ -[switches] -##### General switches +[clouds] +liquid_cloud_optics = true +ice_cloud_optics = true +delta_cloud = false +tilted_columns = false +cloud_mie = false + +[aerosols] +aerosol_optics = false +delta_aerosol = false + +[shortwave] shortwave = true -longwave = true -fluxes = true -cloud-optics = true -delta-cloud = false +# compute ray tracer fluxes using "samples" photons per pixel per spectral quadrature point +raytracing = true +samples = 256 -aerosol-optics = true -delta-aerosol = false +# ray tracing in independent column mode +rt_independent_column = false -##### test_rte_rrtmgp only (CPU & GPU) -output-optical = false -output-bnd-fluxes = false +# compute plane parallel (two-stream) fluxes +plane_parallel = false -##### test_rte_rrtmgp only (GPU) -timings = false +# only compute specified g-point +single_gpt = -1 -##### Forward raytracer only +[longwave] +longwave = true -# switch on ray tracer. Specifiy ray counts under [ints]: samples per pixel for sw, for lw provide n to use 2^n samples -sw-raytracing = true -lw-raytracing = true +# compute ray tracer fluxes using 2**samples photons per spectral quadrature point +raytracing = true +samples = 22 -# Also perform two-stream computations for shortwave. -sw-two-stream = false +# ray tracing in independent column mode +rt_independent_column = false -# Independent column mode -independent-column = false -# Tilted independent column mode -tica = false +# enable scattering +scattering = false -# Enable scattering for longwave. If enabled, 1D solver will also switch from no-scat solution to two-stream formulation -lw-scattering = true +# compute plane parallel fluxes (two-stream or no-scat solution depending on "scattering") +plane_parallel = false -# Only perform ray tracing when the minimum gasous mean free path is larger than (min-mfp-grid-ratio * grid spacing) and use a 1D solutions otherwise. -# Specify threshold ratio under [floats] -min-mfp-grid-ratio = true +# only compute specified g-point +single_gpt = -1 -# liq and ice flags are automatically switch on when cloud-optics is set to true -# to run either liquid or ice cloud optics -liq-cloud-optics = false -ice-cloud-optics = false +# threshold ratio between lowest gaseus mean free path and horizontal grid spacing +# if this ratio falls below the thershold for a g-point, plane parallel fluxes are used for that g-point instead of the raytracer +min_mf_grid_ratio = 0.0 -##### Common raytracer (forward & backward) switches -cloud-mie = false +[backward] +# enable backward ray tracing using "samples" photons per pixel" +bw_raytracing = true +samples = 16 -# Save optical properties and fluxes for a single g-point. Specify under [ints] -single-gpt = false +# raytrace bands that fall within visible spectrum and convert radiances to XYZ tristimulus values +image = true -# Enable profiling -profiling = false +# compute broadband radiances +broadband = true -##### Backward raytracer only -raytracing = true +# output additional cloud information for each camera pixel +cloud_cam = true -# Read variable "surface_type" to determine scattering type and spectral albedo for each surface grid cell. -# Set either to 0 (water; spectrally constant albedo; specular reflection) or a number between 1 (soil-like spectral albedo; lambertian) and 2 (grass-line spectral albedo; lambertian). -lu-albedo = true +# read surface type from "land_use_map" input variable and use corresponding spectral albedo and scattering type (lambertian, specular) +lu_albedo = true -# broadband radiance computation (sum over all g-points) -broadband = true +[solar_angles] +# solar zenith angle in degrees. mu0 from input file is used if sza < 0 +sza = 45 -# image computation: run only g-points in visible spectrum and convert to XYZ tristimulus values -image = true +# solar azimuth angle in degrees. azi from input file is used if azi < 0 +azi = 135 + +[camera] +cam_field_of_view = 50.0 + +# camera type: (0) fish eye camera, (1) rectangular camera, (2) top-of-atmosphere upwelling radiances +cam_type = 1 -[ints] -single-gpt = 1 +# x,y,z positions of virtual camera +cam_px = 0.0 +cam_py = 0.0 +cam_pz = 2500.0 -sw-raytracing = 256 -lw-raytracing = 22 +# width, height (pixels) of virtual camera +cam_nx = 512 +cam_ny = 512 -raytracing = 128 +# yaw, pitch and roll angles (degrees) of the virtual camera +cam_yaw = 135.0 +cam_pitch = 55.0 +cam_roll = 0.0 -[floats] -min-mfp-grid-ratio = 1.0 From 0cfd44bc78fd6a0835d3091c9c195269cdc8fb86 Mon Sep 17 00:00:00 2001 From: Menno Date: Thu, 5 Mar 2026 12:50:15 +0100 Subject: [PATCH 11/11] fix allsky ini --- allsky/allsky.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/allsky/allsky.ini b/allsky/allsky.ini index 8990a87..5277cbd 100644 --- a/allsky/allsky.ini +++ b/allsky/allsky.ini @@ -1,7 +1,7 @@ [clouds] liquid_cloud_optics = true ice_cloud_optics = true -delta_cloud = false +delta_cloud = true tilted_columns = false cloud_mie = false