diff --git a/data/mie_lut_broadband.nc b/data/mie_lut_broadband.nc index 45e90ed..ab349d8 100644 Binary files a/data/mie_lut_broadband.nc and b/data/mie_lut_broadband.nc differ diff --git a/include_test/Radiation_solver_bw.h b/include_test/Radiation_solver_bw.h index c76b86e..3072a1b 100644 --- a/include_test/Radiation_solver_bw.h +++ b/include_test/Radiation_solver_bw.h @@ -145,8 +145,10 @@ class Radiation_solver_shortwave Array& sw_bnd_flux_dn_dir, Array& sw_bnd_flux_net) const; void load_mie_tables( - const std::string& file_name_mie, - const bool switch_broadband); + const std::string& file_name_mie_bb, + const std::string& file_name_mie_im, + const bool switch_broadband, + const bool switch_image); #ifdef __CUDACC__ void solve_gpu( @@ -244,11 +246,16 @@ class Radiation_solver_shortwave std::unique_ptr cloud_optical_props; std::unique_ptr aerosol_optical_props; - Array_gpu mie_phase_angs; + Array_gpu mie_phase_angs_bb; + Array_gpu mie_cdfs_bb; + Array_gpu mie_angs_bb; + Array_gpu mie_phase_bb; + Array_gpu mie_phase_angs; Array_gpu mie_cdfs; Array_gpu mie_angs; Array_gpu mie_phase; + #endif }; #endif diff --git a/src_test/Radiation_solver_bw.cu b/src_test/Radiation_solver_bw.cu index 1e1d8df..6d10d4d 100644 --- a/src_test/Radiation_solver_bw.cu +++ b/src_test/Radiation_solver_bw.cu @@ -733,31 +733,39 @@ Radiation_solver_shortwave::Radiation_solver_shortwave( } void Radiation_solver_shortwave::load_mie_tables( - const std::string& file_name_mie, - const bool switch_broadband) + const std::string& file_name_mie_bb, + const std::string& file_name_mie_im, + const bool switch_broadband, + const bool switch_image) { - Netcdf_file mie_nc(file_name_mie, Netcdf_mode::Read); - const int n_bnd_sw = this->get_n_bnd_gpu(); - const int n_re = mie_nc.get_dimension_size("r_eff"); - const int n_mie = mie_nc.get_dimension_size("n_ang"); - if (switch_broadband) { - Array mie_cdf(mie_nc.get_variable("phase_cdf", {n_bnd_sw, 1, n_mie}), {n_mie, 1, n_bnd_sw}); - Array mie_ang(mie_nc.get_variable("phase_cdf_angle", {n_bnd_sw, 1, n_re, n_mie}), {n_mie, n_re, 1, n_bnd_sw}); - - Array mie_phase(mie_nc.get_variable("phase", {n_bnd_sw, 1, n_re, n_mie}), {n_mie, n_re, 1, n_bnd_sw}); + Netcdf_file mie_nc(file_name_mie_bb, Netcdf_mode::Read); + const int n_bnd_sw = this->get_n_bnd_gpu(); + const int n_re = mie_nc.get_dimension_size("r_eff"); + const int n_mie = mie_nc.get_dimension_size("n_ang"); + + Array mie_cdf(mie_nc.get_variable("phase_cdf", {n_bnd_sw, n_mie}), {n_mie, 1, n_bnd_sw}); + Array mie_ang(mie_nc.get_variable("phase_cdf_angle", {n_bnd_sw, n_re, n_mie}), {n_mie, n_re, 1, n_bnd_sw}); + + Array mie_phase(mie_nc.get_variable("phase", {n_bnd_sw, n_re, n_mie}), {n_mie, n_re, 1, n_bnd_sw}); Array mie_phase_ang(mie_nc.get_variable("phase_angle", {n_mie}), {n_mie}); - this->mie_cdfs = mie_cdf; - this->mie_angs = mie_ang; + this->mie_cdfs_bb = mie_cdf; + this->mie_angs_bb = mie_ang; - this->mie_phase = mie_phase; - this->mie_phase_angs = mie_phase_ang; + this->mie_phase_bb = mie_phase; + this->mie_phase_angs_bb = mie_phase_ang; } - else + + if (switch_image) { + Netcdf_file mie_nc(file_name_mie_im, Netcdf_mode::Read); + const int n_bnd_sw = this->get_n_bnd_gpu(); + const int n_re = mie_nc.get_dimension_size("r_eff"); + const int n_mie = mie_nc.get_dimension_size("n_ang"); const int n_sub = mie_nc.get_dimension_size("sub_band"); + Array mie_cdf(mie_nc.get_variable("phase_cdf", {n_bnd_sw, n_sub, n_mie}), {n_mie, n_sub, n_bnd_sw}); Array mie_ang(mie_nc.get_variable("phase_cdf_angle", {n_bnd_sw, n_sub, n_re, n_mie}), {n_mie, n_re, n_sub, n_bnd_sw}); @@ -1126,9 +1134,9 @@ void Radiation_solver_shortwave::solve_gpu_bb( const int n_gpt = this->kdist_gpu->get_ngpt(); const int n_bnd = this->kdist_gpu->get_nband(); - const int n_mie = (switch_cloud_mie) ? this->mie_angs.dim(1) : 0; - const int n_re = (switch_cloud_mie) ? this->mie_angs.dim(2) : 0; - + const int n_mie = (switch_cloud_mie) ? this->mie_angs_bb.dim(1) : 0; + const int n_re = (switch_cloud_mie) ? this->mie_angs_bb.dim(2) : 0; + const int cam_nx = radiance.dim(1); const int cam_ny = radiance.dim(2); const Bool top_at_1 = p_lay({1, 1}) < p_lay({1, n_lay}); @@ -1285,9 +1293,9 @@ void Radiation_solver_shortwave::solve_gpu_bb( if (switch_cloud_mie) { - mie_cdfs_sub = mie_cdfs.subset({{ {1, n_mie}, {1,1}, {band, band} }}); - mie_angs_sub = mie_angs.subset({{ {1, n_mie}, {1, n_re}, {1,1}, {band, band} }}); - mie_phase_sub = mie_phase.subset({{ {1, n_mie}, {1, n_re}, {1,1}, {band, band} }}); + mie_cdfs_sub = mie_cdfs_bb.subset({{ {1, n_mie}, {1,1}, {band, band} }}); + mie_angs_sub = mie_angs_bb.subset({{ {1, n_mie}, {1, n_re}, {1,1}, {band, band} }}); + mie_phase_sub = mie_phase_bb.subset({{ {1, n_mie}, {1, n_re}, {1,1}, {band, band} }}); } raytracer.trace_rays_bb( @@ -1298,7 +1306,7 @@ void Radiation_solver_shortwave::solve_gpu_bb( mie_cdfs_sub, mie_angs_sub, mie_phase_sub, - mie_phase_angs, + mie_phase_angs_bb, rel, dynamic_cast(*optical_props).get_tau(), dynamic_cast(*optical_props).get_ssa(), diff --git a/src_test/test_rte_rrtmgp_bw.cu b/src_test/test_rte_rrtmgp_bw.cu index 00ab3fa..4929e12 100644 --- a/src_test/test_rte_rrtmgp_bw.cu +++ b/src_test/test_rte_rrtmgp_bw.cu @@ -245,7 +245,8 @@ void solve_radiation(int argc, char** argv) {"output-optical" , { false, "Enable output of optical properties." }}, {"output-bnd-fluxes", { false, "Enable output of band fluxes." }}, {"lu-albedo" , { false, "Compute spectral albedo from land use map" }}, - {"broadband" , { false, "Compute broadband fluxes" }}, + {"image" , { true, "Compute XYZ values to generate RGB images" }}, + {"broadband" , { false, "Compute broadband radiances" }}, {"profiling" , { false, "Perform additional profiling run." }}, {"delta-cloud" , { false, "delta-scaling of cloud optical properties" }}, {"delta-aerosol" , { false, "delta-scaling of aerosol optical properties" }}, @@ -265,12 +266,13 @@ void solve_radiation(int argc, char** argv) const bool switch_output_optical = command_line_options.at("output-optical" ).first; const bool switch_output_bnd_fluxes = command_line_options.at("output-bnd-fluxes").first; const bool switch_lu_albedo = command_line_options.at("lu-albedo" ).first; + const bool switch_image = command_line_options.at("image" ).first; const bool switch_broadband = command_line_options.at("broadband" ).first; const bool switch_profiling = command_line_options.at("profiling" ).first; const bool switch_delta_cloud = command_line_options.at("delta-cloud" ).first; const bool switch_delta_aerosol = command_line_options.at("delta-aerosol" ).first; - const bool switch_cloud_cam = command_line_options.at("cloud-cam" ).first; - const bool switch_raytracing = command_line_options.at("raytracing" ).first; + const bool switch_cloud_cam = command_line_options.at("cloud-cam" ).first; + const bool switch_raytracing = command_line_options.at("raytracing" ).first; if (switch_longwave) { @@ -682,18 +684,16 @@ void solve_radiation(int argc, char** argv) if (switch_broadband) { radiance.set_dims({camera.nx, camera.ny}); - - if (switch_cloud_mie) - rad_sw.load_mie_tables("mie_lut_broadband.nc", switch_broadband); } - else + if (switch_image) { XYZ.set_dims({camera.nx, camera.ny, 3}); - - if (switch_cloud_mie) - rad_sw.load_mie_tables("mie_lut_visualisation.nc", switch_broadband); } + if (switch_cloud_mie) + rad_sw.load_mie_tables("mie_lut_broadband.nc", "mie_lut_visualisation.nc", switch_broadband, switch_image); + + Array_gpu liwp_cam; Array_gpu tauc_cam; Array_gpu dist_cam; @@ -781,7 +781,7 @@ void solve_radiation(int argc, char** argv) cudaEventDestroy(start); cudaEventDestroy(stop); - Status::print_message("Duration shortwave solver: " + std::to_string(duration) + " (ms)"); + Status::print_message("Duration shortwave solver (broadband version): " + std::to_string(duration) + " (ms)"); }; auto run_solver = [&](const bool tune_step) @@ -856,7 +856,7 @@ void solve_radiation(int argc, char** argv) cudaEventDestroy(start); cudaEventDestroy(stop); - Status::print_message("Duration shortwave solver: " + std::to_string(duration) + " (ms)"); + Status::print_message("Duration shortwave solver (image version): " + std::to_string(duration) + " (ms)"); }; if (switch_broadband) @@ -871,7 +871,7 @@ void solve_radiation(int argc, char** argv) cudaProfilerStop(); } } - else + if (switch_image) { // tune step run_solver(true); @@ -893,30 +893,26 @@ void solve_radiation(int argc, char** argv) if (switch_raytracing) { + output_nc.add_dimension("gpt_sw", n_gpt_sw); + output_nc.add_dimension("band_sw", n_bnd_sw); + + 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_broadband) { Array radiance_cpu(radiance); - output_nc.add_dimension("gpt_sw", n_gpt_sw); - output_nc.add_dimension("band_sw", n_bnd_sw); - - 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}); auto nc_var = output_nc.add_variable("radiance", {"y","x"}); nc_var.insert(radiance_cpu.v(), {0, 0}); nc_var.add_attribute("long_name", "shortwave radiance"); nc_var.add_attribute("units", "W m-2 sr-1"); } - else + if (switch_image) { Array xyz_cpu(XYZ); - output_nc.add_dimension("gpt_sw", n_gpt_sw); - output_nc.add_dimension("band_sw", n_bnd_sw); output_nc.add_dimension("n",3); - 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}); - auto nc_xyz = output_nc.add_variable("XYZ", {"n","y","x"}); nc_xyz.insert(xyz_cpu.v(), {0, 0, 0});