diff --git a/allsky/allsky.ini b/allsky/allsky.ini index 85de65b4..8e1fea55 100644 --- a/allsky/allsky.ini +++ b/allsky/allsky.ini @@ -2,34 +2,31 @@ shortwave = true longwave = true fluxes = true -cloud-optics = true -aerosol-optics = false -output-optical = false -output-bnd-fluxes = false -delta-cloud = true -delta-aerosol = false +liquid_cloud_optics = true +ice_cloud_optics = true +aerosol_optics = false +delta_cloud = true +delta_aerosol = false -# gpu-only (test_rte_rrtmgp_gpu) +# gpu_only (test_rte_rrtmgp_gpu) timings = false # raytracer (test_rte_rrtmgp_rt_gpu) -sw-two-stream = false -sw-raytracing = true -lw-raytracing = true -independent-column = false -liq-cloud-optics = false -ice-cloud-optics = false -lw-scattering = false -cloud-mie = false -single-gpt = false +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 +min_mfp_grid_ratio = true tica = false [ints] -sw-raytracing = 256 -lw-raytracing = 22 -single-gpt = 1 +sw_raytracing = 256 +lw_raytracing = 22 +single_gpt = 1 [floats] -min-mfp-grid-ratio = 1.0 +min_mfp_grid_ratio = 1.0 diff --git a/data/aerosol_optics.nc b/data/aerosol_optics.nc deleted file mode 100644 index 719889f3..00000000 Binary files a/data/aerosol_optics.nc and /dev/null differ diff --git a/include_rt/raytracer_functions.h b/include_rt/raytracer_functions.h index dc5d40e2..07eb7517 100644 --- a/include_rt/raytracer_functions.h +++ b/include_rt/raytracer_functions.h @@ -106,11 +106,11 @@ namespace Raytracer_functions } __device__ - inline Float mie_sample_angle(const Float* mie_cdf, const Float* mie_lut, const Float random_number, const Float r_eff, const int n_mie) + inline Float mie_sample_angle(const Float* mie_cdf, const Float* mie_lut, const Float random_number, const Float* r_eff, const int ijk, const int n_mie) { // interpolation over effective radius. Currently, r_eff should range between 2.5 and 21.5 (similar to RRTMGP) - const int r_idx = min(max(int(r_eff-2.5), 0), 18); - const Float r_rest = fmod(r_eff-Float(2.5),Float(1.)); + const int r_idx = min(max(int(r_eff[ijk]-2.5), 0), 18); + const Float r_rest = fmod(r_eff[ijk]-Float(2.5),Float(1.)); const int i = min(max(0, find_index(mie_cdf, n_mie, random_number)), n_mie - 2); @@ -125,11 +125,11 @@ namespace Raytracer_functions } __device__ - inline Float mie_interpolate_phase_table(const Float* mie_phase, const Float* mie_lut, const Float scat_ang, const Float r_eff, const int n_mie) + inline Float mie_interpolate_phase_table(const Float* mie_phase, const Float* mie_lut, const Float scat_ang, const Float* r_eff, const int ijk, const int n_mie) { // interpolation over effective radius. Currently, r_eff should range between 2.5 and 21.5 (similar to RRTMGP) - const int r_idx = min(max(int(r_eff-2.5), 0), 18); - const Float r_rest = fmod(r_eff-Float(2.5),Float(1.)); + const int r_idx = min(max(int(r_eff[ijk]-2.5), 0), 18); + const Float r_rest = fmod(r_eff[ijk]-Float(2.5),Float(1.)); // interpolation between 1800 equally spaced scattering angles between 0 and PI (both inclusive). constexpr Float d_pi = Float(1.74629942e-03); diff --git a/include_rt_kernels/raytracer_kernels_lw.h b/include_rt_kernels/raytracer_kernels_lw.h index cc34dbc5..b25a51c5 100644 --- a/include_rt_kernels/raytracer_kernels_lw.h +++ b/include_rt_kernels/raytracer_kernels_lw.h @@ -18,10 +18,10 @@ constexpr int rt_lw_kernel_grid = 256; constexpr Float k_null_gas_min = Float(1.e-3); +template __global__ void ray_tracer_lw_kernel( const Int rng_offset, - const bool independent_column, const Int photons_to_shoot, const double* __restrict__ alias_prob, const int* __restrict__ alias_idx, diff --git a/include_rt_kernels/raytracer_kernels_sw.h b/include_rt_kernels/raytracer_kernels_sw.h index b1a88574..21243364 100644 --- a/include_rt_kernels/raytracer_kernels_sw.h +++ b/include_rt_kernels/raytracer_kernels_sw.h @@ -18,9 +18,8 @@ constexpr int rt_kernel_grid = 256; constexpr Float k_null_gas_min = Float(1.e-3); -__global__ +template __global__ void ray_tracer_kernel( - const bool independent_column, const Int photons_to_shoot, const Int qrng_grid_x, const Int qrng_grid_y, diff --git a/include_test/radiation_solver.h b/include_test/radiation_solver.h index 7d311270..ded98f77 100644 --- a/include_test/radiation_solver.h +++ b/include_test/radiation_solver.h @@ -46,8 +46,6 @@ class Radiation_solver_longwave void solve( const bool switch_fluxes, const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const Gas_concs& gas_concs, const Array& p_lay, const Array& p_lev, const Array& t_lay, const Array& t_lev, @@ -55,10 +53,7 @@ class Radiation_solver_longwave const Array& t_sfc, const Array& emis_sfc, const Array& lwp, const Array& iwp, const Array& rel, const Array& dei, - Array& tau, Array& lay_source, - Array& lev_source, Array& sfc_source, - Array& lw_flux_up, Array& lw_flux_dn, Array& lw_flux_net, - Array& lw_bnd_flux_up, Array& lw_bnd_flux_dn, Array& lw_bnd_flux_net) const; + Array& lw_flux_up, Array& lw_flux_dn, Array& lw_flux_net) const; int get_n_gpt() const { return this->kdist->get_ngpt(); }; int get_n_bnd() const { return this->kdist->get_nband(); }; @@ -73,8 +68,6 @@ class Radiation_solver_longwave void solve_gpu( const bool switch_fluxes, const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const Gas_concs_gpu& gas_concs, const Array_gpu& p_lay, const Array_gpu& p_lev, const Array_gpu& t_lay, const Array_gpu& t_lev, @@ -82,10 +75,7 @@ class Radiation_solver_longwave const Array_gpu& t_sfc, const Array_gpu& emis_sfc, const Array_gpu& lwp, const Array_gpu& iwp, const Array_gpu& rel, const Array_gpu& dei, - Array_gpu& tau, 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_bnd_flux_up, Array_gpu& lw_bnd_flux_dn, Array_gpu& lw_bnd_flux_net); + Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net); int get_n_gpt_gpu() const { return this->kdist_gpu->get_ngpt(); }; int get_n_bnd_gpu() const { return this->kdist_gpu->get_nband(); }; @@ -140,8 +130,6 @@ class Radiation_solver_shortwave const bool switch_fluxes, const bool switch_cloud_optics, const bool switch_aerosol_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const bool switch_delta_cloud, const bool switch_delta_aerosol, const Gas_concs& gas_concs, @@ -154,12 +142,8 @@ class Radiation_solver_shortwave const Array& rel, const Array& dei, const Array& rh, const Aerosol_concs& aerosol_concs, - Array& tau, Array& ssa, Array& g, - Array& toa_src, Array& sw_flux_up, Array& sw_flux_dn, - Array& sw_flux_dn_dir, Array& sw_flux_net, - Array& sw_bnd_flux_up, Array& sw_bnd_flux_dn, - Array& sw_bnd_flux_dn_dir, Array& sw_bnd_flux_net) const; + Array& sw_flux_dn_dir, Array& sw_flux_net) const; int get_n_gpt() const { return this->kdist->get_ngpt(); }; int get_n_bnd() const { return this->kdist->get_nband(); }; @@ -177,8 +161,6 @@ class Radiation_solver_shortwave const bool switch_fluxes, const bool switch_cloud_optics, const bool switch_aerosol_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const bool switch_delta_cloud, const bool switch_delta_aerosol, const Gas_concs_gpu& gas_concs, @@ -191,12 +173,8 @@ class Radiation_solver_shortwave const Array_gpu& rel, const Array_gpu& dei, const Array_gpu& rh, const Aerosol_concs_gpu& aerosol_concs, - Array_gpu& tau, Array_gpu& ssa, Array_gpu& g, - Array_gpu& toa_src, Array_gpu& sw_flux_up, Array_gpu& sw_flux_dn, - Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net, - Array_gpu& sw_bnd_flux_up, Array_gpu& sw_bnd_flux_dn, - Array_gpu& sw_bnd_flux_dn_dir, Array_gpu& sw_bnd_flux_net); + Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net); int get_n_gpt_gpu() const { return this->kdist_gpu->get_ngpt(); }; int get_n_bnd_gpu() const { return this->kdist_gpu->get_nband(); }; diff --git a/include_test/radiation_solver_bw.h b/include_test/radiation_solver_bw.h index 6a31a76c..a57677dc 100644 --- a/include_test/radiation_solver_bw.h +++ b/include_test/radiation_solver_bw.h @@ -34,35 +34,14 @@ -class Radiation_solver_longwave +class Radiation_solver_bw_longwave { public: - Radiation_solver_longwave( - const Gas_concs& gas_concs, - const std::string& file_name_gas, - const std::string& file_name_cloud); - - Radiation_solver_longwave( + Radiation_solver_bw_longwave( const Gas_concs_gpu& gas_concs, const std::string& file_name_gas, const std::string& file_name_cloud); - void solve( - const bool switch_fluxes, - const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, - const Gas_concs& gas_concs, - const Array& p_lay, const Array& p_lev, - const Array& t_lay, const Array& t_lev, - const Array& col_dry, - const Array& t_sfc, const Array& emis_sfc, - const Array& lwp, const Array& iwp, - const Array& rel, const Array& dei, - Array& tau, Array& lay_source, - Array& lev_source_inc, Array& lev_source_dec, Array& sfc_source, - Array& lw_flux_up, Array& lw_flux_dn, Array& lw_flux_net, - Array& lw_bnd_flux_up, Array& lw_bnd_flux_dn, Array& lw_bnd_flux_net) const; #ifdef __CUDACC__ void solve_gpu( @@ -108,42 +87,15 @@ class Radiation_solver_longwave }; -class Radiation_solver_shortwave +class Radiation_solver_bw_shortwave { public: - Radiation_solver_shortwave( - const Gas_concs& gas_concs, - const std::string& file_name_gas, - const std::string& file_name_cloud, - const std::string& file_name_aerosol); - Radiation_solver_shortwave( + Radiation_solver_bw_shortwave( const Gas_concs_gpu& gas_concs, const std::string& file_name_gas, const std::string& file_name_cloud, const std::string& file_name_aerosol); - void solve( - const bool switch_fluxes, - const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, - const bool switch_delta_cloud, - const bool switch_delta_aerosol, - const Gas_concs& gas_concs, - const Array& p_lay, const Array& p_lev, - const Array& t_lay, const Array& t_lev, - const Array& col_dry, - const Array& sfc_alb, - const Array& tsi_scaling, const Array& mu0, - const Array& lwp, const Array& iwp, - const Array& rel, const Array& dei, - Array& tau, Array& ssa, Array& g, - Array& toa_src, - Array& sw_flux_up, Array& sw_flux_dn, - Array& sw_flux_dn_dir, Array& sw_flux_net, - Array& sw_bnd_flux_up, Array& sw_bnd_flux_dn, - Array& sw_bnd_flux_dn_dir, Array& sw_bnd_flux_net) const; - void load_mie_tables( const std::string& file_name_mie_bb, const std::string& file_name_mie_im, @@ -152,7 +104,6 @@ class Radiation_solver_shortwave #ifdef __CUDACC__ void solve_gpu( - const bool tune_step, const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, diff --git a/include_test/radiation_solver_rt.h b/include_test/radiation_solver_rt.h index 75d0601e..9df1af22 100644 --- a/include_test/radiation_solver_rt.h +++ b/include_test/radiation_solver_rt.h @@ -50,6 +50,8 @@ class Radiation_solver_longwave 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, diff --git a/include_tilt/tilt_utils.h b/include_tilt/tilt_utils.h index 2d49568f..0d9a81cc 100644 --- a/include_tilt/tilt_utils.h +++ b/include_tilt/tilt_utils.h @@ -121,4 +121,4 @@ void tica_tilt(const Float sza, const Float azi, Array& lwp_out, Array& iwp_out, Array& rel_out, Array& dei_out, Array& rh_out, Gas_concs& gas_concs_out, Aerosol_concs& aerosol_concs_out, std::vector gas_names, std::vector aerosol_names, - bool switch_cloud_optics, bool switch_liq_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics); + bool switch_liquid_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics); diff --git a/rcemip/rcemip.ini b/rcemip/rcemip.ini index d1957694..c81c2e49 100644 --- a/rcemip/rcemip.ini +++ b/rcemip/rcemip.ini @@ -2,34 +2,31 @@ shortwave = true longwave = true fluxes = true -cloud-optics = false -aerosol-optics = false -output-optical = false -output-bnd-fluxes = false -delta-cloud = true -delta-aerosol = false +liquid_cloud_optics = false +ice_cloud_optics = false +aerosol_optics = false +delta_cloud = true +delta_aerosol = false -# gpu-only (test_rte_rrtmgp_gpu) +# gpu_only (test_rte_rrtmgp_gpu) timings = false # raytracer (test_rte_rrtmgp_rt_gpu) -sw-two-stream = false -sw-raytracing = true -lw-raytracing = true -independent-column = false -liq-cloud-optics = false -ice-cloud-optics = false -lw-scattering = false -cloud-mie = false -single-gpt = false +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 +min_mfp_grid_ratio = true tica = false [ints] -sw-raytracing = 256 -lw-raytracing = 22 -single-gpt = 1 +sw_raytracing = 256 +lw_raytracing = 22 +single_gpt = 1 [floats] -min-mfp-grid-ratio = 1.0 +min_mfp_grid_ratio = 1.0 diff --git a/rfmip/rfmip.ini b/rfmip/rfmip.ini index d1957694..c81c2e49 100644 --- a/rfmip/rfmip.ini +++ b/rfmip/rfmip.ini @@ -2,34 +2,31 @@ shortwave = true longwave = true fluxes = true -cloud-optics = false -aerosol-optics = false -output-optical = false -output-bnd-fluxes = false -delta-cloud = true -delta-aerosol = false +liquid_cloud_optics = false +ice_cloud_optics = false +aerosol_optics = false +delta_cloud = true +delta_aerosol = false -# gpu-only (test_rte_rrtmgp_gpu) +# gpu_only (test_rte_rrtmgp_gpu) timings = false # raytracer (test_rte_rrtmgp_rt_gpu) -sw-two-stream = false -sw-raytracing = true -lw-raytracing = true -independent-column = false -liq-cloud-optics = false -ice-cloud-optics = false -lw-scattering = false -cloud-mie = false -single-gpt = false +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 +min_mfp_grid_ratio = true tica = false [ints] -sw-raytracing = 256 -lw-raytracing = 22 -single-gpt = 1 +sw_raytracing = 256 +lw_raytracing = 22 +single_gpt = 1 [floats] -min-mfp-grid-ratio = 1.0 +min_mfp_grid_ratio = 1.0 diff --git a/src_cuda_rt/raytracer_lw.cu b/src_cuda_rt/raytracer_lw.cu index 977031d1..c8346fca 100644 --- a/src_cuda_rt/raytracer_lw.cu +++ b/src_cuda_rt/raytracer_lw.cu @@ -317,23 +317,44 @@ void Raytracer_lw::trace_rays( const Int rng_offset = igpt*rt_lw_kernel_grid*rt_lw_kernel_block; - ray_tracer_lw_kernel<<>>( - rng_offset, - switch_independent_column, - photons_per_thread, - alias_prob.ptr(), - alias_idx.ptr(), - n_alias, - k_null_grid.ptr(), - tod_dn_count.ptr(), - tod_up_count.ptr(), - surface_dn_count.ptr(), - surface_up_count.ptr(), - atmos_count.ptr(), - k_ext.ptr(), scat_asy.ptr(), - emis_sfc.ptr(), - grid_size, grid_d, grid_cells, kn_grid); + if (switch_independent_column) + { + ray_tracer_lw_kernel<<>>( + rng_offset, + photons_per_thread, + alias_prob.ptr(), + alias_idx.ptr(), + n_alias, + k_null_grid.ptr(), + tod_dn_count.ptr(), + tod_up_count.ptr(), + surface_dn_count.ptr(), + surface_up_count.ptr(), + atmos_count.ptr(), + k_ext.ptr(), scat_asy.ptr(), + emis_sfc.ptr(), + grid_size, grid_d, grid_cells, kn_grid); + } + else + { + ray_tracer_lw_kernel<<>>( + rng_offset, + photons_per_thread, + alias_prob.ptr(), + alias_idx.ptr(), + n_alias, + k_null_grid.ptr(), + tod_dn_count.ptr(), + tod_up_count.ptr(), + surface_dn_count.ptr(), + surface_up_count.ptr(), + atmos_count.ptr(), + k_ext.ptr(), scat_asy.ptr(), + emis_sfc.ptr(), + grid_size, grid_d, grid_cells, kn_grid); + + } const Float power_per_photon = Float(total_power / (photons_per_thread * rt_lw_kernel_grid * rt_lw_kernel_block)); count_to_flux_2d<<>>( diff --git a/src_cuda_rt/raytracer_sw.cu b/src_cuda_rt/raytracer_sw.cu index e241f75f..ab40fadb 100644 --- a/src_cuda_rt/raytracer_sw.cu +++ b/src_cuda_rt/raytracer_sw.cu @@ -375,30 +375,58 @@ void Raytracer::trace_rays( const int mie_table_size = mie_cdf.size(); const int qrng_gpt_offset = (igpt-1) * rt_kernel_grid_size * rt_kernel_block_size * photons_per_thread; - ray_tracer_kernel<<>>( - switch_independent_column, - photons_per_thread, - qrng_grid_x, - qrng_grid_y, - qrng_gpt_offset, - k_null_grid.ptr(), - tod_dn_count.ptr(), - tod_up_count.ptr(), - surface_down_direct_count.ptr(), - surface_down_diffuse_count.ptr(), - surface_up_count.ptr(), - atmos_direct_count.ptr(), - atmos_diffuse_count.ptr(), - k_ext.ptr(), scat_asy.ptr(), - r_eff.ptr(), - tod_inc_direct, - tod_inc_diffuse, - surface_albedo.ptr(), - grid_size, grid_d, grid_cells, kn_grid, - sun_direction, - this->qrng_vectors_gpu, this->qrng_constants_gpu, - mie_cdf.ptr(), mie_ang.ptr(), mie_table_size); + if (switch_independent_column) + { + ray_tracer_kernel<<>>( + photons_per_thread, + qrng_grid_x, + qrng_grid_y, + qrng_gpt_offset, + k_null_grid.ptr(), + tod_dn_count.ptr(), + tod_up_count.ptr(), + surface_down_direct_count.ptr(), + surface_down_diffuse_count.ptr(), + surface_up_count.ptr(), + atmos_direct_count.ptr(), + atmos_diffuse_count.ptr(), + k_ext.ptr(), scat_asy.ptr(), + r_eff.ptr(), + tod_inc_direct, + tod_inc_diffuse, + surface_albedo.ptr(), + grid_size, grid_d, grid_cells, kn_grid, + sun_direction, + this->qrng_vectors_gpu, this->qrng_constants_gpu, + mie_cdf.ptr(), mie_ang.ptr(), mie_table_size); + } + else + { + ray_tracer_kernel<<>>( + photons_per_thread, + qrng_grid_x, + qrng_grid_y, + qrng_gpt_offset, + k_null_grid.ptr(), + tod_dn_count.ptr(), + tod_up_count.ptr(), + surface_down_direct_count.ptr(), + surface_down_diffuse_count.ptr(), + surface_up_count.ptr(), + atmos_direct_count.ptr(), + atmos_diffuse_count.ptr(), + k_ext.ptr(), scat_asy.ptr(), + r_eff.ptr(), + tod_inc_direct, + tod_inc_diffuse, + surface_albedo.ptr(), + grid_size, grid_d, grid_cells, kn_grid, + sun_direction, + this->qrng_vectors_gpu, this->qrng_constants_gpu, + mie_cdf.ptr(), mie_ang.ptr(), mie_table_size); + + } // convert counts to fluxes const Float toa_src = tod_inc_direct + tod_inc_diffuse; diff --git a/src_kernels_cuda_rt/raytracer_kernels_bw.cu b/src_kernels_cuda_rt/raytracer_kernels_bw.cu index 053739d6..dc60442a 100644 --- a/src_kernels_cuda_rt/raytracer_kernels_bw.cu +++ b/src_kernels_cuda_rt/raytracer_kernels_bw.cu @@ -275,7 +275,8 @@ namespace const Float sun_solid_angle, const Float g, const Float* __restrict__ mie_phase_ang, const Float* __restrict__ mie_phase, - const Float r_eff, + const Float* __restrict__ r_eff, + const int ijk, const int mie_table_size, const Vector& normal, const Phase_kind kind) @@ -289,7 +290,7 @@ namespace else if (kind == Phase_kind::Mie) { // return interpolate_mie_phase_table(mie_phase_ang, mie_phase, max(0.05, acos(cos_angle)), r_eff, mie_table_size) * sun_solid_angle; - return mie_interpolate_phase_table(mie_phase_ang, mie_phase, acos(cos_angle), r_eff, mie_table_size) * sun_solid_angle; + return mie_interpolate_phase_table(mie_phase_ang, mie_phase, acos(cos_angle), r_eff, ijk, mie_table_size) * sun_solid_angle; } else if (kind == Phase_kind::Rayleigh) { @@ -497,7 +498,7 @@ void ray_tracer_kernel_bw( // direct contribution const Phase_kind kind = (scatter_type==0) ? Phase_kind::Rayleigh : Phase_kind::HG; - const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, g, mie_phase_ang_shared, mie_phase, Float(0.), 0, surface_normal, kind); + const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, g, mie_phase_ang_shared, mie_phase, r_eff, 0, 0, surface_normal, kind); const Float trans_sun = transmission_direct_sun(photon,n,rng,sun_direction, k_null_grid,k_ext, bg_tau_cum, z_lev_bg, bg_idx, @@ -585,7 +586,7 @@ void ray_tracer_kernel_bw( : Phase_kind::Lambertian; // SUN SCATTERING GOES HERE - const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, Float(0.), mie_phase_ang_shared, mie_phase, Float(0.), 0, + const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, Float(0.), mie_phase_ang_shared, mie_phase, r_eff, 0, 0, surface_normal, surface_kind); const Float trans_sun = transmission_direct_sun(photon,n,rng,sun_direction, k_null_grid,k_ext, @@ -714,7 +715,7 @@ void ray_tracer_kernel_bw( } const Float cos_scat = scatter_type == 0 ? rayleigh(rng()) : // gases -> rayleigh, 1 ? ( (mie_cdf_table_size > 0) //clouds: Mie or HG - ? cos( mie_sample_angle(mie_cdf_shared, mie_ang, rng(), r_eff[ijk], mie_cdf_table_size) ) + ? cos( mie_sample_angle(mie_cdf_shared, mie_ang, rng(), r_eff, ijk, mie_cdf_table_size) ) : henyey(g, rng())) : henyey(g, rng()); //aerosols const Float sin_scat = max(Float(0.), sqrt(Float(1.) - cos_scat*cos_scat + Float_epsilon)); @@ -725,8 +726,7 @@ void ray_tracer_kernel_bw( ? Phase_kind::Mie : Phase_kind::HG : Phase_kind::HG; - const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, g, mie_phase_ang_shared, mie_phase, r_eff[ijk], mie_phase_table_size, - surface_normal, kind); + const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, g, mie_phase_ang_shared, mie_phase, r_eff, ijk, mie_phase_table_size, surface_normal, kind); const Float trans_sun = transmission_direct_sun(photon,n,rng,sun_direction, k_null_grid,k_ext, bg_tau_cum, z_lev_bg, bg_idx, diff --git a/src_kernels_cuda_rt/raytracer_kernels_lw.cu b/src_kernels_cuda_rt/raytracer_kernels_lw.cu index 7dcd9686..5976d45f 100644 --- a/src_kernels_cuda_rt/raytracer_kernels_lw.cu +++ b/src_kernels_cuda_rt/raytracer_kernels_lw.cu @@ -124,10 +124,9 @@ namespace } -__global__ -void ray_tracer_lw_kernel( +template +__global__ void ray_tracer_lw_kernel( const Int rng_offset, - const bool independent_column, const Int photons_to_shoot, const double* __restrict__ alias_prob, const int* __restrict__ alias_idx, @@ -442,3 +441,31 @@ void ray_tracer_lw_kernel( } } } + +template __global__ void ray_tracer_lw_kernel( + const Int, const Int, const double*, const int*, + int, const Float*, Float*, Float*, + Float*, + Float*, + Float*, + const Float*, + const Optics_scat*, + const Float*, + const Vector, + const Vector, + const Vector, + const Vector); + +template __global__ void ray_tracer_lw_kernel( + const Int, const Int, const double*, const int*, + int, const Float*, Float*, Float*, + Float*, + Float*, + Float*, + const Float*, + const Optics_scat*, + const Float*, + const Vector, + const Vector, + const Vector, + const Vector); diff --git a/src_kernels_cuda_rt/raytracer_kernels_sw.cu b/src_kernels_cuda_rt/raytracer_kernels_sw.cu index 18ed7d3a..9663499a 100644 --- a/src_kernels_cuda_rt/raytracer_kernels_sw.cu +++ b/src_kernels_cuda_rt/raytracer_kernels_sw.cu @@ -119,9 +119,8 @@ namespace } -__global__ +template __global__ void ray_tracer_kernel( - const bool independent_column, const Int photons_to_shoot, const Int qrng_grid_x, const Int qrng_grid_y, @@ -411,7 +410,7 @@ void ray_tracer_kernel( // 0 (gas): rayleigh, 1 (cloud): mie if mie_table_size>0 else HG, 2 (aerosols) HG const Float cos_scat = scatter_type == 0 ? rayleigh(rng()) : // gases -> rayleigh, 1 ? ( (mie_table_size > 0) //clouds: Mie or HG - ? cos( mie_sample_angle(mie_cdf_shared, mie_ang, rng(), r_eff[ijk], mie_table_size) ) + ? cos( mie_sample_angle(mie_cdf_shared, mie_ang, rng(), r_eff, ijk, mie_table_size) ) : henyey(g, rng())) : henyey(g, rng()); //aerosols const Float sin_scat = max(Float(0.), sqrt(Float(1.) - cos_scat*cos_scat + Float_epsilon)); @@ -458,3 +457,15 @@ void ray_tracer_kernel( } } } + +template __global__ void ray_tracer_kernel( + const Int, const Int, const Int, const Int,const Float*, Float*, Float*, Float*, Float*, + Float*, Float*, Float*, const Float*, const Optics_scat*, const Float*, const Float, const Float, + const Float*, const Vector, const Vector, const Vector, const Vector, + const Vector, curandDirectionVectors32_t*, unsigned int*, const Float*, const Float*, const int); + +template __global__ void ray_tracer_kernel( + const Int, const Int, const Int, const Int,const Float*, Float*, Float*, Float*, Float*, + Float*, Float*, Float*, const Float*, const Optics_scat*, const Float*, const Float, const Float, + const Float*, const Vector, const Vector, const Vector, const Vector, + const Vector, curandDirectionVectors32_t*, unsigned int*, const Float*, const Float*, const int); diff --git a/src_test/CMakeLists.txt b/src_test/CMakeLists.txt index 16342d04..6ad4a779 100644 --- a/src_test/CMakeLists.txt +++ b/src_test/CMakeLists.txt @@ -28,11 +28,11 @@ if(USECUDA) add_executable(test_rte_rrtmgp radiation_solver.cu test_rte_rrtmgp.cu) target_link_libraries(test_rte_rrtmgp rte_rrtmgp rte_rrtmgp_cuda ${LIBS} m) - add_executable(test_rte_rrtmgp_rt radiation_solver_rt.cu test_rte_rrtmgp_rt.cu ../src_tilt/tilt_utils.cpp) + add_executable(test_rte_rrtmgp_rt radiation_solver_rt.cu test_rte_rrtmgp_rt.cu radiation_solver_bw.cu ../src_tilt/tilt_utils.cpp) target_link_libraries(test_rte_rrtmgp_rt rte_rrtmgp rte_rrtmgp_cuda rte_rrtmgp_cuda_rt curand ${LIBS} m) - add_executable(test_rte_rrtmgp_bw radiation_solver_bw.cu test_rte_rrtmgp_bw.cu) - target_link_libraries(test_rte_rrtmgp_bw rte_rrtmgp rte_rrtmgp_cuda rte_rrtmgp_cuda_rt curand ${LIBS} m) + #add_executable(test_rte_rrtmgp_bw radiation_solver_bw.cu test_rte_rrtmgp_bw.cu) + #target_link_libraries(test_rte_rrtmgp_bw rte_rrtmgp rte_rrtmgp_cuda rte_rrtmgp_cuda_rt curand ${LIBS} m) add_executable(test_rt_lite test_rt_lite.cu) target_link_libraries(test_rt_lite rte_rrtmgp_cuda rte_rrtmgp_cuda_rt curand ${LIBS} m) diff --git a/src_test/radiation_solver.cpp b/src_test/radiation_solver.cpp index 7ae67f91..1fa0bcfa 100644 --- a/src_test/radiation_solver.cpp +++ b/src_test/radiation_solver.cpp @@ -384,8 +384,6 @@ Radiation_solver_longwave::Radiation_solver_longwave( void Radiation_solver_longwave::solve( const bool switch_fluxes, const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const Gas_concs& gas_concs, const Array& p_lay, const Array& p_lev, const Array& t_lay, const Array& t_lev, @@ -393,10 +391,7 @@ void Radiation_solver_longwave::solve( const Array& t_sfc, const Array& emis_sfc, const Array& lwp, const Array& iwp, const Array& rel, const Array& dei, - Array& tau, Array& lay_source, - Array& lev_source, Array& sfc_source, - Array& lw_flux_up, Array& lw_flux_dn, Array& lw_flux_net, - Array& lw_bnd_flux_up, Array& lw_bnd_flux_dn, Array& lw_bnd_flux_net) const + Array& lw_flux_up, Array& lw_flux_dn, Array& lw_flux_net) const { const int n_col = p_lay.dim(1); const int n_lay = p_lay.dim(2); @@ -445,8 +440,7 @@ void Radiation_solver_longwave::solve( std::unique_ptr& cloud_optical_props_subset_in, Source_func_lw& sources_subset_in, const Array& emis_sfc_subset_in, - Fluxes_broadband& fluxes, - Fluxes_broadband& bnd_fluxes) + Fluxes_broadband& fluxes) { const int n_col_in = col_e_in - col_s_in + 1; Gas_concs gas_concs_subset(gas_concs, col_s_in, n_col_in); @@ -479,52 +473,20 @@ void Radiation_solver_longwave::solve( dei.subset({{ {col_s_in, col_e_in}, {1, n_lay} }}), *cloud_optical_props_subset_in); - // cloud->delta_scale(); - // Add the cloud optical props to the gas optical properties. add_to( dynamic_cast(*optical_props_subset_in), dynamic_cast(*cloud_optical_props_subset_in)); } - // Store the optical properties, if desired. - if (switch_output_optical) - { - for (int igpt=1; igpt<=n_gpt; ++igpt) - for (int ilay=1; ilay<=n_lay; ++ilay) - for (int icol=1; icol<=n_col_in; ++icol) - { - tau ({icol+col_s_in-1, ilay, igpt}) = optical_props_subset_in->get_tau()({icol, ilay, igpt}); - lay_source({icol+col_s_in-1, ilay, igpt}) = sources_subset_in.get_lay_source()({icol, ilay, igpt}); - lev_source({icol+col_s_in-1, ilay, igpt}) = sources_subset_in.get_lev_source()({icol, ilay, igpt}); - } - - for (int igpt=1; igpt<=n_gpt; ++igpt) - for (int icol=1; icol<=n_col_in; ++icol) - lev_source({icol+col_s_in-1, n_lev, igpt}) = sources_subset_in.get_lev_source()({icol, n_lev, igpt}); - - for (int igpt=1; igpt<=n_gpt; ++igpt) - for (int icol=1; icol<=n_col_in; ++icol) - sfc_source({icol+col_s_in-1, igpt}) = sources_subset_in.get_sfc_source()({icol, igpt}); - } - if (!switch_fluxes) return; Array gpt_flux_up; Array gpt_flux_dn; - // Save the output per gpt if postprocessing is desired. - if (switch_output_bnd_fluxes) - { - gpt_flux_up.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn.set_dims({n_col_in, n_lev, n_gpt}); - } - else - { - gpt_flux_up.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); - } + gpt_flux_up.set_dims({n_col_in, n_lev, 1}); + gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); constexpr int n_ang = 1; @@ -537,42 +499,14 @@ void Radiation_solver_longwave::solve( gpt_flux_up, gpt_flux_dn, n_ang); - if (switch_output_bnd_fluxes) - { - // Aggegated fluxes. - fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - lw_flux_up ({icol+col_s_in-1, ilev}) = fluxes.get_flux_up ()({icol, ilev}); - lw_flux_dn ({icol+col_s_in-1, ilev}) = fluxes.get_flux_dn ()({icol, ilev}); - lw_flux_net({icol+col_s_in-1, ilev}) = fluxes.get_flux_net()({icol, ilev}); - } - - // Aggegated fluxes per band - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - for (int ibnd=1; ibnd<=n_bnd; ++ibnd) - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - lw_bnd_flux_up ({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_up ()({icol, ilev, ibnd}); - lw_bnd_flux_dn ({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_dn ()({icol, ilev, ibnd}); - lw_bnd_flux_net({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_net()({icol, ilev, ibnd}); - } - } - else - { - // Copy the data to the output. - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - lw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up({icol, ilev, 1}); - lw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}); - lw_flux_net({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); - } - } + // Copy the data to the output. + for (int ilev=1; ilev<=n_lev; ++ilev) + for (int icol=1; icol<=n_col_in; ++icol) + { + lw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up({icol, ilev, 1}); + lw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}); + lw_flux_net({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); + } }; for (int b=1; b<=n_blocks; ++b) @@ -584,8 +518,6 @@ void Radiation_solver_longwave::solve( std::unique_ptr fluxes_subset = std::make_unique(n_col_block, n_lev); - std::unique_ptr bnd_fluxes_subset = - std::make_unique(n_col_block, n_lev, n_bnd); call_kernels( col_s, col_e, @@ -593,8 +525,7 @@ void Radiation_solver_longwave::solve( cloud_optical_props_subset, *sources_subset, emis_sfc_subset, - *fluxes_subset, - *bnd_fluxes_subset); + *fluxes_subset); } if (n_col_block_residual > 0) @@ -605,8 +536,6 @@ void Radiation_solver_longwave::solve( Array emis_sfc_residual = emis_sfc.subset({{ {1, n_bnd}, {col_s, col_e} }}); std::unique_ptr fluxes_residual = std::make_unique(n_col_block_residual, n_lev); - std::unique_ptr bnd_fluxes_residual = - std::make_unique(n_col_block_residual, n_lev, n_bnd); call_kernels( col_s, col_e, @@ -614,8 +543,7 @@ void Radiation_solver_longwave::solve( cloud_optical_props_residual, *sources_residual, emis_sfc_residual, - *fluxes_residual, - *bnd_fluxes_residual); + *fluxes_residual); } } @@ -646,8 +574,6 @@ void Radiation_solver_shortwave::solve( const bool switch_fluxes, const bool switch_cloud_optics, const bool switch_aerosol_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const bool switch_delta_cloud, const bool switch_delta_aerosol, const Gas_concs& gas_concs, @@ -660,12 +586,8 @@ void Radiation_solver_shortwave::solve( const Array& rel, const Array& dei, const Array& rh, const Aerosol_concs& aerosol_concs, - Array& tau, Array& ssa, Array& g, - Array& toa_src, Array& sw_flux_up, Array& sw_flux_dn, - Array& sw_flux_dn_dir, Array& sw_flux_net, - Array& sw_bnd_flux_up, Array& sw_bnd_flux_dn, - Array& sw_bnd_flux_dn_dir, Array& sw_bnd_flux_net) const + Array& sw_flux_dn_dir, Array& sw_flux_net) const { const int n_col = p_lay.dim(1); const int n_lay = p_lay.dim(2); @@ -714,8 +636,7 @@ void Radiation_solver_shortwave::solve( std::unique_ptr& optical_props_subset_in, std::unique_ptr& cloud_optical_props_subset_in, std::unique_ptr& aerosol_optical_props_subset_in, - Fluxes_broadband& fluxes, - Fluxes_broadband& bnd_fluxes) + Fluxes_broadband& fluxes) { const int n_col_in = col_e_in - col_s_in + 1; Gas_concs gas_concs_subset(gas_concs, col_s_in, n_col_in); @@ -786,23 +707,6 @@ void Radiation_solver_shortwave::solve( } - // Store the optical properties, if desired. - if (switch_output_optical) - { - for (int igpt=1; igpt<=n_gpt; ++igpt) - for (int ilay=1; ilay<=n_lay; ++ilay) - for (int icol=1; icol<=n_col_in; ++icol) - { - tau({icol+col_s_in-1, ilay, igpt}) = optical_props_subset_in->get_tau()({icol, ilay, igpt}); - ssa({icol+col_s_in-1, ilay, igpt}) = optical_props_subset_in->get_ssa()({icol, ilay, igpt}); - g ({icol+col_s_in-1, ilay, igpt}) = optical_props_subset_in->get_g ()({icol, ilay, igpt}); - } - - for (int igpt=1; igpt<=n_gpt; ++igpt) - for (int icol=1; icol<=n_col_in; ++icol) - toa_src({icol+col_s_in-1, igpt}) = toa_src_subset({icol, igpt}); - } - if (!switch_fluxes) return; @@ -810,19 +714,9 @@ void Radiation_solver_shortwave::solve( Array gpt_flux_dn; Array gpt_flux_dn_dir; - // Save the output per gpt if postprocessing is desired. - if (switch_output_bnd_fluxes) - { - gpt_flux_up.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn_dir.set_dims({n_col_in, n_lev, n_gpt}); - } - else - { - gpt_flux_up.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn_dir.set_dims({n_col_in, n_lev, 1}); - } + gpt_flux_up.set_dims({n_col_in, n_lev, 1}); + gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); + gpt_flux_dn_dir.set_dims({n_col_in, n_lev, 1}); Rte_sw::rte_sw( optical_props_subset_in, @@ -836,44 +730,16 @@ void Radiation_solver_shortwave::solve( gpt_flux_dn, gpt_flux_dn_dir); - if (switch_output_bnd_fluxes) - { - // Aggegated fluxes. - fluxes.reduce(gpt_flux_up, gpt_flux_dn, gpt_flux_dn_dir, optical_props_subset_in, top_at_1); - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - sw_flux_up ({icol+col_s_in-1, ilev}) = fluxes.get_flux_up ()({icol, ilev}); - sw_flux_dn ({icol+col_s_in-1, ilev}) = fluxes.get_flux_dn ()({icol, ilev}); - sw_flux_dn_dir({icol+col_s_in-1, ilev}) = fluxes.get_flux_dn_dir()({icol, ilev}); - sw_flux_net ({icol+col_s_in-1, ilev}) = fluxes.get_flux_net ()({icol, ilev}); - } - - // Aggegated fluxes per band - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, gpt_flux_dn_dir, optical_props_subset_in, top_at_1); - - for (int ibnd=1; ibnd<=n_bnd; ++ibnd) - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - sw_bnd_flux_up ({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_up ()({icol, ilev, ibnd}); - sw_bnd_flux_dn ({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_dn ()({icol, ilev, ibnd}); - sw_bnd_flux_dn_dir({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_dn_dir()({icol, ilev, ibnd}); - sw_bnd_flux_net ({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_net ()({icol, ilev, ibnd}); - } - } - else - { - // Copy the data to the output. - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - sw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up ({icol, ilev, 1}); - sw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn ({icol, ilev, 1}); - sw_flux_dn_dir({icol+col_s_in-1, ilev}) = gpt_flux_dn_dir({icol, ilev, 1}); - sw_flux_net ({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); - } - } + // Copy the data to the output. + for (int ilev=1; ilev<=n_lev; ++ilev) + for (int icol=1; icol<=n_col_in; ++icol) + { + sw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up ({icol, ilev, 1}); + sw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn ({icol, ilev, 1}); + sw_flux_dn_dir({icol+col_s_in-1, ilev}) = gpt_flux_dn_dir({icol, ilev, 1}); + sw_flux_net ({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); + } + }; for (int b=1; b<=n_blocks; ++b) @@ -883,16 +749,13 @@ void Radiation_solver_shortwave::solve( std::unique_ptr fluxes_subset = std::make_unique(n_col_block, n_lev); - std::unique_ptr bnd_fluxes_subset = - std::make_unique(n_col_block, n_lev, n_bnd); call_kernels( col_s, col_e, optical_props_subset, cloud_optical_props_subset, aerosol_optical_props_subset, - *fluxes_subset, - *bnd_fluxes_subset); + *fluxes_subset); } if (n_col_block_residual > 0) @@ -902,15 +765,12 @@ void Radiation_solver_shortwave::solve( std::unique_ptr fluxes_residual = std::make_unique(n_col_block_residual, n_lev); - std::unique_ptr bnd_fluxes_residual = - std::make_unique(n_col_block_residual, n_lev, n_bnd); call_kernels( col_s, col_e, optical_props_residual, cloud_optical_props_residual, aerosol_optical_props_residual, - *fluxes_residual, - *bnd_fluxes_residual); + *fluxes_residual); } } diff --git a/src_test/radiation_solver.cu b/src_test/radiation_solver.cu index 3e240be8..75c3fdd1 100644 --- a/src_test/radiation_solver.cu +++ b/src_test/radiation_solver.cu @@ -419,8 +419,6 @@ Radiation_solver_longwave::Radiation_solver_longwave( void Radiation_solver_longwave::solve_gpu( const bool switch_fluxes, const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const Gas_concs_gpu& gas_concs, const Array_gpu& p_lay, const Array_gpu& p_lev, const Array_gpu& t_lay, const Array_gpu& t_lev, @@ -428,10 +426,7 @@ void Radiation_solver_longwave::solve_gpu( const Array_gpu& t_sfc, const Array_gpu& emis_sfc, const Array_gpu& lwp, const Array_gpu& iwp, const Array_gpu& rel, const Array_gpu& dei, - Array_gpu& tau, 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_bnd_flux_up, Array_gpu& lw_bnd_flux_dn, Array_gpu& lw_bnd_flux_net) + Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net) { const int n_col = p_lay.dim(1); const int n_lay = p_lay.dim(2); @@ -472,8 +467,7 @@ void Radiation_solver_longwave::solve_gpu( std::unique_ptr& cloud_optical_props_subset_in, Source_func_lw_gpu& sources_subset_in, const Array_gpu& emis_sfc_subset_in, - Fluxes_broadband_gpu& fluxes, - Fluxes_broadband_gpu& bnd_fluxes) + Fluxes_broadband_gpu& fluxes) { const int n_col_in = col_e_in - col_s_in + 1; Gas_concs_gpu gas_concs_subset(gas_concs, col_s_in, n_col_in); @@ -506,51 +500,18 @@ void Radiation_solver_longwave::solve_gpu( dei.subset({{ {col_s_in, col_e_in}, {1, n_lay} }}), *cloud_optical_props_subset_in); - //cloud->delta_scale(); - // Add the cloud optical props to the gas optical properties. add_to( dynamic_cast(*optical_props_subset_in), dynamic_cast(*cloud_optical_props_subset_in)); } - // Store the optical properties, if desired. - if (switch_output_optical) - { - Subset_kernels_cuda::get_from_subset( - n_col, n_lay, n_gpt, n_col_in, col_s_in, tau.ptr(), lay_source.ptr(), lev_source.ptr(), - optical_props_subset_in->get_tau().ptr(), sources_subset_in.get_lay_source().ptr(), - sources_subset_in.get_lev_source().ptr()); - - Subset_kernels_cuda::get_from_subset( - n_col, n_gpt, n_col_in, col_s_in, sfc_source.ptr(), sources_subset_in.get_sfc_source().ptr()); - } - if (!switch_fluxes) return; Array_gpu gpt_flux_up({n_col_in, n_lev, n_gpt}); Array_gpu gpt_flux_dn({n_col_in, n_lev, n_gpt}); - - // CvH The structure below is valid if broadband flux solvers are implemented - /* - Array_gpu gpt_flux_up; - Array_gpu gpt_flux_dn; - - // Save the output per gpt if postprocessing is desired. - if (switch_output_bnd_fluxes) - { - gpt_flux_up.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn.set_dims({n_col_in, n_lev, n_gpt}); - } - else - { - gpt_flux_up.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); - } - */ - constexpr int n_ang = 1; rte_lw.rte_lw( @@ -570,47 +531,6 @@ void Radiation_solver_longwave::solve_gpu( fluxes.get_flux_up().ptr(), fluxes.get_flux_dn().ptr(), fluxes.get_flux_net().ptr()); - if (switch_output_bnd_fluxes) - { - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_bnd, n_col_in, col_s_in, lw_bnd_flux_up.ptr(), lw_bnd_flux_dn.ptr(), lw_bnd_flux_net.ptr(), - bnd_fluxes.get_bnd_flux_up().ptr(), bnd_fluxes.get_bnd_flux_dn().ptr(), bnd_fluxes.get_bnd_flux_net().ptr()); - - } - - // CvH The structure below is valid if broadband flux solvers are implemented - /* - if (switch_output_bnd_fluxes) - { - // Aggegated fluxes. - fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - // Copy the data to the output. - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_col_in, col_s_in, 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()); - - // Aggegated fluxes per band - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_bnd, n_col_in, col_s_in, lw_bnd_flux_up.ptr(), lw_bnd_flux_dn.ptr(), lw_bnd_flux_net.ptr(), - bnd_fluxes.get_bnd_flux_up().ptr(), bnd_fluxes.get_bnd_flux_dn().ptr(), bnd_fluxes.get_bnd_flux_net().ptr()); - } - else - { - // Copy the data to the output. - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - lw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up({icol, ilev, 1}); - lw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}); - lw_flux_net({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); - } - } - */ }; for (int b=1; b<=n_blocks; ++b) @@ -622,8 +542,6 @@ void Radiation_solver_longwave::solve_gpu( std::unique_ptr fluxes_subset = std::make_unique(n_col_block, n_lev); - std::unique_ptr bnd_fluxes_subset = - std::make_unique(n_col_block, n_lev, n_bnd); call_kernels( col_s, col_e, @@ -631,8 +549,7 @@ void Radiation_solver_longwave::solve_gpu( cloud_optical_props_subset, *sources_subset, emis_sfc_subset, - *fluxes_subset, - *bnd_fluxes_subset); + *fluxes_subset); } if (n_col_block_residual > 0) @@ -643,8 +560,6 @@ void Radiation_solver_longwave::solve_gpu( Array_gpu emis_sfc_residual = emis_sfc.subset({{ {1, n_bnd}, {col_s, col_e} }}); std::unique_ptr fluxes_residual = std::make_unique(n_col_block_residual, n_lev); - std::unique_ptr bnd_fluxes_residual = - std::make_unique(n_col_block_residual, n_lev, n_bnd); call_kernels( col_s, col_e, @@ -652,8 +567,7 @@ void Radiation_solver_longwave::solve_gpu( cloud_optical_props_residual, *sources_residual, emis_sfc_residual, - *fluxes_residual, - *bnd_fluxes_residual); + *fluxes_residual); } } @@ -684,8 +598,6 @@ void Radiation_solver_shortwave::solve_gpu( const bool switch_fluxes, const bool switch_cloud_optics, const bool switch_aerosol_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const bool switch_delta_cloud, const bool switch_delta_aerosol, const Gas_concs_gpu& gas_concs, @@ -698,12 +610,8 @@ void Radiation_solver_shortwave::solve_gpu( const Array_gpu& rel, const Array_gpu& dei, const Array_gpu& rh, const Aerosol_concs_gpu& aerosol_concs, - Array_gpu& tau, Array_gpu& ssa, Array_gpu& g, - Array_gpu& toa_src, Array_gpu& sw_flux_up, Array_gpu& sw_flux_dn, - Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net, - Array_gpu& sw_bnd_flux_up, Array_gpu& sw_bnd_flux_dn, - Array_gpu& sw_bnd_flux_dn_dir, Array_gpu& sw_bnd_flux_net) + Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net) { const int n_col = p_lay.dim(1); const int n_lay = p_lay.dim(2); @@ -745,8 +653,7 @@ void Radiation_solver_shortwave::solve_gpu( std::unique_ptr& optical_props_subset_in, std::unique_ptr& cloud_optical_props_subset_in, std::unique_ptr& aerosol_optical_props_subset_in, - Fluxes_broadband_gpu& fluxes, - Fluxes_broadband_gpu& bnd_fluxes) + Fluxes_broadband_gpu& fluxes) { const int n_col_in = col_e_in - col_s_in + 1; Gas_concs_gpu gas_concs_subset(gas_concs, col_s_in, n_col_in); @@ -768,8 +675,10 @@ void Radiation_solver_shortwave::solve_gpu( optical_props_subset_in, toa_src_subset, col_dry_subset); + auto tsi_scaling_subset = tsi_scaling.subset({{ {col_s_in, col_e_in} }}); scaling_to_subset(n_col_in, n_gpt, toa_src_subset, tsi_scaling_subset); + if (switch_cloud_optics) { Array cld_mask_liq({n_col_in, n_lay}); @@ -809,16 +718,6 @@ void Radiation_solver_shortwave::solve_gpu( dynamic_cast(*aerosol_optical_props_subset_in)); } - // Store the optical properties, if desired. - if (switch_output_optical) - { - Subset_kernels_cuda::get_from_subset( - n_col, n_lay, n_gpt, n_col_in, col_s_in, tau.ptr(), ssa.ptr(), g.ptr(), optical_props_subset_in->get_tau().ptr(), - optical_props_subset_in->get_ssa().ptr(), optical_props_subset_in->get_g().ptr()); - - Subset_kernels_cuda::get_from_subset( - n_col, n_gpt, n_col_in, col_s_in, toa_src.ptr(), toa_src_subset.ptr()); - } if (!switch_fluxes) return; @@ -827,26 +726,6 @@ void Radiation_solver_shortwave::solve_gpu( Array_gpu gpt_flux_dn({n_col_in, n_lev, n_gpt}); Array_gpu gpt_flux_dn_dir({n_col_in, n_lev, n_gpt}); - // CvH The structure below is valid if broadband flux solvers are implemented - /* - Array_gpu gpt_flux_up; - Array_gpu gpt_flux_dn; - Array_gpu gpt_flux_dn_dir; - - if (switch_output_bnd_fluxes) - { - gpt_flux_up.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn_dir.set_dims({n_col_in, n_lev, n_gpt}); - } - else - { - gpt_flux_up.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn_dir.set_dims({n_col_in, n_lev, 1}); - } - */ - rte_sw.rte_sw( optical_props_subset_in, top_at_1, @@ -866,50 +745,8 @@ void Radiation_solver_shortwave::solve_gpu( n_col, n_lev, n_col_in, col_s_in, 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_output_bnd_fluxes) - { - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_bnd, n_col_in, col_s_in, sw_bnd_flux_up.ptr(), sw_bnd_flux_dn.ptr(), sw_bnd_flux_dn_dir.ptr(), sw_bnd_flux_net.ptr(), - bnd_fluxes.get_bnd_flux_up().ptr(), bnd_fluxes.get_bnd_flux_dn().ptr(), bnd_fluxes.get_bnd_flux_dn_dir().ptr(), bnd_fluxes.get_bnd_flux_net().ptr()); - } - - // CvH The structure below is valid if broadband flux solvers are implemented - /* - if (switch_output_bnd_fluxes) - { - // Aggegated fluxes. - fluxes.reduce(gpt_flux_up, gpt_flux_dn, gpt_flux_dn_dir, optical_props_subset_in, top_at_1); - - // Copy the data to the output. - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_col_in, col_s_in, sw_flux_up, sw_flux_dn, sw_flux_dn_dir, sw_flux_net, - fluxes.get_flux_up(), fluxes.get_flux_dn(), fluxes.get_flux_dn_dir(), fluxes.get_flux_net()); - - // Aggegated fluxes per band - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_bnd, n_col_in, col_s_in, sw_bnd_flux_up, sw_bnd_flux_dn, sw_bnd_flux_dn_dir, sw_bnd_flux_net, - bnd_fluxes.get_bnd_flux_up(), bnd_fluxes.get_bnd_flux_dn(), bnd_fluxes.get_bnd_flux_dn_dir(), bnd_fluxes.get_bnd_flux_net()); - } - else - { - // Copy the data to the output. - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - sw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up ({icol, ilev, 1}); - sw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn ({icol, ilev, 1}); - sw_flux_dn_dir({icol+col_s_in-1, ilev}) = gpt_flux_dn_dir({icol, ilev, 1}); - sw_flux_net ({icol+col_s_in-1, ilev}) = gpt_flux_dn ({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); - } - } - */ }; - for (int b=1; b<=n_blocks; ++b) { const int col_s = (b-1) * n_col_block + 1; @@ -917,16 +754,13 @@ void Radiation_solver_shortwave::solve_gpu( std::unique_ptr fluxes_subset = std::make_unique(n_col_block, n_lev); - std::unique_ptr bnd_fluxes_subset = - std::make_unique(n_col_block, n_lev, n_bnd); call_kernels( col_s, col_e, optical_props_subset, cloud_optical_props_subset, aerosol_optical_props_subset, - *fluxes_subset, - *bnd_fluxes_subset); + *fluxes_subset); } if (n_col_block_residual > 0) @@ -936,15 +770,12 @@ void Radiation_solver_shortwave::solve_gpu( std::unique_ptr fluxes_residual = std::make_unique(n_col_block_residual, n_lev); - std::unique_ptr bnd_fluxes_residual = - std::make_unique(n_col_block_residual, n_lev, n_bnd); call_kernels( col_s, col_e, optical_props_residual, cloud_optical_props_residual, aerosol_optical_props_residual, - *fluxes_residual, - *bnd_fluxes_residual); + *fluxes_residual); } } diff --git a/src_test/radiation_solver_bw.cu b/src_test/radiation_solver_bw.cu index 94a4be71..2f6570b0 100644 --- a/src_test/radiation_solver_bw.cu +++ b/src_test/radiation_solver_bw.cu @@ -477,7 +477,7 @@ namespace } -Radiation_solver_longwave::Radiation_solver_longwave( +Radiation_solver_bw_longwave::Radiation_solver_bw_longwave( const Gas_concs_gpu& gas_concs, const std::string& file_name_gas, const std::string& file_name_cloud) @@ -491,7 +491,7 @@ Radiation_solver_longwave::Radiation_solver_longwave( } -void Radiation_solver_longwave::solve_gpu( +void Radiation_solver_bw_longwave::solve_gpu( const bool switch_fluxes, const bool switch_cloud_optics, const bool switch_output_optical, @@ -508,7 +508,7 @@ void Radiation_solver_longwave::solve_gpu( Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net, Array_gpu& lw_bnd_flux_up, Array_gpu& lw_bnd_flux_dn, Array_gpu& lw_bnd_flux_net) { - throw std::runtime_error("Longwave raytracing is not implemented"); + throw std::runtime_error("Longwave backward raytracing is not implemented"); /* const int n_col = p_lay.dim(1); @@ -720,7 +720,7 @@ Float xyz_irradiance( } -Radiation_solver_shortwave::Radiation_solver_shortwave( +Radiation_solver_bw_shortwave::Radiation_solver_bw_shortwave( const Gas_concs_gpu& gas_concs, const std::string& file_name_gas, const std::string& file_name_cloud, @@ -737,7 +737,7 @@ Radiation_solver_shortwave::Radiation_solver_shortwave( load_and_init_aerosol_optics(file_name_aerosol)); } -void Radiation_solver_shortwave::load_mie_tables( +void Radiation_solver_bw_shortwave::load_mie_tables( const std::string& file_name_mie_bb, const std::string& file_name_mie_im, const bool switch_broadband, @@ -790,8 +790,7 @@ void Radiation_solver_shortwave::load_mie_tables( } -void Radiation_solver_shortwave::solve_gpu( - const bool tune_step, +void Radiation_solver_bw_shortwave::solve_gpu( const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, @@ -886,7 +885,7 @@ void Radiation_solver_shortwave::solve_gpu( Array_gpu albedo; if (!switch_lu_albedo) albedo = sfc_alb.subset({{ {band, band}, {1, n_col}}}); - if (!tune_step && (! (band == 10 || band == 11 || band ==12))) continue; + if (! (band == 10 || band == 11 || band ==12)) continue; const Float solar_source_band = kdist_gpu->band_source(band_limits_gpt({1,band}), band_limits_gpt({2,band})); @@ -929,8 +928,6 @@ void Radiation_solver_shortwave::solve_gpu( gas_optics_subset(col_s, n_col_residual); } - if (tune_step) return; - constexpr bool do_scattering = true; toa_src.fill(toa_src_temp({1}) * tsi_scaling({1})); @@ -1112,7 +1109,7 @@ void Radiation_solver_shortwave::solve_gpu( } } -void Radiation_solver_shortwave::solve_gpu_bb( +void Radiation_solver_bw_shortwave::solve_gpu_bb( const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, diff --git a/src_test/radiation_solver_rt.cu b/src_test/radiation_solver_rt.cu index 9f77d26e..7278181a 100644 --- a/src_test/radiation_solver_rt.cu +++ b/src_test/radiation_solver_rt.cu @@ -483,6 +483,8 @@ void Radiation_solver_longwave::solve_gpu( 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, @@ -643,7 +645,9 @@ void Radiation_solver_longwave::solve_gpu( dei, switch_lw_scattering, *cloud_optical_props); - // cloud->delta_scale(); + + if (switch_lw_scattering && switch_delta_cloud) + cloud_optical_props->delta_scale(); } @@ -670,7 +674,9 @@ void Radiation_solver_longwave::solve_gpu( rh, p_lev, switch_lw_scattering, *aerosol_optical_props); - // aerosol->delta_scale(); + + if (switch_lw_scattering && switch_delta_aerosol) + aerosol_optical_props->delta_scale(); } // Add the cloud optical props to the gas optical properties. diff --git a/src_test/test_rte_rrtmgp.cpp b/src_test/test_rte_rrtmgp.cpp index 734a693d..2f8a1aa4 100644 --- a/src_test/test_rte_rrtmgp.cpp +++ b/src_test/test_rte_rrtmgp.cpp @@ -137,15 +137,14 @@ 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_cloud_optics = get_ini_value(settings, "switches", "cloud-optics", false); - const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol-optics", false); - const bool switch_output_optical = get_ini_value(settings, "switches", "output-optical", false); - const bool switch_output_bnd_fluxes = get_ini_value(settings, "switches", "output-bnd-fluxes", 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_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); ////// READ THE ATMOSPHERIC DATA ////// @@ -202,19 +201,35 @@ void solve_radiation(int argc, char** argv) Array rel; Array dei; + const bool switch_cloud_optics = switch_liquid_cloud_optics || switch_ice_cloud_optics; if (switch_cloud_optics) { lwp.set_dims({n_col, n_lay}); - lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); - + rel.set_dims({n_col, n_lay}); iwp.set_dims({n_col, n_lay}); - iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); + dei.set_dims({n_col, n_lay}); - rel.set_dims({n_col, n_lay}); - rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); + if (switch_liquid_cloud_optics) + { + lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); + rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); + } + else + { + lwp.fill(Float(0.)); + rel.fill(Float(0.)); + } - dei.set_dims({n_col, n_lay}); - dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); + if (switch_ice_cloud_optics) + { + iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); + dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); + } + else + { + iwp.fill(Float(0.)); + dei.fill(Float(0.)); + } } Array rh; @@ -270,20 +285,6 @@ void solve_radiation(int argc, char** argv) Array emis_sfc(input_nc.get_variable("emis_sfc", {n_col_y, n_col_x, n_bnd_lw}), {n_bnd_lw, n_col}); Array t_sfc(input_nc.get_variable("t_sfc", {n_col_y, n_col_x}), {n_col}); - // Create output arrays. - Array lw_tau; - Array lay_source; - Array lev_source; - Array sfc_source; - - if (switch_output_optical) - { - lw_tau .set_dims({n_col, n_lay, n_gpt_lw}); - lay_source.set_dims({n_col, n_lay, n_gpt_lw}); - lev_source.set_dims({n_col, n_lev, n_gpt_lw}); - sfc_source.set_dims({n_col, n_gpt_lw}); - } - Array lw_flux_up; Array lw_flux_dn; Array lw_flux_net; @@ -295,18 +296,6 @@ void solve_radiation(int argc, char** argv) lw_flux_net.set_dims({n_col, n_lev}); } - Array lw_bnd_flux_up; - Array lw_bnd_flux_dn; - Array lw_bnd_flux_net; - - if (switch_output_bnd_fluxes) - { - lw_bnd_flux_up .set_dims({n_col, n_lev, n_bnd_lw}); - lw_bnd_flux_dn .set_dims({n_col, n_lev, n_bnd_lw}); - lw_bnd_flux_net.set_dims({n_col, n_lev, n_bnd_lw}); - } - - // Solve the radiation. Status::print_message("Solving the longwave radiation."); @@ -315,8 +304,6 @@ void solve_radiation(int argc, char** argv) rad_lw.solve( switch_fluxes, switch_cloud_optics, - switch_output_optical, - switch_output_bnd_fluxes, gas_concs, p_lay, p_lev, t_lay, t_lev, @@ -324,9 +311,7 @@ void solve_radiation(int argc, char** argv) t_sfc, emis_sfc, lwp, iwp, rel, dei, - lw_tau, lay_source, lev_source, sfc_source, - lw_flux_up, lw_flux_dn, lw_flux_net, - lw_bnd_flux_up, lw_bnd_flux_dn, lw_bnd_flux_net); + lw_flux_up, lw_flux_dn, lw_flux_net); auto time_end = std::chrono::high_resolution_clock::now(); auto duration = std::chrono::duration(time_end-time_start).count(); @@ -343,25 +328,6 @@ 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_output_optical) - { - 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().v(), {0, 0}); - - auto nc_lw_tau = output_nc.add_variable("lw_tau", {"gpt_lw", "lay", "y", "x"}); - nc_lw_tau.insert(lw_tau.v(), {0, 0, 0, 0}); - - auto nc_lay_source = output_nc.add_variable("lay_source", {"gpt_lw", "lay", "y", "x"}); - auto nc_lev_source = output_nc.add_variable("lev_source", {"gpt_lw", "lev", "y", "x"}); - - auto nc_sfc_source = output_nc.add_variable("sfc_source", {"gpt_lw", "y", "x"}); - - nc_lay_source.insert(lay_source.v(), {0, 0, 0, 0}); - nc_lev_source.insert(lev_source.v(), {0, 0, 0, 0}); - - nc_sfc_source.insert(sfc_source.v(), {0, 0, 0}); - } - if (switch_fluxes) { auto nc_lw_flux_up = output_nc.add_variable("lw_flux_up" , {"lev", "y", "x"}); @@ -372,16 +338,6 @@ void solve_radiation(int argc, char** argv) nc_lw_flux_dn .insert(lw_flux_dn .v(), {0, 0, 0}); nc_lw_flux_net.insert(lw_flux_net.v(), {0, 0, 0}); - if (switch_output_bnd_fluxes) - { - auto nc_lw_bnd_flux_up = output_nc.add_variable("lw_bnd_flux_up" , {"band_lw", "lev", "y", "x"}); - auto nc_lw_bnd_flux_dn = output_nc.add_variable("lw_bnd_flux_dn" , {"band_lw", "lev", "y", "x"}); - auto nc_lw_bnd_flux_net = output_nc.add_variable("lw_bnd_flux_net", {"band_lw", "lev", "y", "x"}); - - nc_lw_bnd_flux_up .insert(lw_bnd_flux_up .v(), {0, 0, 0, 0}); - nc_lw_bnd_flux_dn .insert(lw_bnd_flux_dn .v(), {0, 0, 0, 0}); - nc_lw_bnd_flux_net.insert(lw_bnd_flux_net.v(), {0, 0, 0, 0}); - } } } @@ -422,20 +378,6 @@ void solve_radiation(int argc, char** argv) tsi_scaling({icol}) = Float(1.); } - // Create output arrays. - Array sw_tau; - Array ssa; - Array g; - Array toa_source; - - if (switch_output_optical) - { - sw_tau .set_dims({n_col, n_lay, n_gpt_sw}); - ssa .set_dims({n_col, n_lay, n_gpt_sw}); - g .set_dims({n_col, n_lay, n_gpt_sw}); - toa_source.set_dims({n_col, n_gpt_sw}); - } - Array sw_flux_up; Array sw_flux_dn; Array sw_flux_dn_dir; @@ -449,20 +391,6 @@ void solve_radiation(int argc, char** argv) sw_flux_net .set_dims({n_col, n_lev}); } - Array sw_bnd_flux_up; - Array sw_bnd_flux_dn; - Array sw_bnd_flux_dn_dir; - Array sw_bnd_flux_net; - - if (switch_output_bnd_fluxes) - { - sw_bnd_flux_up .set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_dn .set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_dn_dir.set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_net .set_dims({n_col, n_lev, n_bnd_sw}); - } - - // Solve the radiation. Status::print_message("Solving the shortwave radiation."); @@ -472,8 +400,6 @@ void solve_radiation(int argc, char** argv) switch_fluxes, switch_cloud_optics, switch_aerosol_optics, - switch_output_optical, - switch_output_bnd_fluxes, switch_delta_cloud, switch_delta_aerosol, gas_concs, @@ -486,12 +412,8 @@ void solve_radiation(int argc, char** argv) rel, dei, rh, aerosol_concs, - sw_tau, ssa, g, - toa_source, sw_flux_up, sw_flux_dn, - sw_flux_dn_dir, sw_flux_net, - sw_bnd_flux_up, sw_bnd_flux_dn, - sw_bnd_flux_dn_dir, sw_bnd_flux_net); + sw_flux_dn_dir, sw_flux_net); auto time_end = std::chrono::high_resolution_clock::now(); auto duration = std::chrono::duration(time_end-time_start).count(); @@ -508,23 +430,6 @@ 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_output_optical) - { - 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().v(), {0, 0}); - - auto nc_sw_tau = output_nc.add_variable("sw_tau", {"gpt_sw", "lay", "y", "x"}); - auto nc_ssa = output_nc.add_variable("ssa" , {"gpt_sw", "lay", "y", "x"}); - auto nc_g = output_nc.add_variable("g" , {"gpt_sw", "lay", "y", "x"}); - - nc_sw_tau.insert(sw_tau.v(), {0, 0, 0, 0}); - nc_ssa .insert(ssa .v(), {0, 0, 0, 0}); - nc_g .insert(g .v(), {0, 0, 0, 0}); - - auto nc_toa_source = output_nc.add_variable("toa_source", {"gpt_sw", "y", "x"}); - nc_toa_source.insert(toa_source.v(), {0, 0, 0}); - } - if (switch_fluxes) { auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); @@ -548,31 +453,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 (switch_output_bnd_fluxes) - { - auto nc_sw_bnd_flux_up = output_nc.add_variable("sw_bnd_flux_up" , {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_dn = output_nc.add_variable("sw_bnd_flux_dn" , {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_dn_dir = output_nc.add_variable("sw_bnd_flux_dn_dir", {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_net = output_nc.add_variable("sw_bnd_flux_net" , {"band_sw", "lev", "y", "x"}); - - nc_sw_bnd_flux_up .insert(sw_bnd_flux_up .v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_dn .insert(sw_bnd_flux_dn .v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_dn_dir.insert(sw_bnd_flux_dn_dir.v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_net .insert(sw_bnd_flux_net .v(), {0, 0, 0, 0}); - - nc_sw_bnd_flux_up.add_attribute("long_name","Upwelling shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_up.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_dn.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_dn_dir.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_net.add_attribute("long_name","Net shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_net.add_attribute("units","W m-2"); - } } } diff --git a/src_test/test_rte_rrtmgp.cu b/src_test/test_rte_rrtmgp.cu index b385bf65..0bb86ca1 100644 --- a/src_test/test_rte_rrtmgp.cu +++ b/src_test/test_rte_rrtmgp.cu @@ -163,16 +163,14 @@ 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_cloud_optics = get_ini_value(settings, "switches", "cloud-optics", false); - const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol-optics", false); - const bool switch_output_optical = get_ini_value(settings, "switches", "output-optical", false); - const bool switch_output_bnd_fluxes = get_ini_value(settings, "switches", "output-bnd-fluxes", false); - const bool switch_timings = get_ini_value(settings, "switches", "timings", 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_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); ////// READ THE ATMOSPHERIC DATA ////// @@ -185,6 +183,7 @@ void solve_radiation(int argc, char** argv) const int n_lay = input_nc.get_dimension_size("lay"); const int n_lev = input_nc.get_dimension_size("lev"); const int n_col = n_col_x * n_col_y; + // 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}); Array t_lay(input_nc.get_variable("t_lay", {n_lay, n_col_y, n_col_x}), {n_col, n_lay}); @@ -228,19 +227,35 @@ void solve_radiation(int argc, char** argv) Array rel; Array dei; + const bool switch_cloud_optics = switch_liquid_cloud_optics || switch_ice_cloud_optics; if (switch_cloud_optics) { lwp.set_dims({n_col, n_lay}); - lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); - + rel.set_dims({n_col, n_lay}); iwp.set_dims({n_col, n_lay}); - iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); + dei.set_dims({n_col, n_lay}); - rel.set_dims({n_col, n_lay}); - rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); + if (switch_liquid_cloud_optics) + { + lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); + rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); + } + else + { + lwp.fill(Float(0.)); + rel.fill(Float(0.)); + } - dei.set_dims({n_col, n_lay}); - dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); + if (switch_ice_cloud_optics) + { + iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); + dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); + } + else + { + iwp.fill(Float(0.)); + dei.fill(Float(0.)); + } } Array rh; @@ -317,20 +332,6 @@ void solve_radiation(int argc, char** argv) Array emis_sfc(input_nc.get_variable("emis_sfc", {n_col_y, n_col_x, n_bnd_lw}), {n_bnd_lw, n_col}); Array t_sfc(input_nc.get_variable("t_sfc", {n_col_y, n_col_x}), {n_col}); - // Create output arrays. - Array_gpu lw_tau; - Array_gpu lay_source; - Array_gpu lev_source; - Array_gpu sfc_source; - - if (switch_output_optical) - { - lw_tau .set_dims({n_col, n_lay, n_gpt_lw}); - lay_source.set_dims({n_col, n_lay, n_gpt_lw}); - lev_source.set_dims({n_col, n_lay, n_gpt_lw}); - sfc_source.set_dims({n_col, n_gpt_lw}); - } - Array_gpu lw_flux_up; Array_gpu lw_flux_dn; Array_gpu lw_flux_net; @@ -342,18 +343,6 @@ void solve_radiation(int argc, char** argv) lw_flux_net.set_dims({n_col, n_lev}); } - Array_gpu lw_bnd_flux_up; - Array_gpu lw_bnd_flux_dn; - Array_gpu lw_bnd_flux_net; - - if (switch_output_bnd_fluxes) - { - lw_bnd_flux_up .set_dims({n_col, n_lev, n_bnd_lw}); - lw_bnd_flux_dn .set_dims({n_col, n_lev, n_bnd_lw}); - lw_bnd_flux_net.set_dims({n_col, n_lev, n_bnd_lw}); - } - - // Solve the radiation. Status::print_message("Solving the longwave radiation."); @@ -383,8 +372,6 @@ void solve_radiation(int argc, char** argv) rad_lw.solve_gpu( switch_fluxes, switch_cloud_optics, - switch_output_optical, - switch_output_bnd_fluxes, gas_concs_gpu, p_lay_gpu, p_lev_gpu, t_lay_gpu, t_lev_gpu, @@ -392,9 +379,7 @@ void solve_radiation(int argc, char** argv) t_sfc_gpu, emis_sfc_gpu, lwp_gpu, iwp_gpu, rel_gpu, dei_gpu, - lw_tau, lay_source, lev_source, sfc_source, - lw_flux_up, lw_flux_dn, lw_flux_net, - lw_bnd_flux_up, lw_bnd_flux_dn, lw_bnd_flux_net); + lw_flux_up, lw_flux_dn, lw_flux_net); cudaEventRecord(stop, 0); cudaEventSynchronize(stop); @@ -410,30 +395,11 @@ void solve_radiation(int argc, char** argv) // Tuning step; run_solver(); - // Profiling step; - cudaProfilerStart(); - run_solver(); - cudaProfilerStop(); - - if (switch_timings) - { - constexpr int n_measures=10; - for (int n=0; n lw_tau_cpu(lw_tau); - Array lay_source_cpu(lay_source); - Array sfc_source_cpu(sfc_source); - Array lev_source_cpu(lev_source); 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_bnd_flux_up_cpu(lw_bnd_flux_up); - Array lw_bnd_flux_dn_cpu(lw_bnd_flux_dn); - Array lw_bnd_flux_net_cpu(lw_bnd_flux_net); output_nc.add_dimension("gpt_lw", n_gpt_lw); output_nc.add_dimension("band_lw", n_bnd_lw); @@ -441,25 +407,6 @@ 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_output_optical) - { - 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}); - - auto nc_lw_tau = output_nc.add_variable("lw_tau", {"gpt_lw", "lay", "y", "x"}); - nc_lw_tau.insert(lw_tau_cpu.v(), {0, 0, 0, 0}); - - auto nc_lay_source = output_nc.add_variable("lay_source", {"gpt_lw", "lay", "y", "x"}); - auto nc_lev_source = output_nc.add_variable("lev_source", {"gpt_lw", "lev", "y", "x"}); - - auto nc_sfc_source = output_nc.add_variable("sfc_source", {"gpt_lw", "y", "x"}); - - nc_lay_source.insert(lay_source_cpu.v(), {0, 0, 0, 0}); - nc_lev_source.insert(lev_source_cpu.v(), {0, 0, 0, 0}); - - nc_sfc_source.insert(sfc_source_cpu.v(), {0, 0, 0}); - } - if (switch_fluxes) { auto nc_lw_flux_up = output_nc.add_variable("lw_flux_up" , {"lev", "y", "x"}); @@ -469,17 +416,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 (switch_output_bnd_fluxes) - { - auto nc_lw_bnd_flux_up = output_nc.add_variable("lw_bnd_flux_up" , {"band_lw", "lev", "y", "x"}); - auto nc_lw_bnd_flux_dn = output_nc.add_variable("lw_bnd_flux_dn" , {"band_lw", "lev", "y", "x"}); - auto nc_lw_bnd_flux_net = output_nc.add_variable("lw_bnd_flux_net", {"band_lw", "lev", "y", "x"}); - - nc_lw_bnd_flux_up .insert(lw_bnd_flux_up_cpu.v(), {0, 0, 0, 0}); - nc_lw_bnd_flux_dn .insert(lw_bnd_flux_dn_cpu.v(), {0, 0, 0, 0}); - nc_lw_bnd_flux_net.insert(lw_bnd_flux_net_cpu.v(), {0, 0, 0, 0}); - } } } @@ -492,7 +428,8 @@ void solve_radiation(int argc, char** argv) Gas_concs_gpu gas_concs_gpu(gas_concs); - Radiation_solver_shortwave rad_sw(gas_concs_gpu, switch_cloud_optics, switch_aerosol_optics, "coefficients_sw.nc", "cloud_coefficients_sw.nc", "aerosol_optics_sw.nc"); + Radiation_solver_shortwave rad_sw(gas_concs_gpu, switch_cloud_optics, switch_aerosol_optics, + "coefficients_sw.nc", "cloud_coefficients_sw.nc", "aerosol_optics_sw.nc"); // Read the boundary conditions. const int n_bnd_sw = rad_sw.get_n_bnd_gpu(); @@ -522,20 +459,6 @@ void solve_radiation(int argc, char** argv) tsi_scaling({icol}) = Float(1.); } - // Create output arrays. - Array_gpu sw_tau; - Array_gpu ssa; - Array_gpu g; - Array_gpu toa_source; - - if (switch_output_optical) - { - sw_tau .set_dims({n_col, n_lay, n_gpt_sw}); - ssa .set_dims({n_col, n_lay, n_gpt_sw}); - g .set_dims({n_col, n_lay, n_gpt_sw}); - toa_source.set_dims({n_col, n_gpt_sw}); - } - Array_gpu sw_flux_up; Array_gpu sw_flux_dn; Array_gpu sw_flux_dn_dir; @@ -549,20 +472,6 @@ void solve_radiation(int argc, char** argv) sw_flux_net .set_dims({n_col, n_lev}); } - Array_gpu sw_bnd_flux_up; - Array_gpu sw_bnd_flux_dn; - Array_gpu sw_bnd_flux_dn_dir; - Array_gpu sw_bnd_flux_net; - - if (switch_output_bnd_fluxes) - { - sw_bnd_flux_up .set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_dn .set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_dn_dir.set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_net .set_dims({n_col, n_lev, n_bnd_sw}); - } - - // Solve the radiation. Status::print_message("Solving the shortwave radiation."); @@ -597,8 +506,6 @@ void solve_radiation(int argc, char** argv) switch_fluxes, switch_cloud_optics, switch_aerosol_optics, - switch_output_optical, - switch_output_bnd_fluxes, switch_delta_cloud, switch_delta_aerosol, gas_concs_gpu, @@ -611,12 +518,8 @@ void solve_radiation(int argc, char** argv) rel_gpu, dei_gpu, rh_gpu, aerosol_concs_gpu, - sw_tau, ssa, g, - toa_source, sw_flux_up, sw_flux_dn, - sw_flux_dn_dir, sw_flux_net, - sw_bnd_flux_up, sw_bnd_flux_dn, - sw_bnd_flux_dn_dir, sw_bnd_flux_net); + sw_flux_dn_dir, sw_flux_net); cudaEventRecord(stop, 0); cudaEventSynchronize(stop); @@ -632,32 +535,12 @@ void solve_radiation(int argc, char** argv) // Tuning step; run_solver(); - // Profiling step; - cudaProfilerStart(); - run_solver(); - cudaProfilerStop(); - - if (switch_timings) - { - constexpr int n_measures=10; - for (int n=0; n sw_tau_cpu(sw_tau); - Array ssa_cpu(ssa); - Array g_cpu(g); - Array toa_source_cpu(toa_source); Array sw_flux_up_cpu(sw_flux_up); 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_bnd_flux_up_cpu(sw_bnd_flux_up); - Array sw_bnd_flux_dn_cpu(sw_bnd_flux_dn); - Array sw_bnd_flux_dn_dir_cpu(sw_bnd_flux_dn_dir); - Array sw_bnd_flux_net_cpu(sw_bnd_flux_net); output_nc.add_dimension("gpt_sw", n_gpt_sw); output_nc.add_dimension("band_sw", n_bnd_sw); @@ -665,23 +548,6 @@ 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_output_optical) - { - 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}); - - auto nc_sw_tau = output_nc.add_variable("sw_tau", {"gpt_sw", "lay", "y", "x"}); - auto nc_ssa = output_nc.add_variable("ssa" , {"gpt_sw", "lay", "y", "x"}); - auto nc_g = output_nc.add_variable("g" , {"gpt_sw", "lay", "y", "x"}); - - nc_sw_tau.insert(sw_tau_cpu.v(), {0, 0, 0, 0}); - nc_ssa .insert(ssa_cpu .v(), {0, 0, 0, 0}); - nc_g .insert(g_cpu .v(), {0, 0, 0, 0}); - - auto nc_toa_source = output_nc.add_variable("toa_source", {"gpt_sw", "y", "x"}); - nc_toa_source.insert(toa_source_cpu.v(), {0, 0, 0}); - } - if (switch_fluxes) { auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); @@ -706,30 +572,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 (switch_output_bnd_fluxes) - { - auto nc_sw_bnd_flux_up = output_nc.add_variable("sw_bnd_flux_up" , {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_dn = output_nc.add_variable("sw_bnd_flux_dn" , {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_dn_dir = output_nc.add_variable("sw_bnd_flux_dn_dir", {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_net = output_nc.add_variable("sw_bnd_flux_net" , {"band_sw", "lev", "y", "x"}); - - nc_sw_bnd_flux_up .insert(sw_bnd_flux_up_cpu .v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_dn .insert(sw_bnd_flux_dn_cpu .v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_dn_dir.insert(sw_bnd_flux_dn_dir_cpu.v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_net .insert(sw_bnd_flux_net_cpu .v(), {0, 0, 0, 0}); - - nc_sw_bnd_flux_up.add_attribute("long_name","Upwelling shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_up.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_dn.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_dn_dir.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_net.add_attribute("long_name","Net shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_net.add_attribute("units","W m-2"); - } } } diff --git a/src_test/test_rte_rrtmgp_bw.cu b/src_test/test_rte_rrtmgp_bw.cu index 55020210..ccd8e346 100644 --- a/src_test/test_rte_rrtmgp_bw.cu +++ b/src_test/test_rte_rrtmgp_bw.cu @@ -228,24 +228,23 @@ void solve_radiation(int argc, char** argv) const Vector kn_grid = {input_nc.get_variable("ngrid_x"), input_nc.get_variable("ngrid_y"), input_nc.get_variable("ngrid_z")}; - - // Reading camera data - Netcdf_group cam_in = input_nc.get_group("camera-settings"); Camera camera; - camera.fov = cam_in.get_variable("fov"); - camera.cam_type = int(cam_in.get_variable("cam_type")); - camera.position = {cam_in.get_variable("px"), - cam_in.get_variable("py"), - cam_in.get_variable("pz")}; - - camera.nx = int(cam_in.get_variable("nx")); - camera.ny = int(cam_in.get_variable("ny")); + 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); - camera.setup_rotation_matrix(cam_in.get_variable("yaw"), - cam_in.get_variable("pitch"), - cam_in.get_variable("roll")); - camera.setup_normal_camera(camera); + 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}); @@ -386,7 +385,7 @@ void solve_radiation(int argc, char** argv) Gas_concs_gpu gas_concs_gpu(gas_concs); - Radiation_solver_longwave rad_lw(gas_concs_gpu, "coefficients_lw.nc", "cloud_coefficients_lw.nc"); + Radiation_solver_bw_longwave rad_lw(gas_concs_gpu, "coefficients_lw.nc", "cloud_coefficients_lw.nc"); // Read the boundary conditions. const int n_bnd_lw = rad_lw.get_n_bnd_gpu(); @@ -573,7 +572,7 @@ void solve_radiation(int argc, char** argv) Gas_concs_gpu gas_concs_gpu(gas_concs); - Radiation_solver_shortwave rad_sw(gas_concs_gpu, "coefficients_sw.nc", "cloud_coefficients_sw.nc","aerosol_optics_sw.nc"); + Radiation_solver_bw_shortwave rad_sw(gas_concs_gpu, "coefficients_sw.nc", "cloud_coefficients_sw.nc","aerosol_optics_sw.nc"); // Read the boundary conditions. const int n_bnd_sw = rad_sw.get_n_bnd_gpu(); @@ -582,6 +581,12 @@ void solve_radiation(int argc, char** argv) Array mu0(input_nc.get_variable("mu0", {n_col_y, n_col_x}), {n_col}); Array azi(input_nc.get_variable("azi", {n_col_y, n_col_x}), {n_col}); + if (input_sza < 0) + mu0.fill(cos(input_sza / Float(180.0) * M_PI)); + + if (input_azi < 0) + azi.fill(input_azi / Float(180.0) * M_PI); + Array sfc_alb(input_nc.get_variable("sfc_alb_dir", {n_col_y, n_col_x, n_bnd_sw}), {n_bnd_sw, n_col}); Array tsi_scaling({n_col}); @@ -877,14 +882,14 @@ void solve_radiation(int argc, char** argv) nc_azi.insert(azi({1})/M_PI * Float(180.), {0}); // camera position and direction - Netcdf_group output_cam = output_nc.add_group("camera-settings"); - - std::string cam_vars[] = {"yaw","pitch","roll","px","py","pz"}; - for (auto &&cam_var : cam_vars) - { - auto nc_cam_out = output_cam.add_variable(cam_var); - nc_cam_out.insert(cam_in.get_variable(cam_var), {0}); - } + // Netcdf_group output_cam = output_nc.add_group("camera-settings"); + + // std::string cam_vars[] = {"yaw","pitch","roll","px","py","pz"}; + // for (auto &&cam_var : cam_vars) + // { + // auto nc_cam_out = output_cam.add_variable(cam_var); + // nc_cam_out.insert(cam_in.get_variable(cam_var), {0}); + // } } Status::print_message("###### Finished RTE+RRTMGP solver ######"); diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index 37d1b98e..4dd5911b 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -28,6 +28,8 @@ #include "raytracer_kernels_sw.h" #include "radiation_solver_rt.h" #include "aerosol_optics_rt.h" +#include "raytracer_kernels_bw.h" +#include "radiation_solver_bw.h" #include "gas_concs.h" #include "types.h" #include "mem_pool_gpu.h" @@ -168,23 +170,28 @@ void solve_radiation(int argc, char** argv) 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_cloud_optics = get_ini_value(settings, "switches", "cloud-optics", false); - bool switch_liq_cloud_optics = get_ini_value(settings, "switches", "liq-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_profiling = get_ini_value(settings, "switches", "profiling", 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); + 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; @@ -201,7 +208,7 @@ void solve_radiation(int argc, char** argv) 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, "ints", "sw_raytracing", 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 "; @@ -209,25 +216,26 @@ void solve_radiation(int argc, char** argv) } } + 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, "ints", "lw_raytracing", Int(22)); lw_photon_count = 1 << lw_photon_power; } - if (switch_cloud_optics) + if (switch_longwave && switch_bw_raytracing) { - switch_liq_cloud_optics = true; - switch_ice_cloud_optics = true; + Status::print_warning("No longwave radiation implemented in the backward ray tracer"); } - if (switch_liq_cloud_optics || switch_ice_cloud_optics) - switch_cloud_optics = true; if (switch_tica) switch_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!"); + if (switch_cloud_mie && switch_ice_cloud_optics) { std::string error = "Thou shall not use mie tables as long as ice optics are enabled"; @@ -236,7 +244,7 @@ void solve_radiation(int argc, char** argv) 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.); + 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"); @@ -245,6 +253,9 @@ 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."); @@ -269,20 +280,40 @@ void solve_radiation(int argc, char** argv) Array grid_y(input_nc.get_variable("y", {n_col_y}), {n_col_y}); Array grid_yh(input_nc.get_variable("yh", {n_col_y+1}), {n_col_y+1}); Array grid_z(input_nc.get_variable("z", {n_z_in}), {n_z_in}); + Array z_lev(input_nc.get_variable("z_lev", {n_lev}), {n_lev}); const Vector grid_d = {grid_xh({2}) - grid_xh({1}), grid_yh({2}) - grid_yh({1}), grid_z({2}) - grid_z({1})}; const Vector kn_grid = {input_nc.get_variable("ngrid_x"), 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}); Array t_lay(input_nc.get_variable("t_lay", {n_lay, n_col_y, n_col_x}), {n_col, n_lay}); 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) { std::string error = "col_dry is not supported in tica mode"; @@ -297,6 +328,19 @@ void solve_radiation(int argc, char** argv) } + // read land use map if present, used for choosing between spectral and lambertian reflection and for spectral albedo + // 0: water, 1: "grass", 2: "soil", 3: "concrete". Interpolating between 1 and 2 is currently possible + Array land_use_map({n_col}); + if (input_nc.variable_exists("land_use_map") && switch_lu_albedo) + { + land_use_map = std::move(input_nc.get_variable("land_use_map", {n_col_y, n_col_x})); + } + else + { + // default to grass with some soil + land_use_map.fill(Float(1.3)); + } + // Fetch the col_dry in case present. Array col_dry; if (input_nc.variable_exists("col_dry")) @@ -334,26 +378,24 @@ void solve_radiation(int argc, char** argv) Array rel; Array dei; - if (switch_cloud_optics) - { + const bool switch_cloud_optics = switch_liquid_cloud_optics || switch_ice_cloud_optics; - if (switch_liq_cloud_optics) - { - lwp.set_dims({n_col, n_lay}); - lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); + if (switch_liquid_cloud_optics) + { + lwp.set_dims({n_col, n_lay}); + lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); - rel.set_dims({n_col, n_lay}); - rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); - } + rel.set_dims({n_col, n_lay}); + rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); + } - if (switch_ice_cloud_optics) - { - iwp.set_dims({n_col, n_lay}); - iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); + if (switch_ice_cloud_optics) + { + iwp.set_dims({n_col, n_lay}); + iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); - dei.set_dims({n_col, n_lay}); - dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); - } + dei.set_dims({n_col, n_lay}); + dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); } Array rh; @@ -386,6 +428,12 @@ void solve_radiation(int argc, char** argv) mu0 = input_nc.get_variable("mu0", {n_col_y, n_col_x}); azi = input_nc.get_variable("azi", {n_col_y, n_col_x}); + if (input_sza < 0) + mu0.fill(cos(input_sza / Float(180.0) * M_PI)); + + if (input_azi < 0) + azi.fill(input_azi / Float(180.0) * M_PI); + if (switch_tica) { @@ -463,7 +511,7 @@ void solve_radiation(int argc, char** argv) lwp_out, iwp_out, rel_out, dei_out, rh_out, gas_concs_out, aerosol_concs_out, gas_names, aerosol_names, - switch_cloud_optics, switch_liq_cloud_optics, switch_ice_cloud_optics, switch_aerosol_optics + switch_liquid_cloud_optics, switch_ice_cloud_optics, switch_aerosol_optics ); lwp_out.expand_dims({n_col, n_lay}); @@ -498,6 +546,7 @@ void solve_radiation(int argc, char** argv) Netcdf_file output_nc(case_name + "_output.nc", Netcdf_mode::Create); output_nc.add_dimension("col", n_col); + output_nc.add_dimension("x", n_col_x); output_nc.add_dimension("y", n_col_y); output_nc.add_dimension("lay", n_lay); @@ -525,13 +574,27 @@ void solve_radiation(int argc, char** argv) Netcdf_file coef_nc_lw("coefficients_lw.nc", Netcdf_mode::Read); nbnds = std::max(coef_nc_lw.get_dimension_size("bnd"), nbnds); ngpts = std::max(coef_nc_lw.get_dimension_size("gpt"), ngpts); + + output_nc.add_dimension("gpt_lw", ngpts); + output_nc.add_dimension("band_lw", nbnds); + } - if (switch_shortwave) + if (switch_shortwave || switch_bw_raytracing) { Netcdf_file coef_nc_sw("coefficients_sw.nc", Netcdf_mode::Read); nbnds = std::max(coef_nc_sw.get_dimension_size("bnd"), nbnds); ngpts = std::max(coef_nc_sw.get_dimension_size("gpt"), ngpts); + + output_nc.add_dimension("gpt_sw", ngpts); + output_nc.add_dimension("band_sw", nbnds); + } + + if (switch_bw_raytracing) + { + output_nc.add_dimension("px", camera.nx); + output_nc.add_dimension("py", camera.ny); } + configure_memory_pool(n_lay, n_col, 1024, ngpts, nbnds); @@ -659,6 +722,8 @@ void solve_radiation(int argc, char** argv) switch_lw_raytracing, switch_cloud_optics, switch_aerosol_optics, + switch_delta_cloud, + switch_delta_aerosol, switch_single_gpt, switch_lw_scattering, switch_independent_column, @@ -699,14 +764,6 @@ void solve_radiation(int argc, char** argv) // Tuning step; run_solver(); - // Profiling step; - if (switch_profiling) - { - cudaProfilerStart(); - run_solver(); - cudaProfilerStop(); - } - //// Store the output. Status::print_message("Storing the longwave output."); Array lw_tau_tot_cpu(lw_tau_tot); @@ -741,9 +798,6 @@ void solve_radiation(int argc, char** argv) Array rt_flux_sfc_dn_cpu(rt_flux_sfc_dn); Array rt_flux_abs_cpu(rt_flux_abs); - output_nc.add_dimension("gpt_lw", n_gpt_lw); - output_nc.add_dimension("band_lw", n_bnd_lw); - 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}); @@ -835,7 +889,6 @@ void solve_radiation(int argc, char** argv) // Initialize the solver. Status::print_message("Initializing the shortwave solver."); - Gas_concs_gpu gas_concs_gpu(gas_concs); Radiation_solver_shortwave rad_sw(gas_concs_gpu, "coefficients_sw.nc", "cloud_coefficients_sw.nc", "aerosol_optics_sw.nc"); @@ -1037,14 +1090,6 @@ void solve_radiation(int argc, char** argv) // Tuning step; run_solver(); - // Profiling step; - if (switch_profiling) - { - cudaProfilerStart(); - run_solver(); - cudaProfilerStop(); - } - // Store the output. Status::print_message("Storing the shortwave output."); Array sw_tot_tau_cpu(sw_tot_tau); @@ -1072,9 +1117,6 @@ void solve_radiation(int argc, char** argv) Array rt_flux_abs_dir_cpu(rt_flux_abs_dir); Array rt_flux_abs_dif_cpu(rt_flux_abs_dif); - 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}); @@ -1217,6 +1259,297 @@ void solve_radiation(int argc, char** argv) } } + ////// RUN THE BACKWARD SOLVER ////// + if (switch_bw_raytracing || switch_cloud_cam) + { + // Initialize the solver. + Status::print_message("Initializing the shortwave backward solver."); + + Gas_concs_gpu gas_concs_gpu(gas_concs); + Radiation_solver_bw_shortwave rad_sw(gas_concs_gpu, "coefficients_sw.nc", "cloud_coefficients_sw.nc","aerosol_optics_sw.nc"); + + // Read the boundary conditions. + const int n_bnd_sw = rad_sw.get_n_bnd_gpu(); + const int n_gpt_sw = rad_sw.get_n_gpt_gpu(); + + Array mu0(input_nc.get_variable("mu0", {n_col_y, n_col_x}), {n_col}); + Array azi(input_nc.get_variable("azi", {n_col_y, n_col_x}), {n_col}); + + if (input_sza < 0) + mu0.fill(cos(input_sza / Float(180.0) * M_PI)); + + if (input_azi < 0) + azi.fill(input_azi / Float(180.0) * M_PI); + + Array sfc_alb(input_nc.get_variable("sfc_alb_dir", {n_col_y, n_col_x, n_bnd_sw}), {n_bnd_sw, n_col}); + + Array tsi_scaling({n_col}); + if (input_nc.variable_exists("tsi")) + { + Array tsi(input_nc.get_variable("tsi", {n_col_y, n_col_x}), {n_col}); + const Float tsi_ref = rad_sw.get_tsi_gpu(); + for (int icol=1; icol<=n_col; ++icol) + tsi_scaling({icol}) = tsi({icol}) / tsi_ref; + } + else if (input_nc.variable_exists("tsi_scaling")) + { + Float tsi_scaling_in = input_nc.get_variable("tsi_scaling"); + for (int icol=1; icol<=n_col; ++icol) + tsi_scaling({icol}) = tsi_scaling_in; + } + else + { + for (int icol=1; icol<=n_col; ++icol) + tsi_scaling({icol}) = Float(1.); + } + + Array_gpu XYZ; + Array_gpu radiance; + + if (switch_broadband) + { + radiance.set_dims({camera.nx, camera.ny}); + } + if (switch_image) + { + XYZ.set_dims({camera.nx, camera.ny, 3}); + } + + 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; + Array_gpu zen_cam; + + if (switch_cloud_cam) + { + liwp_cam.set_dims({camera.nx, camera.ny}); + tauc_cam.set_dims({camera.nx, camera.ny}); + dist_cam.set_dims({camera.nx, camera.ny}); + zen_cam.set_dims({camera.nx, camera.ny}); + } + + const Vector grid_cells = {n_col_x, n_col_y, n_z_in}; + + auto run_solver_bb = [&]() + { + Array_gpu p_lay_gpu(p_lay); + Array_gpu p_lev_gpu(p_lev); + Array_gpu t_lay_gpu(t_lay); + Array_gpu t_lev_gpu(t_lev); + Array_gpu z_lev_gpu(z_lev); + Array_gpu col_dry_gpu(col_dry); + Array_gpu sfc_alb_gpu(sfc_alb); + Array_gpu tsi_scaling_gpu(tsi_scaling); + Array_gpu mu0_gpu(mu0); + Array_gpu azi_gpu(azi); + Array_gpu lwp_gpu(lwp); + Array_gpu iwp_gpu(iwp); + Array_gpu rel_gpu(rel); + Array_gpu dei_gpu(dei); + Array_gpu rh_gpu(rh); + Aerosol_concs_gpu aerosol_concs_gpu(aerosol_concs); + + Array_gpu land_use_map_gpu(land_use_map); + + cudaDeviceSynchronize(); + cudaEvent_t start; + cudaEvent_t stop; + cudaEventCreate(&start); + cudaEventCreate(&stop); + + cudaEventRecord(start, 0); + + rad_sw.solve_gpu_bb( + switch_cloud_optics, + switch_cloud_mie, + switch_aerosol_optics, + switch_lu_albedo, + switch_delta_cloud, + switch_delta_aerosol, + switch_cloud_cam, + switch_bw_raytracing, + grid_cells, + grid_d, + kn_grid, + photons_per_pixel, + gas_concs_gpu, + p_lay_gpu, p_lev_gpu, + t_lay_gpu, t_lev_gpu, + z_lev_gpu, + col_dry_gpu, + sfc_alb_gpu, + tsi_scaling_gpu, + mu0_gpu, azi_gpu, + lwp_gpu, iwp_gpu, + rel_gpu, dei_gpu, + land_use_map_gpu, + rh_gpu, + aerosol_concs, + camera, + radiance, + liwp_cam, + tauc_cam, + dist_cam, + zen_cam); + + cudaEventRecord(stop, 0); + cudaEventSynchronize(stop); + float duration = 0.f; + cudaEventElapsedTime(&duration, start, stop); + + cudaEventDestroy(start); + cudaEventDestroy(stop); + + Status::print_message("Duration shortwave backward solver (broadband version): " + std::to_string(duration) + " (ms)"); + }; + + auto run_solver = [&]() + { + Array_gpu p_lay_gpu(p_lay); + Array_gpu p_lev_gpu(p_lev); + Array_gpu t_lay_gpu(t_lay); + Array_gpu t_lev_gpu(t_lev); + Array_gpu z_lev_gpu(z_lev); + Array_gpu col_dry_gpu(col_dry); + Array_gpu sfc_alb_gpu(sfc_alb); + Array_gpu tsi_scaling_gpu(tsi_scaling); + Array_gpu mu0_gpu(mu0); + Array_gpu azi_gpu(azi); + Array_gpu lwp_gpu(lwp); + Array_gpu iwp_gpu(iwp); + Array_gpu rel_gpu(rel); + Array_gpu dei_gpu(dei); + + Array_gpu rh_gpu(rh); + Aerosol_concs_gpu aerosol_concs_gpu(aerosol_concs); + + Array_gpu land_use_map_gpu(land_use_map); + + cudaDeviceSynchronize(); + cudaEvent_t start; + cudaEvent_t stop; + cudaEventCreate(&start); + cudaEventCreate(&stop); + + cudaEventRecord(start, 0); + + rad_sw.solve_gpu( + switch_cloud_optics, + switch_cloud_mie, + switch_aerosol_optics, + switch_lu_albedo, + switch_delta_cloud, + switch_delta_aerosol, + switch_cloud_cam, + switch_bw_raytracing, + grid_cells, + grid_d, + kn_grid, + photons_per_pixel, + gas_concs_gpu, + p_lay_gpu, p_lev_gpu, + t_lay_gpu, t_lev_gpu, + z_lev_gpu, + col_dry_gpu, + sfc_alb_gpu, + tsi_scaling_gpu, + mu0_gpu, azi_gpu, + lwp_gpu, iwp_gpu, + rel_gpu, dei_gpu, + land_use_map_gpu, + rh_gpu, + aerosol_concs, + camera, + XYZ, + liwp_cam, + tauc_cam, + dist_cam, + zen_cam); + + cudaEventRecord(stop, 0); + cudaEventSynchronize(stop); + float duration = 0.f; + cudaEventElapsedTime(&duration, start, stop); + + cudaEventDestroy(start); + cudaEventDestroy(stop); + + + Status::print_message("Duration shortwave backward solver (image version): " + std::to_string(duration) + " (ms)"); + }; + + if (switch_broadband) + { + run_solver_bb(); + } + if (switch_image) + { + run_solver(); + } + + // Store the output. + Status::print_message("Storing the backward output."); + + if (switch_bw_raytracing) + { + if (switch_broadband) + { + Array radiance_cpu(radiance); + + auto nc_var = output_nc.add_variable("radiance", {"py","px"}); + 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"); + } + if (switch_image) + { + Array xyz_cpu(XYZ); + output_nc.add_dimension("n",3); + + auto nc_xyz = output_nc.add_variable("XYZ", {"n","py","px"}); + nc_xyz.insert(xyz_cpu.v(), {0, 0, 0}); + + nc_xyz.add_attribute("long_name", "X Y Z tristimulus values"); + } + } + + if (switch_cloud_cam) + { + + Array liwp_cam_cpu(liwp_cam); + Array tauc_cam_cpu(tauc_cam); + Array dist_cam_cpu(dist_cam); + Array zen_cam_cpu(zen_cam); + + auto nc_var_liwp = output_nc.add_variable("liq_ice_wp_cam", {"py","px"}); + nc_var_liwp.insert(liwp_cam_cpu.v(), {0, 0}); + nc_var_liwp.add_attribute("long_name", "accumulated liquid+ice water path"); + + auto nc_var_tauc = output_nc.add_variable("tau_cld_cam", {"py","px"}); + nc_var_tauc.insert(tauc_cam_cpu.v(), {0, 0}); + nc_var_tauc.add_attribute("long_name", "accumulated cloud optical depth (441-615nm band)"); + + auto nc_var_dist = output_nc.add_variable("dist_cld_cam", {"py","px"}); + nc_var_dist.insert(dist_cam_cpu.v(), {0, 0}); + nc_var_dist.add_attribute("long_name", "distance to first cloudy cell"); + + auto nc_var_csza = output_nc.add_variable("zen_cam", {"py","px"}); + nc_var_csza.insert(zen_cam_cpu.v(), {0, 0}); + nc_var_csza.add_attribute("long_name", "zenith angle of camera pixel"); + } + + auto nc_mu0 = output_nc.add_variable("sza"); + nc_mu0.insert(acos(mu0({1}))/M_PI * Float(180.), {0}); + + auto nc_azi = output_nc.add_variable("azi"); + nc_azi.insert(azi({1})/M_PI * Float(180.), {0}); + + + } + Status::print_message("###### Finished RTE+RRTMGP solver ######"); } diff --git a/src_tilt/tilt_utils.cpp b/src_tilt/tilt_utils.cpp index fa78aacb..116349e4 100644 --- a/src_tilt/tilt_utils.cpp +++ b/src_tilt/tilt_utils.cpp @@ -227,12 +227,12 @@ void restore_bkg_profile_bundle(const int n_col_x, const int n_col_y, Array* lwp, Array* iwp, Array* rel, Array* dei, Array* rh, Gas_concs& gas_concs, Aerosol_concs& aerosol_concs, std::vector gas_names, std::vector aerosol_names, - bool switch_liq_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics + bool switch_liquid_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics ) { const int n_col = n_col_x*n_col_y; - if (switch_liq_cloud_optics) + if (switch_liquid_cloud_optics) { restore_bkg_profile(n_col_x, n_col_y, n_lay, n_z_in, bkg_start_z, lwp_copy->v(), lwp->v()); lwp_copy->expand_dims({n_col, n_lay_tot}); @@ -333,7 +333,7 @@ void post_process_output(const std::vector& col_results, Array* iwp_out, Array* rel_out, Array* dei_out, - const bool switch_liq_cloud_optics, + const bool switch_liquid_cloud_optics, const bool switch_ice_cloud_optics) { const int total_cols = n_col_x * n_col_y; @@ -343,7 +343,7 @@ void post_process_output(const std::vector& col_results, const ColumnResult& col = col_results[col_idx]; for (int j = 0; j < n_z; ++j) { int out_idx = col_idx + j * total_cols; - if (switch_liq_cloud_optics) { + if (switch_liquid_cloud_optics) { lwp_out->v()[out_idx] = col.lwp.v()[j]; rel_out->v()[out_idx] = col.rel.v()[j]; } @@ -822,7 +822,7 @@ void tica_tilt( Array& lwp_out, Array& iwp_out, Array& rel_out, Array& dei_out, Array& rh_out, Gas_concs& gas_concs_out, Aerosol_concs& aerosol_concs_out, std::vector gas_names, std::vector aerosol_names, - bool switch_cloud_optics, bool switch_liq_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics + bool switch_liquid_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics ) { // if t lev all 0, interpolate from t lay @@ -869,7 +869,7 @@ void tica_tilt( gas_concs_out, gas_names, aerosol_concs_out, aerosol_names, switch_aerosol_optics); ////// SETUP FOR RANDOM START POINT TILTING ////// - if (switch_cloud_optics) + if (switch_liquid_cloud_optics || switch_ice_cloud_optics) { std::vector> x_y_start_arr(n_col); std::vector> by_col_paths(n_col); @@ -954,7 +954,7 @@ void tica_tilt( for (int ilay = 1; ilay <= n_zh_in; ++ilay) { Float dz = zh({ilay + 1}) - zh({ilay}); - if (switch_liq_cloud_optics) + if (switch_liquid_cloud_optics) { (lwp_norm_reshaped)({ilay, icol}) = (lwp)({icol, ilay})/dz; (rel_reshaped)({ilay, icol}) = (rel)({icol, ilay}); @@ -967,7 +967,7 @@ void tica_tilt( } } - if (switch_liq_cloud_optics){ + if (switch_liquid_cloud_optics){ Array lwp_compress; lwp_compress.set_dims({n_z_in}); Array rel_compress; @@ -1151,7 +1151,7 @@ void tica_tilt( post_process_output(col_results, n_col_x, n_col_y, n_z_in, n_zh_in, &lwp_out, &iwp_out, &rel_out, &dei_out, - switch_liq_cloud_optics, switch_ice_cloud_optics); + switch_liquid_cloud_optics, switch_ice_cloud_optics); } @@ -1165,7 +1165,7 @@ void tica_tilt( &lwp, &iwp, &rel, &dei, &rh, gas_concs, aerosol_concs, gas_names, aerosol_names, - switch_liq_cloud_optics, switch_ice_cloud_optics, switch_aerosol_optics + switch_liquid_cloud_optics, switch_ice_cloud_optics, switch_aerosol_optics ); }