diff --git a/allsky/allsky.ini b/allsky/allsky.ini index 8e1fea5..5277cbd 100644 --- a/allsky/allsky.ini +++ b/allsky/allsky.ini @@ -1,32 +1,29 @@ -[switches] -shortwave = true -longwave = true -fluxes = true +[clouds] liquid_cloud_optics = true ice_cloud_optics = true -aerosol_optics = false delta_cloud = true +tilted_columns = false +cloud_mie = false + +[aerosols] +aerosol_optics = false delta_aerosol = false -# gpu_only (test_rte_rrtmgp_gpu) -timings = false +[shortwave] +shortwave = true + +# compute ray tracer fluxes using "samples" photons per pixel per spectral quadrature point +raytracing = false + +# compute plane parallel (two-stream) fluxes +plane_parallel = true + +[longwave] +longwave = true + +# compute ray tracer fluxes using 2**samples photons per spectral quadrature point +raytracing = false + +# compute plane parallel fluxes (two-stream or no-scat solution depending on "scattering") +plane_parallel = true -# raytracer (test_rte_rrtmgp_rt_gpu) -sw_two_stream = false -sw_raytracing = true -lw_raytracing = true -independent_column = false -lw_scattering = false -cloud_mie = false -single_gpt = false -profiling = false -min_mfp_grid_ratio = true -tica = false - -[ints] -sw_raytracing = 256 -lw_raytracing = 22 -single_gpt = 1 - -[floats] -min_mfp_grid_ratio = 1.0 diff --git a/include_test/radiation_solver_rt.h b/include_test/radiation_solver_rt.h index 9df1af2..d26029c 100644 --- a/include_test/radiation_solver_rt.h +++ b/include_test/radiation_solver_rt.h @@ -46,13 +46,12 @@ class Radiation_solver_longwave #ifdef __CUDACC__ void solve_gpu( - const bool switch_fluxes, const bool switch_raytracing, + const bool switch_plane_parallel, const bool switch_cloud_optics, const bool switch_aerosol_optics, const bool switch_delta_cloud, const bool switch_delta_aerosol, - const bool switch_single_gpt, const bool switch_lw_scattering, const bool switch_independent_column, const int single_gpt, @@ -75,7 +74,6 @@ class Radiation_solver_longwave Array_gpu& aer_tau_out,Array_gpu& aer_ssa_out, Array_gpu& aer_asy_out, Array_gpu& lay_source, Array_gpu& lev_source, Array_gpu& sfc_source, Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net, - Array_gpu& lw_gpt_flux_up, Array_gpu& lw_gpt_flux_dn, Array_gpu& lw_gpt_flux_net, Array_gpu& rt_flux_tod_up, Array_gpu& rt_flux_tod_dn, Array_gpu& rt_flux_sfc_up, Array_gpu& rt_flux_sfc_dn, Array_gpu& rt_flux_abs); @@ -122,14 +120,12 @@ class Radiation_solver_shortwave #ifdef __CUDACC__ void solve_gpu( - const bool switch_fluxes, - const bool switch_twostream, const bool switch_raytracing, + const bool switch_plane_parallel, const bool switch_independent_column, const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, - const bool switch_single_gpt, const bool switch_delta_cloud, const bool switch_delta_aerosol, const bool switch_attenuate_tica, @@ -154,8 +150,6 @@ class Radiation_solver_shortwave Array_gpu& aer_tau_out, Array_gpu& aer_ssa_out, Array_gpu& aer_asy_out, Array_gpu& sw_flux_up, Array_gpu& sw_flux_dn, Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net, - Array_gpu& sw_gpt_flux_up, Array_gpu& sw_gpt_flux_dn, - Array_gpu& sw_gpt_flux_dn_dir, Array_gpu& sw_gpt_flux_net, Array_gpu& rt_flux_tod_up, Array_gpu& rt_flux_sfc_dir, Array_gpu& rt_flux_sfc_dif, diff --git a/les_cloudfield/test.ini b/les_cloudfield/test.ini index f35ae36..40d58df 100644 --- a/les_cloudfield/test.ini +++ b/les_cloudfield/test.ini @@ -1,77 +1,94 @@ -[switches] -##### General switches +[clouds] +liquid_cloud_optics = true +ice_cloud_optics = true +delta_cloud = false +tilted_columns = false +cloud_mie = false + +[aerosols] +aerosol_optics = false +delta_aerosol = false + +[shortwave] shortwave = true -longwave = true -fluxes = true -cloud-optics = true -delta-cloud = false +# compute ray tracer fluxes using "samples" photons per pixel per spectral quadrature point +raytracing = true +samples = 256 -aerosol-optics = true -delta-aerosol = false +# ray tracing in independent column mode +rt_independent_column = false -##### test_rte_rrtmgp only (CPU & GPU) -output-optical = false -output-bnd-fluxes = false +# compute plane parallel (two-stream) fluxes +plane_parallel = false -##### test_rte_rrtmgp only (GPU) -timings = false +# only compute specified g-point +single_gpt = -1 -##### Forward raytracer only +[longwave] +longwave = true -# switch on ray tracer. Specifiy ray counts under [ints]: samples per pixel for sw, for lw provide n to use 2^n samples -sw-raytracing = true -lw-raytracing = true +# compute ray tracer fluxes using 2**samples photons per spectral quadrature point +raytracing = true +samples = 22 -# Also perform two-stream computations for shortwave. -sw-two-stream = false +# ray tracing in independent column mode +rt_independent_column = false -# Independent column mode -independent-column = false -# Tilted independent column mode -tica = false +# enable scattering +scattering = false -# Enable scattering for longwave. If enabled, 1D solver will also switch from no-scat solution to two-stream formulation -lw-scattering = true +# compute plane parallel fluxes (two-stream or no-scat solution depending on "scattering") +plane_parallel = false -# Only perform ray tracing when the minimum gasous mean free path is larger than (min-mfp-grid-ratio * grid spacing) and use a 1D solutions otherwise. -# Specify threshold ratio under [floats] -min-mfp-grid-ratio = true +# only compute specified g-point +single_gpt = -1 -# liq and ice flags are automatically switch on when cloud-optics is set to true -# to run either liquid or ice cloud optics -liq-cloud-optics = false -ice-cloud-optics = false +# threshold ratio between lowest gaseus mean free path and horizontal grid spacing +# if this ratio falls below the thershold for a g-point, plane parallel fluxes are used for that g-point instead of the raytracer +min_mf_grid_ratio = 0.0 -##### Common raytracer (forward & backward) switches -cloud-mie = false +[backward] +# enable backward ray tracing using "samples" photons per pixel" +bw_raytracing = true +samples = 16 -# Save optical properties and fluxes for a single g-point. Specify under [ints] -single-gpt = false +# raytrace bands that fall within visible spectrum and convert radiances to XYZ tristimulus values +image = true -# Enable profiling -profiling = false +# compute broadband radiances +broadband = true -##### Backward raytracer only -raytracing = true +# output additional cloud information for each camera pixel +cloud_cam = true -# Read variable "surface_type" to determine scattering type and spectral albedo for each surface grid cell. -# Set either to 0 (water; spectrally constant albedo; specular reflection) or a number between 1 (soil-like spectral albedo; lambertian) and 2 (grass-line spectral albedo; lambertian). -lu-albedo = true +# read surface type from "land_use_map" input variable and use corresponding spectral albedo and scattering type (lambertian, specular) +lu_albedo = true -# broadband radiance computation (sum over all g-points) -broadband = true +[solar_angles] +# solar zenith angle in degrees. mu0 from input file is used if sza < 0 +sza = 45 -# image computation: run only g-points in visible spectrum and convert to XYZ tristimulus values -image = true +# solar azimuth angle in degrees. azi from input file is used if azi < 0 +azi = 135 + +[camera] +cam_field_of_view = 50.0 + +# camera type: (0) fish eye camera, (1) rectangular camera, (2) top-of-atmosphere upwelling radiances +cam_type = 1 -[ints] -single-gpt = 1 +# x,y,z positions of virtual camera +cam_px = 0.0 +cam_py = 0.0 +cam_pz = 2500.0 -sw-raytracing = 256 -lw-raytracing = 22 +# width, height (pixels) of virtual camera +cam_nx = 512 +cam_ny = 512 -raytracing = 128 +# yaw, pitch and roll angles (degrees) of the virtual camera +cam_yaw = 135.0 +cam_pitch = 55.0 +cam_roll = 0.0 -[floats] -min-mfp-grid-ratio = 1.0 diff --git a/rcemip/rcemip.ini b/rcemip/rcemip.ini index c81c2e4..0229057 100644 --- a/rcemip/rcemip.ini +++ b/rcemip/rcemip.ini @@ -1,32 +1,29 @@ -[switches] -shortwave = true -longwave = true -fluxes = true +[clouds] liquid_cloud_optics = false ice_cloud_optics = false +delta_cloud = false +tilted_columns = false +cloud_mie = false + +[aerosols] aerosol_optics = false -delta_cloud = true delta_aerosol = false -# gpu_only (test_rte_rrtmgp_gpu) -timings = false +[shortwave] +shortwave = true + +# compute ray tracer fluxes using "samples" photons per pixel per spectral quadrature point +raytracing = false + +# compute plane parallel (two-stream) fluxes +plane_parallel = true + +[longwave] +longwave = true + +# compute ray tracer fluxes using 2**samples photons per spectral quadrature point +raytracing = false + +# compute plane parallel fluxes (two-stream or no-scat solution depending on "scattering") +plane_parallel = true -# raytracer (test_rte_rrtmgp_rt_gpu) -sw_two_stream = false -sw_raytracing = true -lw_raytracing = true -independent_column = false -lw_scattering = false -cloud_mie = false -single_gpt = false -profiling = false -min_mfp_grid_ratio = true -tica = false - -[ints] -sw_raytracing = 256 -lw_raytracing = 22 -single_gpt = 1 - -[floats] -min_mfp_grid_ratio = 1.0 diff --git a/rfmip/rfmip.ini b/rfmip/rfmip.ini index c81c2e4..0229057 100644 --- a/rfmip/rfmip.ini +++ b/rfmip/rfmip.ini @@ -1,32 +1,29 @@ -[switches] -shortwave = true -longwave = true -fluxes = true +[clouds] liquid_cloud_optics = false ice_cloud_optics = false +delta_cloud = false +tilted_columns = false +cloud_mie = false + +[aerosols] aerosol_optics = false -delta_cloud = true delta_aerosol = false -# gpu_only (test_rte_rrtmgp_gpu) -timings = false +[shortwave] +shortwave = true + +# compute ray tracer fluxes using "samples" photons per pixel per spectral quadrature point +raytracing = false + +# compute plane parallel (two-stream) fluxes +plane_parallel = true + +[longwave] +longwave = true + +# compute ray tracer fluxes using 2**samples photons per spectral quadrature point +raytracing = false + +# compute plane parallel fluxes (two-stream or no-scat solution depending on "scattering") +plane_parallel = true -# raytracer (test_rte_rrtmgp_rt_gpu) -sw_two_stream = false -sw_raytracing = true -lw_raytracing = true -independent_column = false -lw_scattering = false -cloud_mie = false -single_gpt = false -profiling = false -min_mfp_grid_ratio = true -tica = false - -[ints] -sw_raytracing = 256 -lw_raytracing = 22 -single_gpt = 1 - -[floats] -min_mfp_grid_ratio = 1.0 diff --git a/src_test/radiation_solver_rt.cu b/src_test/radiation_solver_rt.cu index 7278181..fa9e9ec 100644 --- a/src_test/radiation_solver_rt.cu +++ b/src_test/radiation_solver_rt.cu @@ -479,13 +479,12 @@ Radiation_solver_longwave::Radiation_solver_longwave( void Radiation_solver_longwave::solve_gpu( - const bool switch_fluxes, const bool switch_raytracing, + const bool switch_plane_parallel, const bool switch_cloud_optics, const bool switch_aerosol_optics, const bool switch_delta_cloud, const bool switch_delta_aerosol, - const bool switch_single_gpt, const bool switch_lw_scattering, const bool switch_independent_column, const int single_gpt, @@ -508,7 +507,6 @@ void Radiation_solver_longwave::solve_gpu( Array_gpu& aer_tau_out, Array_gpu& aer_ssa_out, Array_gpu& aer_asy_out, Array_gpu& lay_source, Array_gpu& lev_source, Array_gpu& sfc_source, Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net, - Array_gpu& lw_gpt_flux_up, Array_gpu& lw_gpt_flux_dn, Array_gpu& lw_gpt_flux_net, Array_gpu& rt_flux_tod_up, Array_gpu& rt_flux_tod_dn, Array_gpu& rt_flux_sfc_up, Array_gpu& rt_flux_sfc_dn, Array_gpu& rt_flux_abs) { @@ -533,15 +531,11 @@ void Radiation_solver_longwave::solve_gpu( Gas_optics_rrtmgp_rt::get_col_dry(col_dry, gas_concs.get_vmr("h2o"), p_lev); } - if (switch_fluxes) + if (switch_plane_parallel) { Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, n_col, lw_flux_up.ptr()); Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, n_col, lw_flux_dn.ptr()); Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, n_col, lw_flux_net.ptr()); - - if (switch_raytracing) { } - - } const Array& band_limits_gpt(this->kdist_gpu->get_band_lims_gpoint()); @@ -551,6 +545,8 @@ void Radiation_solver_longwave::solve_gpu( for (int igpt=1; igpt<=n_gpt; ++igpt) { + if ((single_gpt > 0) && (igpt != single_gpt)) continue; + int band = 0; for (int ibnd=1; ibnd<=n_bnd; ++ibnd) { @@ -606,7 +602,6 @@ void Radiation_solver_longwave::solve_gpu( gas_optics_subset(col_s, n_col_residual); } - // Find maximum gasous optical depth to and compute lowest mean free path on the clearsky atmosphere // Allocate temporary storage @@ -691,8 +686,8 @@ void Radiation_solver_longwave::solve_gpu( Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_g().ptr()); } - // Store the optical properties, if desired. - if (switch_single_gpt && igpt == single_gpt) + // Store the optical properties, only in single_gpt mode + if (single_gpt > 0) { lay_source = (*sources).get_lay_source(); lev_source = (*sources).get_lev_source(); @@ -713,13 +708,15 @@ void Radiation_solver_longwave::solve_gpu( } - if (switch_fluxes) - { - constexpr int n_ang = 1; + constexpr int n_ang = 1; + + const bool use_raytracer = switch_raytracing && (lowest_gas_mean_free_path / grid_d_xy_min) > min_mfp_grid_ratio; - std::unique_ptr fluxes = - std::make_unique(grid_cells.x, grid_cells.y, grid_cells.z, n_lev); + std::unique_ptr fluxes = + std::make_unique(grid_cells.x, grid_cells.y, grid_cells.z, n_lev); + if (switch_plane_parallel || !use_raytracer) + { rte_lw.rte_lw( optical_props, top_at_1, @@ -733,78 +730,73 @@ void Radiation_solver_longwave::solve_gpu( (*fluxes).net_flux(); // Copy the data to the output. - Gpt_combine_kernels_cuda_rt::add_from_gpoint( + if (switch_plane_parallel) + { + Gpt_combine_kernels_cuda_rt::add_from_gpoint( n_col, n_lev, lw_flux_up.ptr(), lw_flux_dn.ptr(), lw_flux_net.ptr(), (*fluxes).get_flux_up().ptr(), (*fluxes).get_flux_dn().ptr(), (*fluxes).get_flux_net().ptr()); - - if (switch_single_gpt && igpt == single_gpt) - { - lw_gpt_flux_up = (*fluxes).get_flux_up(); - lw_gpt_flux_dn = (*fluxes).get_flux_dn(); - lw_gpt_flux_net = (*fluxes).get_flux_net(); } + } - if (switch_raytracing) + if (switch_raytracing) + { + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_tod_dn().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_tod_up().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_sfc_dif().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_sfc_up().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, grid_cells.z, (*fluxes).get_flux_abs_dif().ptr()); + + if (use_raytracer) { - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_tod_dn().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_tod_up().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_sfc_dif().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, (*fluxes).get_flux_sfc_up().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.x, grid_cells.y, grid_cells.z, (*fluxes).get_flux_abs_dif().ptr()); + ++monte_carlo_gpoints; - if ( (lowest_gas_mean_free_path / grid_d_xy_min) > min_mfp_grid_ratio) - { - ++monte_carlo_gpoints; - - raytracer_lw.trace_rays( - igpt, - switch_independent_column, - ray_count, - grid_cells, - grid_d, - kn_grid, - dynamic_cast(*optical_props).get_tau(), - dynamic_cast(*optical_props).get_ssa(), - dynamic_cast(*cloud_optical_props).get_tau(), - dynamic_cast(*cloud_optical_props).get_ssa(), - dynamic_cast(*cloud_optical_props).get_g(), - dynamic_cast(*aerosol_optical_props).get_tau(), - dynamic_cast(*aerosol_optical_props).get_ssa(), - dynamic_cast(*aerosol_optical_props).get_g(), - (*sources).get_lay_source(), - (*sources).get_sfc_source(), - emis_sfc, - (*fluxes).get_flux_dn()({1, grid_cells.z}), - (*fluxes).get_flux_tod_dn(), - (*fluxes).get_flux_tod_up(), - (*fluxes).get_flux_sfc_dif(), - (*fluxes).get_flux_sfc_up(), - (*fluxes).get_flux_abs_dif()); - } - else - { - convert_1d_to_rt_output( - n_col, n_lay, grid_cells.z, grid_d.z, - (*fluxes).get_flux_up(), - (*fluxes).get_flux_dn(), - (*fluxes).get_flux_net(), + raytracer_lw.trace_rays( + igpt, + switch_independent_column, + ray_count, + grid_cells, + grid_d, + kn_grid, + dynamic_cast(*optical_props).get_tau(), + dynamic_cast(*optical_props).get_ssa(), + dynamic_cast(*cloud_optical_props).get_tau(), + dynamic_cast(*cloud_optical_props).get_ssa(), + dynamic_cast(*cloud_optical_props).get_g(), + dynamic_cast(*aerosol_optical_props).get_tau(), + dynamic_cast(*aerosol_optical_props).get_ssa(), + dynamic_cast(*aerosol_optical_props).get_g(), + (*sources).get_lay_source(), + (*sources).get_sfc_source(), + emis_sfc, + (*fluxes).get_flux_dn()({1, grid_cells.z}), (*fluxes).get_flux_tod_dn(), (*fluxes).get_flux_tod_up(), (*fluxes).get_flux_sfc_dif(), (*fluxes).get_flux_sfc_up(), (*fluxes).get_flux_abs_dif()); + } + else + { + convert_1d_to_rt_output( + n_col, n_lay, grid_cells.z, grid_d.z, + (*fluxes).get_flux_up(), + (*fluxes).get_flux_dn(), + (*fluxes).get_flux_net(), + (*fluxes).get_flux_tod_dn(), + (*fluxes).get_flux_tod_up(), + (*fluxes).get_flux_sfc_dif(), + (*fluxes).get_flux_sfc_up(), + (*fluxes).get_flux_abs_dif()); - } - - Gpt_combine_kernels_cuda_rt::add_from_gpoint( - grid_cells.x, grid_cells.y, - rt_flux_tod_dn.ptr(), rt_flux_tod_up.ptr(), rt_flux_sfc_dn.ptr(), rt_flux_sfc_up.ptr(), - (*fluxes).get_flux_tod_dn().ptr(), (*fluxes).get_flux_tod_up().ptr(), (*fluxes).get_flux_sfc_dif().ptr(), (*fluxes).get_flux_sfc_up().ptr()); + } - Gpt_combine_kernels_cuda_rt::add_from_gpoint( - n_col, grid_cells.z, rt_flux_abs.ptr(), (*fluxes).get_flux_abs_dif().ptr()); + Gpt_combine_kernels_cuda_rt::add_from_gpoint( + grid_cells.x, grid_cells.y, + rt_flux_tod_dn.ptr(), rt_flux_tod_up.ptr(), rt_flux_sfc_dn.ptr(), rt_flux_sfc_up.ptr(), + (*fluxes).get_flux_tod_dn().ptr(), (*fluxes).get_flux_tod_up().ptr(), (*fluxes).get_flux_sfc_dif().ptr(), (*fluxes).get_flux_sfc_up().ptr()); - } + Gpt_combine_kernels_cuda_rt::add_from_gpoint( + n_col, grid_cells.z, rt_flux_abs.ptr(), (*fluxes).get_flux_abs_dif().ptr()); } Tools_gpu::free_gpu(max_tau_gas_g); @@ -846,14 +838,12 @@ void Radiation_solver_shortwave::load_mie_tables( } void Radiation_solver_shortwave::solve_gpu( - const bool switch_fluxes, - const bool switch_twostream, const bool switch_raytracing, + const bool switch_plane_parallel, const bool switch_independent_column, const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, - const bool switch_single_gpt, const bool switch_delta_cloud, const bool switch_delta_aerosol, const bool switch_attenuate_tica, @@ -878,8 +868,6 @@ void Radiation_solver_shortwave::solve_gpu( Array_gpu& aer_tau_out, Array_gpu& aer_ssa_out, Array_gpu& aer_asy_out, Array_gpu& sw_flux_up, Array_gpu& sw_flux_dn, Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net, - Array_gpu& sw_gpt_flux_up, Array_gpu& sw_gpt_flux_dn, - Array_gpu& sw_gpt_flux_dn_dir, Array_gpu& sw_gpt_flux_net, Array_gpu& rt_flux_tod_up, Array_gpu& rt_flux_sfc_dir, Array_gpu& rt_flux_sfc_dif, @@ -914,24 +902,22 @@ void Radiation_solver_shortwave::solve_gpu( Array_gpu mie_cdfs_sub; Array_gpu mie_angs_sub; - if (switch_fluxes) + if (switch_plane_parallel) { - if (switch_twostream) - { - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_up.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn_dir.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_net.ptr()); - } - if (switch_raytracing) - { - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_tod_up.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_dir.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_dif.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_up.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.z, grid_cells.y, grid_cells.x, rt_flux_abs_dir.ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.z, grid_cells.y, grid_cells.x, rt_flux_abs_dif.ptr()); - } + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_up.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_dn_dir.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_lev, grid_cells.y, grid_cells.x, sw_flux_net.ptr()); + } + + if (switch_raytracing) + { + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_tod_up.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_dir.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_dif.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.y, grid_cells.x, rt_flux_sfc_up.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.z, grid_cells.y, grid_cells.x, rt_flux_abs_dir.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(grid_cells.z, grid_cells.y, grid_cells.x, rt_flux_abs_dif.ptr()); } const Array& band_limits_gpt(this->kdist_gpu->get_band_lims_gpoint()); @@ -940,6 +926,8 @@ void Radiation_solver_shortwave::solve_gpu( for (int igpt=1; igpt<=n_gpt; ++igpt) { + if ((single_gpt > 0) && (igpt != single_gpt)) continue; + int band = 0; for (int ibnd=1; ibnd<=n_bnd; ++ibnd) { @@ -1063,8 +1051,8 @@ void Radiation_solver_shortwave::solve_gpu( Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_g().ptr()); } - // Store the optical properties, if desired - if (switch_single_gpt && igpt == single_gpt) + // Store the optical properties, only in single_gpt mode + if (single_gpt > 0) { tot_tau_out = optical_props->get_tau(); tot_ssa_out = optical_props->get_ssa(); @@ -1076,97 +1064,86 @@ void Radiation_solver_shortwave::solve_gpu( aer_asy_out = aerosol_optical_props->get_g(); } - if (switch_fluxes) - { - std::unique_ptr fluxes = - std::make_unique(grid_cells.x, grid_cells.y, grid_cells.z, n_lev); + std::unique_ptr fluxes = + std::make_unique(grid_cells.x, grid_cells.y, grid_cells.z, n_lev); - if (switch_twostream) - { + if (switch_plane_parallel) + { + rte_sw.rte_sw( + optical_props, + top_at_1, + mu0, + toa_src, + sfc_alb_dir.subset({{ {band, band}, {1, n_col}}}), + sfc_alb_dif.subset({{ {band, band}, {1, n_col}}}), + Array_gpu(), // Add an empty array, no inc_flux. + (*fluxes).get_flux_up(), + (*fluxes).get_flux_dn(), + (*fluxes).get_flux_dn_dir()); + } - rte_sw.rte_sw( - optical_props, - top_at_1, - mu0, - toa_src, - sfc_alb_dir.subset({{ {band, band}, {1, n_col}}}), - sfc_alb_dif.subset({{ {band, band}, {1, n_col}}}), - Array_gpu(), // Add an empty array, no inc_flux. - (*fluxes).get_flux_up(), - (*fluxes).get_flux_dn(), - (*fluxes).get_flux_dn_dir()); - } + if (switch_raytracing) + { + Float zenith_angle = std::acos(mu0({1})); + Float azimuth_angle = azi({1}); // sun approximately from south - if (switch_raytracing) + if (switch_cloud_mie) { - Float zenith_angle = std::acos(mu0({1})); - Float azimuth_angle = azi({1}); // sun approximately from south - - if (switch_cloud_mie) - { - mie_cdfs_sub = mie_cdfs.subset({{ {1, n_mie}, {band, band} }}); - mie_angs_sub = mie_angs.subset({{ {1, n_mie}, {1, n_re}, {band, band} }}); - } - - raytracer.trace_rays( - igpt, - switch_independent_column, - ray_count, - grid_cells, - grid_d, - kn_grid, - mie_cdfs_sub, - mie_angs_sub, - dynamic_cast(*optical_props).get_tau(), - dynamic_cast(*optical_props).get_ssa(), - dynamic_cast(*cloud_optical_props).get_tau(), - dynamic_cast(*cloud_optical_props).get_ssa(), - dynamic_cast(*cloud_optical_props).get_g(), - dynamic_cast(*aerosol_optical_props).get_tau(), - dynamic_cast(*aerosol_optical_props).get_ssa(), - dynamic_cast(*aerosol_optical_props).get_g(), - rel, - sfc_alb_dir.subset({{ {band, band}, {1, n_col}}}), - zenith_angle, - azimuth_angle, - toa_src({1}) * mu0({1}), - Float(0.), - (*fluxes).get_flux_tod_dn(), - (*fluxes).get_flux_tod_up(), - (*fluxes).get_flux_sfc_dir(), - (*fluxes).get_flux_sfc_dif(), - (*fluxes).get_flux_sfc_up(), - (*fluxes).get_flux_abs_dir(), - (*fluxes).get_flux_abs_dif()); + mie_cdfs_sub = mie_cdfs.subset({{ {1, n_mie}, {band, band} }}); + mie_angs_sub = mie_angs.subset({{ {1, n_mie}, {1, n_re}, {band, band} }}); } - (*fluxes).net_flux(); - if (switch_twostream) - { - Gpt_combine_kernels_cuda_rt::add_from_gpoint( - n_col, n_lev, sw_flux_up.ptr(), sw_flux_dn.ptr(), sw_flux_dn_dir.ptr(), sw_flux_net.ptr(), - (*fluxes).get_flux_up().ptr(), (*fluxes).get_flux_dn().ptr(), (*fluxes).get_flux_dn_dir().ptr(), (*fluxes).get_flux_net().ptr()); - } + raytracer.trace_rays( + igpt, + switch_independent_column, + ray_count, + grid_cells, + grid_d, + kn_grid, + mie_cdfs_sub, + mie_angs_sub, + dynamic_cast(*optical_props).get_tau(), + dynamic_cast(*optical_props).get_ssa(), + dynamic_cast(*cloud_optical_props).get_tau(), + dynamic_cast(*cloud_optical_props).get_ssa(), + dynamic_cast(*cloud_optical_props).get_g(), + dynamic_cast(*aerosol_optical_props).get_tau(), + dynamic_cast(*aerosol_optical_props).get_ssa(), + dynamic_cast(*aerosol_optical_props).get_g(), + rel, + sfc_alb_dir.subset({{ {band, band}, {1, n_col}}}), + zenith_angle, + azimuth_angle, + toa_src({1}) * mu0({1}), + Float(0.), + (*fluxes).get_flux_tod_dn(), + (*fluxes).get_flux_tod_up(), + (*fluxes).get_flux_sfc_dir(), + (*fluxes).get_flux_sfc_dif(), + (*fluxes).get_flux_sfc_up(), + (*fluxes).get_flux_abs_dir(), + (*fluxes).get_flux_abs_dif()); + } - if (switch_raytracing) - { - Gpt_combine_kernels_cuda_rt::add_from_gpoint( - grid_cells.x, grid_cells.y, rt_flux_tod_up.ptr(), rt_flux_sfc_dir.ptr(), rt_flux_sfc_dif.ptr(), rt_flux_sfc_up.ptr(), - (*fluxes).get_flux_tod_up().ptr(), (*fluxes).get_flux_sfc_dir().ptr(), (*fluxes).get_flux_sfc_dif().ptr(), (*fluxes).get_flux_sfc_up().ptr()); + (*fluxes).net_flux(); + if (switch_plane_parallel) + { + Gpt_combine_kernels_cuda_rt::add_from_gpoint( + n_col, n_lev, sw_flux_up.ptr(), sw_flux_dn.ptr(), sw_flux_dn_dir.ptr(), sw_flux_net.ptr(), + (*fluxes).get_flux_up().ptr(), (*fluxes).get_flux_dn().ptr(), (*fluxes).get_flux_dn_dir().ptr(), (*fluxes).get_flux_net().ptr()); + } - Gpt_combine_kernels_cuda_rt::add_from_gpoint( - n_col, grid_cells.z, rt_flux_abs_dir.ptr(), rt_flux_abs_dif.ptr(), - (*fluxes).get_flux_abs_dir().ptr(), (*fluxes).get_flux_abs_dif().ptr()); - } + if (switch_raytracing) + { + Gpt_combine_kernels_cuda_rt::add_from_gpoint( + grid_cells.x, grid_cells.y, rt_flux_tod_up.ptr(), rt_flux_sfc_dir.ptr(), rt_flux_sfc_dif.ptr(), rt_flux_sfc_up.ptr(), + (*fluxes).get_flux_tod_up().ptr(), (*fluxes).get_flux_sfc_dir().ptr(), (*fluxes).get_flux_sfc_dif().ptr(), (*fluxes).get_flux_sfc_up().ptr()); - if (switch_single_gpt && igpt == single_gpt) - { - sw_gpt_flux_up = (*fluxes).get_flux_up(); - sw_gpt_flux_dn = (*fluxes).get_flux_dn(); - sw_gpt_flux_dn_dir = (*fluxes).get_flux_dn_dir(); - sw_gpt_flux_net = (*fluxes).get_flux_net(); - } + Gpt_combine_kernels_cuda_rt::add_from_gpoint( + n_col, grid_cells.z, rt_flux_abs_dir.ptr(), rt_flux_abs_dif.ptr(), + (*fluxes).get_flux_abs_dir().ptr(), (*fluxes).get_flux_abs_dif().ptr()); } + previous_band = band; } } diff --git a/src_test/test_rte_rrtmgp.cpp b/src_test/test_rte_rrtmgp.cpp index 2f8a1aa..dfb7c5b 100644 --- a/src_test/test_rte_rrtmgp.cpp +++ b/src_test/test_rte_rrtmgp.cpp @@ -137,14 +137,18 @@ void solve_radiation(int argc, char** argv) const auto settings = toml::parse(case_name + ".ini"); ////// FLOW CONTROL SWITCHES ////// - const bool switch_shortwave = get_ini_value(settings, "switches", "shortwave", true); - const bool switch_longwave = get_ini_value(settings, "switches", "longwave", true); - const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); - const bool switch_liquid_cloud_optics = get_ini_value(settings, "switches", "liquid_cloud_optics", false); - const bool switch_ice_cloud_optics = get_ini_value(settings, "switches", "ice_cloud_optics", false); - const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol_optics", false); - const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta_cloud", true); - const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta_aerosol", false); + const bool switch_liquid_cloud_optics = get_ini_value(settings, "clouds", "liquid_cloud_optics", false); + const bool switch_ice_cloud_optics = get_ini_value(settings, "clouds", "ice_cloud_optics", false); + const bool switch_delta_cloud = get_ini_value(settings, "clouds", "delta_cloud", false); + + const bool switch_shortwave = get_ini_value(settings, "shortwave", "shortwave", true); + const bool switch_sw_fluxes = get_ini_value(settings, "shortwave", "plane_parallel", true); + + const bool switch_longwave = get_ini_value(settings, "longwave", "longwave", true); + const bool switch_lw_fluxes = get_ini_value(settings, "longwave", "plane_parallel", true); + + const bool switch_aerosol_optics = get_ini_value(settings, "aerosols", "aerosol_optics", false); + const bool switch_delta_aerosol = get_ini_value(settings, "aerosols", "delta_aerosol", false); ////// READ THE ATMOSPHERIC DATA ////// @@ -289,7 +293,7 @@ void solve_radiation(int argc, char** argv) Array lw_flux_dn; Array lw_flux_net; - if (switch_fluxes) + if (switch_lw_fluxes) { lw_flux_up .set_dims({n_col, n_lev}); lw_flux_dn .set_dims({n_col, n_lev}); @@ -302,7 +306,7 @@ void solve_radiation(int argc, char** argv) auto time_start = std::chrono::high_resolution_clock::now(); rad_lw.solve( - switch_fluxes, + switch_lw_fluxes, switch_cloud_optics, gas_concs, p_lay, p_lev, @@ -328,7 +332,7 @@ void solve_radiation(int argc, char** argv) auto nc_lw_band_lims_wvn = output_nc.add_variable("lw_band_lims_wvn", {"band_lw", "pair"}); nc_lw_band_lims_wvn.insert(rad_lw.get_band_lims_wavenumber().v(), {0, 0}); - if (switch_fluxes) + if (switch_lw_fluxes) { auto nc_lw_flux_up = output_nc.add_variable("lw_flux_up" , {"lev", "y", "x"}); auto nc_lw_flux_dn = output_nc.add_variable("lw_flux_dn" , {"lev", "y", "x"}); @@ -383,7 +387,7 @@ void solve_radiation(int argc, char** argv) Array sw_flux_dn_dir; Array sw_flux_net; - if (switch_fluxes) + if (switch_sw_fluxes) { sw_flux_up .set_dims({n_col, n_lev}); sw_flux_dn .set_dims({n_col, n_lev}); @@ -397,7 +401,7 @@ void solve_radiation(int argc, char** argv) auto time_start = std::chrono::high_resolution_clock::now(); rad_sw.solve( - switch_fluxes, + switch_sw_fluxes, switch_cloud_optics, switch_aerosol_optics, switch_delta_cloud, @@ -430,7 +434,7 @@ void solve_radiation(int argc, char** argv) auto nc_sw_band_lims_wvn = output_nc.add_variable("sw_band_lims_wvn", {"band_sw", "pair"}); nc_sw_band_lims_wvn.insert(rad_sw.get_band_lims_wavenumber().v(), {0, 0}); - if (switch_fluxes) + if (switch_sw_fluxes) { auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); auto nc_sw_flux_dn = output_nc.add_variable("sw_flux_dn" , {"lev", "y", "x"}); diff --git a/src_test/test_rte_rrtmgp.cu b/src_test/test_rte_rrtmgp.cu index 0bb86ca..57c1aca 100644 --- a/src_test/test_rte_rrtmgp.cu +++ b/src_test/test_rte_rrtmgp.cu @@ -163,15 +163,18 @@ void solve_radiation(int argc, char** argv) const auto settings = toml::parse(case_name + ".ini"); ////// FLOW CONTROL SWITCHES ////// - const bool switch_shortwave = get_ini_value(settings, "switches", "shortwave", true); - const bool switch_longwave = get_ini_value(settings, "switches", "longwave", true); - const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); - const bool switch_liquid_cloud_optics = get_ini_value(settings, "switches", "liquid_cloud_optics", false); - const bool switch_ice_cloud_optics = get_ini_value(settings, "switches", "ice_cloud_optics", false); - const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol_optics", false); - const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta_cloud", true); - const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta_aerosol", false); + const bool switch_liquid_cloud_optics = get_ini_value(settings, "clouds", "liquid_cloud_optics", false); + const bool switch_ice_cloud_optics = get_ini_value(settings, "clouds", "ice_cloud_optics", false); + const bool switch_delta_cloud = get_ini_value(settings, "clouds", "delta_cloud", false); + const bool switch_shortwave = get_ini_value(settings, "shortwave", "shortwave", true); + const bool switch_sw_fluxes = get_ini_value(settings, "shortwave", "plane_parallel", true); + + const bool switch_longwave = get_ini_value(settings, "longwave", "longwave", true); + const bool switch_lw_fluxes = get_ini_value(settings, "longwave", "plane_parallel", true); + + const bool switch_aerosol_optics = get_ini_value(settings, "aerosols", "aerosol_optics", false); + const bool switch_delta_aerosol = get_ini_value(settings, "aerosols", "delta_aerosol", false); ////// READ THE ATMOSPHERIC DATA ////// Status::print_message("Reading atmospheric input data from NetCDF."); @@ -336,7 +339,7 @@ void solve_radiation(int argc, char** argv) Array_gpu lw_flux_dn; Array_gpu lw_flux_net; - if (switch_fluxes) + if (switch_lw_fluxes) { lw_flux_up .set_dims({n_col, n_lev}); lw_flux_dn .set_dims({n_col, n_lev}); @@ -370,7 +373,7 @@ void solve_radiation(int argc, char** argv) cudaEventRecord(start, 0); rad_lw.solve_gpu( - switch_fluxes, + switch_lw_fluxes, switch_cloud_optics, gas_concs_gpu, p_lay_gpu, p_lev_gpu, @@ -407,7 +410,7 @@ void solve_radiation(int argc, char** argv) auto nc_lw_band_lims_wvn = output_nc.add_variable("lw_band_lims_wvn", {"band_lw", "pair"}); nc_lw_band_lims_wvn.insert(rad_lw.get_band_lims_wavenumber_gpu().v(), {0, 0}); - if (switch_fluxes) + if (switch_lw_fluxes) { auto nc_lw_flux_up = output_nc.add_variable("lw_flux_up" , {"lev", "y", "x"}); auto nc_lw_flux_dn = output_nc.add_variable("lw_flux_dn" , {"lev", "y", "x"}); @@ -464,7 +467,7 @@ void solve_radiation(int argc, char** argv) Array_gpu sw_flux_dn_dir; Array_gpu sw_flux_net; - if (switch_fluxes) + if (switch_sw_fluxes) { sw_flux_up .set_dims({n_col, n_lev}); sw_flux_dn .set_dims({n_col, n_lev}); @@ -503,7 +506,7 @@ void solve_radiation(int argc, char** argv) cudaEventRecord(start, 0); rad_sw.solve_gpu( - switch_fluxes, + switch_sw_fluxes, switch_cloud_optics, switch_aerosol_optics, switch_delta_cloud, @@ -548,7 +551,7 @@ void solve_radiation(int argc, char** argv) auto nc_sw_band_lims_wvn = output_nc.add_variable("sw_band_lims_wvn", {"band_sw", "pair"}); nc_sw_band_lims_wvn.insert(rad_sw.get_band_lims_wavenumber_gpu().v(), {0, 0}); - if (switch_fluxes) + if (switch_sw_fluxes) { auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); auto nc_sw_flux_dn = output_nc.add_variable("sw_flux_dn" , {"lev", "y", "x"}); diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index 4dd5911..8d597dd 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -166,49 +166,101 @@ void solve_radiation(int argc, char** argv) const auto settings = toml::parse(case_name + ".ini"); - ////// FLOW CONTROL SWITCHES ////// - const bool switch_shortwave = get_ini_value(settings, "switches", "shortwave", true); - const bool switch_longwave = get_ini_value(settings, "switches", "longwave", true); - const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); - bool switch_sw_twostream = get_ini_value(settings, "switches", "sw_two_stream", false); - bool switch_sw_raytracing = get_ini_value(settings, "switches", "sw_raytracing", true); - bool switch_lw_raytracing = get_ini_value(settings, "switches", "lw_raytracing", true); - bool switch_independent_column = get_ini_value(settings, "switches", "independent_column", false); - bool switch_liquid_cloud_optics = get_ini_value(settings, "switches", "liquid_cloud_optics", false); - bool switch_ice_cloud_optics = get_ini_value(settings, "switches", "ice_cloud_optics", false); - const bool switch_lw_scattering = get_ini_value(settings, "switches", "lw_scattering", false); - const bool switch_min_mfp_grid_ratio= get_ini_value(settings, "switches", "min_mfp_grid_ratio", true); - const bool switch_cloud_mie = get_ini_value(settings, "switches", "cloud_mie", false); - const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol_optics", false); - const bool switch_single_gpt = get_ini_value(settings, "switches", "single_gpt", false); - const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta_cloud", false); - const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta_aerosol", false); - const bool switch_tica = get_ini_value(settings, "switches", "tica", false); - - const bool switch_bw_raytracing = get_ini_value(settings, "switches", "bw_raytracing", true); - const bool switch_lu_albedo = get_ini_value(settings, "switches", "lu_albedo", false); - const bool switch_image = get_ini_value(settings, "switches", "image", true); - const bool switch_broadband = get_ini_value(settings, "switches", "broadband", false); - const bool switch_cloud_cam = get_ini_value(settings, "switches", "cloud_cam", false); - - const int photons_per_pixel = get_ini_value(settings, "ints", "bw_raytracing", 1); - - if (!switch_shortwave) - switch_sw_raytracing = false; - - if (!switch_longwave) - switch_lw_raytracing = false; - - if (switch_shortwave && !switch_sw_twostream && !switch_sw_raytracing) + ////// READ INI FILE ////// + const bool switch_liquid_cloud_optics = get_ini_value(settings, "clouds", "liquid_cloud_optics", false); + const bool switch_ice_cloud_optics = get_ini_value(settings, "clouds", "ice_cloud_optics", false); + const bool switch_delta_cloud = get_ini_value(settings, "clouds", "delta_cloud", false); + const bool switch_tilted_columns = get_ini_value(settings, "clouds", "tilted_columns", false); + + // use mie scattering. Currently only for shortwave and only for liquid cloud droplets (requires ice_cloud_optics = false) + const bool switch_cloud_mie = get_ini_value(settings, "clouds", "cloud_mie", false); + + const bool switch_aerosol_optics = get_ini_value(settings, "aerosols", "aerosol_optics", false); + const bool switch_delta_aerosol = get_ini_value(settings, "aerosols", "delta_aerosol", false); + + const bool switch_shortwave = get_ini_value(settings, "shortwave", "shortwave", true); + // compute and output ray tracer fluxes + const bool switch_sw_raytracing = get_ini_value(settings, "shortwave", "raytracing", true); + // compute and output plane parallel 1D fluxes (two-stream) + const bool switch_sw_plane_parallel = get_ini_value(settings, "shortwave", "plane_parallel", true); + // solve ray tracer in independent column mode + bool switch_sw_independent_column = get_ini_value(settings, "shortwave", "rt_independent_column", false); + // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) + const int sw_single_gpt = get_ini_value(settings, "shortwave", "single_gpt", 0); + + const bool switch_longwave = get_ini_value(settings, "longwave", "longwave", true); + // compute and output ray tracer fluxes + const bool switch_lw_raytracing = get_ini_value(settings, "longwave", "raytracing", true); + // compute and output plane parallel 1D fluxes (two-stream or no-scattering solution) + const bool switch_lw_plane_parallel = get_ini_value(settings, "longwave", "plane_parallel", true); + // enable scattering in longwave solver (only possible in ray tracer executable) + const bool switch_lw_scattering = get_ini_value(settings, "longwave", "scattering", false); + // solve ray tracer in independent column mode + bool switch_lw_independent_column = get_ini_value(settings, "longwave", "rt_independent_column", false); + // only solve a single g-point and also output optical properties for that g-point. Defaults to -1 (broadband) + const int lw_single_gpt = get_ini_value(settings, "longwave", "single_gpt", 0); + // minimum ratio between the lowest gaseous mean free path and the horizontal grid spacing at which ray tracer is still used. Set to 0. to use ray tracer for all g-points + const Float min_mfp_grid_ratio = get_ini_value(settings, "longwave", "min_mfp_grid_ratio", Float(1.0)); + + const bool switch_bw_raytracing = get_ini_value(settings, "backward", "bw_raytracing", true); + // read surface type from land_use_map variable and compute spectral albedo and reflection type (lambertian/specular) accordingly + const bool switch_lu_albedo = get_ini_value(settings, "backward", "lu_albedo", false); + // solve visible bands and convert to XYZ tristimulus values + const bool switch_image = get_ini_value(settings, "backward", "image", true); + // solve broadband radiances + const bool switch_broadband = get_ini_value(settings, "backward", "broadband", false); + // output additional cloud statistics for each camera pixel + const bool switch_cloud_cam = get_ini_value(settings, "backward", "cloud_cam", false); + + // if >0, overwrite zenith angle from input netcdf file + const Float input_sza = get_ini_value(settings, "solar_angles", "sza", -1.0); + // if >0, overwrite azimuth angle from input netcdf file + const Float input_azi = get_ini_value(settings, "solar_angles", "azi", -1.0); + + Camera camera; + if (switch_bw_raytracing) + { + camera.fov = get_ini_value(settings, "camera", "cam_field_of_view", 80.); + + // camera type: (0) fish eye camera, (1) rectangular camera, (2) top-of-atmosphere upwelling radiances + camera.cam_type = get_ini_value(settings, "camera", "cam_type", 0); + + // x,y,z positions of virtual camera + camera.position = {get_ini_value(settings, "camera", "cam_px", 0.), + get_ini_value(settings, "camera", "cam_py", 0.), + get_ini_value(settings, "camera", "cam_pz", 100.)}; + + // width, height (pixels) of virtual camera or number of zenith and azimuth angles of fish camera + camera.nx = get_ini_value(settings, "camera", "cam_nx", 0); + camera.ny = get_ini_value(settings, "camera", "cam_ny", 0); + + // yaw, pitch and roll angles (degrees) of the virtual camera + camera.setup_rotation_matrix(get_ini_value(settings, "camera", "cam_yaw", 0.), + get_ini_value(settings, "camera", "cam_pitch", 0.), + get_ini_value(settings, "camera", "cam_roll", 0.)); + camera.setup_normal_camera(camera); + + camera.npix = Int(camera.nx * camera.ny); + } + + + if (switch_shortwave && !(switch_sw_plane_parallel || switch_sw_raytracing)) + { + std::string error = "With shortwave enabled, need to run shortwave plane parallel solver and/or ray tracer "; + throw std::runtime_error(error); + } + + if (switch_longwave && !(switch_lw_plane_parallel || switch_lw_raytracing)) { - std::string error = "With shortwave enable, need to run either two-stream solver or ray tracer "; + std::string error = "With longwave enabled, need to run longwave plane parallel solver and/or ray tracer "; throw std::runtime_error(error); } + // read samples counts (if applicable) Int sw_photons_per_pixel; if (switch_sw_raytracing) { - sw_photons_per_pixel = get_ini_value(settings, "ints", "sw_raytracing", Int(256)); + sw_photons_per_pixel = get_ini_value(settings, "shortwave", "samples", Int(256)); if (Float(int(std::log2(Float(sw_photons_per_pixel)))) != std::log2(Float(sw_photons_per_pixel))) { std::string error = "number of photons per pixel should be a power of 2 "; @@ -216,12 +268,15 @@ void solve_radiation(int argc, char** argv) } } + Int bw_photons_per_pixel; + if (switch_bw_raytracing) + bw_photons_per_pixel = get_ini_value(settings, "backward", "samples", Int(1)); Int lw_photon_power; Int lw_photon_count; if (switch_lw_raytracing) { - lw_photon_power = get_ini_value(settings, "ints", "lw_raytracing", Int(22)); + lw_photon_power = get_ini_value(settings, "longwave", "samples", Int(22)); lw_photon_count = 1 << lw_photon_power; } @@ -230,8 +285,11 @@ void solve_radiation(int argc, char** argv) Status::print_warning("No longwave radiation implemented in the backward ray tracer"); } - if (switch_tica) - switch_independent_column = true; + if (switch_tilted_columns) + { + switch_sw_independent_column = true; + switch_lw_independent_column = true; + } if (switch_cloud_cam && !(switch_liquid_cloud_optics || switch_ice_cloud_optics)) throw std::runtime_error("Enable cloud optics (liquid and/or ice) when cloud-cam is switch on!"); @@ -242,19 +300,14 @@ void solve_radiation(int argc, char** argv) throw std::runtime_error(error); } - int single_gpt = get_ini_value(settings, "ints", "single-gpt", 1); - - const Float min_mfp_grid_ratio = switch_min_mfp_grid_ratio ? get_ini_value(settings, "floats", "min_mfp_grid_ratio", Float(1.)) : Float(0.); - if (switch_sw_raytracing) Status::print_message("Shortwave: using "+ std::to_string(sw_photons_per_pixel) + " rays per pixel per g-point"); if (switch_lw_raytracing) Status::print_message("Longwave: using 2**"+std::to_string(lw_photon_power) + " ("+std::to_string(lw_photon_count) + ") rays per g-point"); - - const Float input_sza = get_ini_value(settings, "floats", "sza", -1.0); - const Float input_azi = get_ini_value(settings, "floats", "azi", -1.0); + if (switch_bw_raytracing) + Status::print_message("Backward: using "+ std::to_string(bw_photons_per_pixel) + " rays per pixel per g-point"); ////// READ THE ATMOSPHERIC DATA ////// Status::print_message("Reading atmospheric input data from NetCDF."); @@ -287,26 +340,6 @@ void solve_radiation(int argc, char** argv) input_nc.get_variable("ngrid_y"), input_nc.get_variable("ngrid_z")}; - Camera camera; - if (switch_bw_raytracing) - { - camera.fov = get_ini_value(settings, "bw_camera", "field-of-view", 80.); - camera.cam_type = get_ini_value(settings, "bw_camera", "cam-type", 0); - camera.position = {get_ini_value(settings, "bw_camera", "px", 0.), - get_ini_value(settings, "bw_camera", "py", 0.), - get_ini_value(settings, "bw_camera", "pz", 100.)}; - camera.nx = get_ini_value(settings, "bw_camera", "nx", 0); - camera.ny = get_ini_value(settings, "bw_camera", "ny", 0); - camera.setup_rotation_matrix(get_ini_value(settings, "bw_camera", "yaw", 0.), - get_ini_value(settings, "bw_camera", "pitch", 0.), - get_ini_value(settings, "bw_camera", "roll", 0.)); - camera.setup_normal_camera(camera); - - camera.npix = Int(camera.nx * camera.ny); - - const Float input_sza = get_ini_value(settings, "floats", "sza", -1.0); - const Float input_azi = get_ini_value(settings, "floats", "azi", -1.0); - } // Read the atmospheric fields. Array p_lay(input_nc.get_variable("p_lay", {n_lay, n_col_y, n_col_x}), {n_col, n_lay}); @@ -314,14 +347,14 @@ void solve_radiation(int argc, char** argv) Array p_lev(input_nc.get_variable("p_lev", {n_lev, n_col_y, n_col_x}), {n_col, n_lev}); Array t_lev(input_nc.get_variable("t_lev", {n_lev, n_col_y, n_col_x}), {n_col, n_lev}); - if (input_nc.variable_exists("col_dry") && switch_tica) + if (input_nc.variable_exists("col_dry") && switch_tilted_columns) { std::string error = "col_dry is not supported in tica mode"; throw std::runtime_error(error); } - if (input_nc.variable_exists("tsi") && switch_tica) + if (input_nc.variable_exists("tsi") && switch_tilted_columns) { std::string error = "tsi is overwritten in tica mode"; throw std::runtime_error(error); @@ -435,7 +468,7 @@ void solve_radiation(int argc, char** argv) azi.fill(input_azi / Float(180.0) * M_PI); - if (switch_tica) + if (switch_tilted_columns) { tica_sza = acos(mu0.v()[0]); tica_azi = azi.v()[0]; @@ -628,7 +661,7 @@ void solve_radiation(int argc, char** argv) Array_gpu lev_source; Array_gpu sfc_source; - if (switch_single_gpt) + if (lw_single_gpt > 0) { lw_tau_tot .set_dims({n_col, n_lay}); lw_tau_cld .set_dims({n_col, n_lay}); @@ -650,24 +683,13 @@ void solve_radiation(int argc, char** argv) Array_gpu lw_flux_dn; Array_gpu lw_flux_net; - if (switch_fluxes) + if (switch_lw_plane_parallel) { lw_flux_up .set_dims({n_col, n_lev}); lw_flux_dn .set_dims({n_col, n_lev}); lw_flux_net.set_dims({n_col, n_lev}); } - Array_gpu lw_gpt_flux_up; - Array_gpu lw_gpt_flux_dn; - Array_gpu lw_gpt_flux_net; - - if (switch_single_gpt) - { - lw_gpt_flux_up .set_dims({n_col, n_lev}); - lw_gpt_flux_dn .set_dims({n_col, n_lev}); - lw_gpt_flux_net.set_dims({n_col, n_lev}); - } - Array_gpu rt_flux_tod_up; Array_gpu rt_flux_tod_dn; Array_gpu rt_flux_sfc_up; @@ -718,16 +740,15 @@ void solve_radiation(int argc, char** argv) cudaEventRecord(start, 0); rad_lw.solve_gpu( - switch_fluxes, switch_lw_raytracing, + switch_lw_plane_parallel, switch_cloud_optics, switch_aerosol_optics, switch_delta_cloud, switch_delta_aerosol, - switch_single_gpt, switch_lw_scattering, - switch_independent_column, - single_gpt, + switch_lw_independent_column, + lw_single_gpt, min_mfp_grid_ratio, lw_photon_count, grid_cells, @@ -746,7 +767,6 @@ void solve_radiation(int argc, char** argv) lw_tau_aer, lw_ssa_aer, lw_asy_aer, lay_source, lev_source, sfc_source, lw_flux_up, lw_flux_dn, lw_flux_net, - lw_gpt_flux_up, lw_gpt_flux_dn, lw_gpt_flux_net, rt_flux_tod_up, rt_flux_tod_dn, rt_flux_sfc_up, rt_flux_sfc_dn, rt_flux_abs); @@ -788,9 +808,6 @@ void solve_radiation(int argc, char** argv) Array lw_flux_up_cpu(lw_flux_up); Array lw_flux_dn_cpu(lw_flux_dn); Array lw_flux_net_cpu(lw_flux_net); - Array lw_gpt_flux_up_cpu(lw_gpt_flux_up); - Array lw_gpt_flux_dn_cpu(lw_gpt_flux_dn); - Array lw_gpt_flux_net_cpu(lw_gpt_flux_net); Array rt_flux_tod_up_cpu(rt_flux_tod_up); Array rt_flux_tod_dn_cpu(rt_flux_tod_dn); @@ -801,7 +818,7 @@ void solve_radiation(int argc, char** argv) auto nc_lw_band_lims_wvn = output_nc.add_variable("lw_band_lims_wvn", {"band_lw", "pair"}); nc_lw_band_lims_wvn.insert(rad_lw.get_band_lims_wavenumber_gpu().v(), {0, 0}); - if (switch_single_gpt) + if (lw_single_gpt > 0) { auto nc_lw_band_lims_gpt = output_nc.add_variable("lw_band_lims_gpt", {"band_lw", "pair"}); nc_lw_band_lims_gpt.insert(rad_lw.get_band_lims_gpoint_gpu().v(), {0, 0}); @@ -844,7 +861,7 @@ void solve_radiation(int argc, char** argv) nc_sfc_source.insert(sfc_source_cpu.v(), {0, 0}); } - if (switch_fluxes) + if (switch_lw_plane_parallel) { auto nc_lw_flux_up = output_nc.add_variable("lw_flux_up" , {"lev", "y", "x"}); auto nc_lw_flux_dn = output_nc.add_variable("lw_flux_dn" , {"lev", "y", "x"}); @@ -853,32 +870,21 @@ 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_single_gpt) - { - auto nc_lw_gpt_flux_up = output_nc.add_variable("lw_gpt_flux_up" , {"lev", "y", "x"}); - auto nc_lw_gpt_flux_dn = output_nc.add_variable("lw_gpt_flux_dn" , {"lev", "y", "x"}); - auto nc_lw_gpt_flux_net = output_nc.add_variable("lw_gpt_flux_net", {"lev", "y", "x"}); - - nc_lw_gpt_flux_up .insert(lw_gpt_flux_up_cpu.v(), {0, 0, 0}); - nc_lw_gpt_flux_dn .insert(lw_gpt_flux_dn_cpu.v(), {0, 0, 0}); - nc_lw_gpt_flux_net.insert(lw_gpt_flux_net_cpu.v(), {0, 0, 0}); - } - - if (switch_lw_raytracing) - { - auto rt_flux_tod_up = output_nc.add_variable("rt_lw_flux_tod_up" , { "y", "x"}); - auto rt_flux_tod_dn = output_nc.add_variable("rt_lw_flux_tod_dn" , { "y", "x"}); - auto rt_flux_sfc_up = output_nc.add_variable("rt_lw_flux_sfc_up" , { "y", "x"}); - auto rt_flux_sfc_dn = output_nc.add_variable("rt_lw_flux_sfc_dn" , { "y", "x"}); - auto rt_flux_abs = output_nc.add_variable("rt_lw_flux_abs" , {"z", "y", "x"}); - - rt_flux_tod_up.insert(rt_flux_tod_up_cpu.v(), {0, 0}); - rt_flux_tod_dn.insert(rt_flux_tod_dn_cpu.v(), {0, 0}); - rt_flux_sfc_up.insert(rt_flux_sfc_up_cpu.v(), {0, 0}); - rt_flux_sfc_dn.insert(rt_flux_sfc_dn_cpu.v(), {0, 0}); - rt_flux_abs .insert(rt_flux_abs_cpu.v(), {0, 0, 0}); - } + if (switch_lw_raytracing) + { + auto rt_flux_tod_up = output_nc.add_variable("rt_lw_flux_tod_up" , { "y", "x"}); + auto rt_flux_tod_dn = output_nc.add_variable("rt_lw_flux_tod_dn" , { "y", "x"}); + auto rt_flux_sfc_up = output_nc.add_variable("rt_lw_flux_sfc_up" , { "y", "x"}); + auto rt_flux_sfc_dn = output_nc.add_variable("rt_lw_flux_sfc_dn" , { "y", "x"}); + auto rt_flux_abs = output_nc.add_variable("rt_lw_flux_abs" , {"z", "y", "x"}); + + rt_flux_tod_up.insert(rt_flux_tod_up_cpu.v(), {0, 0}); + rt_flux_tod_dn.insert(rt_flux_tod_dn_cpu.v(), {0, 0}); + rt_flux_sfc_up.insert(rt_flux_sfc_up_cpu.v(), {0, 0}); + rt_flux_sfc_dn.insert(rt_flux_sfc_dn_cpu.v(), {0, 0}); + rt_flux_abs .insert(rt_flux_abs_cpu.v(), {0, 0, 0}); } } @@ -925,7 +931,7 @@ void solve_radiation(int argc, char** argv) tsi_scaling({icol}) = Float(1.); } - if (switch_tica) + if (switch_tilted_columns) { for (int icol=1; icol<=n_col; ++icol) tsi_scaling({icol}) = std::cos(tica_sza); @@ -941,7 +947,7 @@ void solve_radiation(int argc, char** argv) Array_gpu sw_aer_ssa; Array_gpu sw_aer_asy; - if (switch_single_gpt) + if (sw_single_gpt > 0) { sw_tot_tau .set_dims({n_col, n_lay}); sw_tot_ssa .set_dims({n_col, n_lay}); @@ -966,42 +972,24 @@ void solve_radiation(int argc, char** argv) Array_gpu rt_flux_abs_dif; - if (switch_fluxes) + if(switch_sw_plane_parallel) { - if(switch_sw_twostream) - { - sw_flux_up .set_dims({n_col, n_lev}); - sw_flux_dn .set_dims({n_col, n_lev}); - sw_flux_dn_dir.set_dims({n_col, n_lev}); - sw_flux_net .set_dims({n_col, n_lev}); - } - - if (switch_sw_raytracing) - { - rt_flux_tod_up .set_dims({n_col_x, n_col_y}); - rt_flux_sfc_dir.set_dims({n_col_x, n_col_y}); - rt_flux_sfc_dif.set_dims({n_col_x, n_col_y}); - rt_flux_sfc_up .set_dims({n_col_x, n_col_y}); - rt_flux_abs_dir.set_dims({n_col_x, n_col_y, n_z}); - rt_flux_abs_dif.set_dims({n_col_x, n_col_y, n_z}); - } - + sw_flux_up .set_dims({n_col, n_lev}); + sw_flux_dn .set_dims({n_col, n_lev}); + sw_flux_dn_dir.set_dims({n_col, n_lev}); + sw_flux_net .set_dims({n_col, n_lev}); } - Array_gpu sw_gpt_flux_up; - Array_gpu sw_gpt_flux_dn; - Array_gpu sw_gpt_flux_dn_dir; - Array_gpu sw_gpt_flux_net; - - if (switch_single_gpt) + if (switch_sw_raytracing) { - sw_gpt_flux_up .set_dims({n_col, n_lev}); - sw_gpt_flux_dn .set_dims({n_col, n_lev}); - sw_gpt_flux_dn_dir.set_dims({n_col, n_lev}); - sw_gpt_flux_net .set_dims({n_col, n_lev}); + rt_flux_tod_up .set_dims({n_col_x, n_col_y}); + rt_flux_sfc_dir.set_dims({n_col_x, n_col_y}); + rt_flux_sfc_dif.set_dims({n_col_x, n_col_y}); + rt_flux_sfc_up .set_dims({n_col_x, n_col_y}); + rt_flux_abs_dir.set_dims({n_col_x, n_col_y, n_z}); + rt_flux_abs_dif.set_dims({n_col_x, n_col_y, n_z}); } - // Solve the radiation. Status::print_message("Solving the shortwave radiation."); @@ -1036,18 +1024,16 @@ void solve_radiation(int argc, char** argv) cudaEventRecord(start, 0); rad_sw.solve_gpu( - switch_fluxes, - switch_sw_twostream, switch_sw_raytracing, - switch_independent_column, + switch_sw_plane_parallel, + switch_sw_independent_column, switch_cloud_optics, switch_cloud_mie, switch_aerosol_optics, - switch_single_gpt, switch_delta_cloud, switch_delta_aerosol, - switch_tica, - single_gpt, + switch_tilted_columns, + sw_single_gpt, sw_photons_per_pixel, grid_cells, grid_d, @@ -1067,8 +1053,6 @@ void solve_radiation(int argc, char** argv) sw_aer_tau, sw_aer_ssa, sw_aer_asy, sw_flux_up, sw_flux_dn, sw_flux_dn_dir, sw_flux_net, - sw_gpt_flux_up, sw_gpt_flux_dn, - sw_gpt_flux_dn_dir, sw_gpt_flux_net, rt_flux_tod_up, rt_flux_sfc_dir, rt_flux_sfc_dif, @@ -1105,10 +1089,6 @@ void solve_radiation(int argc, char** argv) Array sw_flux_dn_cpu(sw_flux_dn); Array sw_flux_dn_dir_cpu(sw_flux_dn_dir); Array sw_flux_net_cpu(sw_flux_net); - Array sw_gpt_flux_up_cpu(sw_gpt_flux_up); - Array sw_gpt_flux_dn_cpu(sw_gpt_flux_dn); - Array sw_gpt_flux_dn_dir_cpu(sw_gpt_flux_dn_dir); - Array sw_gpt_flux_net_cpu(sw_gpt_flux_net); Array rt_flux_tod_up_cpu(rt_flux_tod_up); Array rt_flux_sfc_dir_cpu(rt_flux_sfc_dir); @@ -1120,7 +1100,7 @@ void solve_radiation(int argc, char** argv) auto nc_sw_band_lims_wvn = output_nc.add_variable("sw_band_lims_wvn", {"band_sw", "pair"}); nc_sw_band_lims_wvn.insert(rad_sw.get_band_lims_wavenumber_gpu().v(), {0, 0}); - if (switch_single_gpt) + if (sw_single_gpt > 0) { auto nc_sw_band_lims_gpt = output_nc.add_variable("sw_band_lims_gpt", {"band_sw", "pair"}); nc_sw_band_lims_gpt.insert(rad_sw.get_band_lims_gpoint_gpu().v(), {0, 0}); @@ -1143,14 +1123,14 @@ void solve_radiation(int argc, char** argv) nc_aer_ssa.insert(sw_aer_ssa_cpu.v(), {0, 0, 0}); nc_aer_asy.insert(sw_aer_asy_cpu.v(), {0, 0, 0}); - nc_tot_tau.add_attribute("long_name","Total optical depth at g-point "+std::to_string(single_gpt)); - nc_tot_ssa.add_attribute("long_name","Total single scattering albedo at g-point "+std::to_string(single_gpt)); - nc_cld_tau.add_attribute("long_name","Cloud optical depth at g-point "+std::to_string(single_gpt)); - nc_cld_ssa.add_attribute("long_name","Cloud single scattering albedo at g-point "+std::to_string(single_gpt)); - nc_cld_asy.add_attribute("long_name","Cloud asymmetry parameter at g-point "+std::to_string(single_gpt)); - nc_aer_tau.add_attribute("long_name","Aerosol optical depth at g-point "+std::to_string(single_gpt)); - nc_aer_ssa.add_attribute("long_name","Aerosol single scattering albedo at g-point "+std::to_string(single_gpt)); - nc_aer_asy.add_attribute("long_name","Aerosol asymmetry parameter at g-point "+std::to_string(single_gpt)); + nc_tot_tau.add_attribute("long_name","Total optical depth at g-point "+std::to_string(sw_single_gpt)); + nc_tot_ssa.add_attribute("long_name","Total single scattering albedo at g-point "+std::to_string(sw_single_gpt)); + nc_cld_tau.add_attribute("long_name","Cloud optical depth at g-point "+std::to_string(sw_single_gpt)); + nc_cld_ssa.add_attribute("long_name","Cloud single scattering albedo at g-point "+std::to_string(sw_single_gpt)); + nc_cld_asy.add_attribute("long_name","Cloud asymmetry parameter at g-point "+std::to_string(sw_single_gpt)); + nc_aer_tau.add_attribute("long_name","Aerosol optical depth at g-point "+std::to_string(sw_single_gpt)); + nc_aer_ssa.add_attribute("long_name","Aerosol single scattering albedo at g-point "+std::to_string(sw_single_gpt)); + nc_aer_asy.add_attribute("long_name","Aerosol asymmetry parameter at g-point "+std::to_string(sw_single_gpt)); nc_tot_tau.add_attribute("units", "-"); nc_tot_ssa.add_attribute("units", "-"); @@ -1163,99 +1143,65 @@ void solve_radiation(int argc, char** argv) } - if (switch_fluxes) + if (switch_sw_plane_parallel) { - if (switch_sw_twostream) - { - auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); - auto nc_sw_flux_dn = output_nc.add_variable("sw_flux_dn" , {"lev", "y", "x"}); - auto nc_sw_flux_dn_dir = output_nc.add_variable("sw_flux_dn_dir", {"lev", "y", "x"}); - auto nc_sw_flux_net = output_nc.add_variable("sw_flux_net" , {"lev", "y", "x"}); - - nc_sw_flux_up .insert(sw_flux_up_cpu .v(), {0, 0, 0}); - nc_sw_flux_dn .insert(sw_flux_dn_cpu .v(), {0, 0, 0}); - nc_sw_flux_dn_dir.insert(sw_flux_dn_dir_cpu.v(), {0, 0, 0}); - nc_sw_flux_net .insert(sw_flux_net_cpu .v(), {0, 0, 0}); - - nc_sw_flux_up.add_attribute("long_name","Upwelling shortwave fluxes (TwoStream solver)"); - nc_sw_flux_up.add_attribute("units","W m-2"); - - nc_sw_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes (TwoStream solver)"); - nc_sw_flux_dn.add_attribute("units","W m-2"); - - nc_sw_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes (TwoStream solver)"); - nc_sw_flux_dn_dir.add_attribute("units","W m-2"); - - nc_sw_flux_net.add_attribute("long_name","Net shortwave fluxes (TwoStream solver)"); - nc_sw_flux_net.add_attribute("units","W m-2"); - - } - - if (switch_sw_raytracing) - { - auto nc_rt_flux_tod_up = output_nc.add_variable("rt_flux_tod_up", {"y","x"}); - auto nc_rt_flux_sfc_dir = output_nc.add_variable("rt_flux_sfc_dir", {"y","x"}); - auto nc_rt_flux_sfc_dif = output_nc.add_variable("rt_flux_sfc_dif", {"y","x"}); - auto nc_rt_flux_sfc_up = output_nc.add_variable("rt_flux_sfc_up", {"y","x"}); - auto nc_rt_flux_abs_dir = output_nc.add_variable("rt_flux_abs_dir", {"z","y","x"}); - auto nc_rt_flux_abs_dif = output_nc.add_variable("rt_flux_abs_dif", {"z","y","x"}); + auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); + auto nc_sw_flux_dn = output_nc.add_variable("sw_flux_dn" , {"lev", "y", "x"}); + auto nc_sw_flux_dn_dir = output_nc.add_variable("sw_flux_dn_dir", {"lev", "y", "x"}); + auto nc_sw_flux_net = output_nc.add_variable("sw_flux_net" , {"lev", "y", "x"}); - nc_rt_flux_tod_up .insert(rt_flux_tod_up_cpu .v(), {0,0}); - nc_rt_flux_sfc_dir.insert(rt_flux_sfc_dir_cpu.v(), {0,0}); - nc_rt_flux_sfc_dif.insert(rt_flux_sfc_dif_cpu.v(), {0,0}); - nc_rt_flux_sfc_up .insert(rt_flux_sfc_up_cpu .v(), {0,0}); + nc_sw_flux_up .insert(sw_flux_up_cpu .v(), {0, 0, 0}); + nc_sw_flux_dn .insert(sw_flux_dn_cpu .v(), {0, 0, 0}); + nc_sw_flux_dn_dir.insert(sw_flux_dn_dir_cpu.v(), {0, 0, 0}); + nc_sw_flux_net .insert(sw_flux_net_cpu .v(), {0, 0, 0}); - nc_rt_flux_abs_dir.insert(rt_flux_abs_dir_cpu.v(), {0,0,0}); - nc_rt_flux_abs_dif.insert(rt_flux_abs_dif_cpu.v(), {0,0,0}); + nc_sw_flux_up.add_attribute("long_name","Upwelling shortwave fluxes (TwoStream solver)"); + nc_sw_flux_up.add_attribute("units","W m-2"); - nc_rt_flux_tod_up.add_attribute("long_name","Upwelling shortwave top-of-domain fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_tod_up.add_attribute("units","W m-2"); + nc_sw_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes (TwoStream solver)"); + nc_sw_flux_dn.add_attribute("units","W m-2"); - nc_rt_flux_sfc_dir.add_attribute("long_name","Downwelling direct shortwave surface fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_sfc_dir.add_attribute("units","W m-2"); + nc_sw_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes (TwoStream solver)"); + nc_sw_flux_dn_dir.add_attribute("units","W m-2"); - nc_rt_flux_sfc_dif.add_attribute("long_name","Downwelling diffuse shortwave surface fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_sfc_dif.add_attribute("units","W m-2"); + nc_sw_flux_net.add_attribute("long_name","Net shortwave fluxes (TwoStream solver)"); + nc_sw_flux_net.add_attribute("units","W m-2"); + } - nc_rt_flux_sfc_up.add_attribute("long_name","Upwelling shortwave surface fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_sfc_up.add_attribute("units","W m-2"); + if (switch_sw_raytracing) + { + auto nc_rt_flux_tod_up = output_nc.add_variable("rt_flux_tod_up", {"y","x"}); + auto nc_rt_flux_sfc_dir = output_nc.add_variable("rt_flux_sfc_dir", {"y","x"}); + auto nc_rt_flux_sfc_dif = output_nc.add_variable("rt_flux_sfc_dif", {"y","x"}); + auto nc_rt_flux_sfc_up = output_nc.add_variable("rt_flux_sfc_up", {"y","x"}); + auto nc_rt_flux_abs_dir = output_nc.add_variable("rt_flux_abs_dir", {"z","y","x"}); + auto nc_rt_flux_abs_dif = output_nc.add_variable("rt_flux_abs_dif", {"z","y","x"}); - nc_rt_flux_abs_dir.add_attribute("long_name","Absorbed direct shortwave fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_abs_dir.add_attribute("units","W m-3"); + nc_rt_flux_tod_up .insert(rt_flux_tod_up_cpu .v(), {0,0}); + nc_rt_flux_sfc_dir.insert(rt_flux_sfc_dir_cpu.v(), {0,0}); + nc_rt_flux_sfc_dif.insert(rt_flux_sfc_dif_cpu.v(), {0,0}); + nc_rt_flux_sfc_up .insert(rt_flux_sfc_up_cpu .v(), {0,0}); - nc_rt_flux_abs_dif.add_attribute("long_name","Absorbed diffuse shortwave fluxes (Monte Carlo ray tracer)"); - nc_rt_flux_abs_dif.add_attribute("units","W m-3"); + nc_rt_flux_abs_dir.insert(rt_flux_abs_dir_cpu.v(), {0,0,0}); + nc_rt_flux_abs_dif.insert(rt_flux_abs_dif_cpu.v(), {0,0,0}); - } + nc_rt_flux_tod_up.add_attribute("long_name","Upwelling shortwave top-of-domain fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_tod_up.add_attribute("units","W m-2"); + nc_rt_flux_sfc_dir.add_attribute("long_name","Downwelling direct shortwave surface fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_sfc_dir.add_attribute("units","W m-2"); - if (switch_single_gpt) - { - if (switch_sw_twostream) - { - auto nc_sw_gpt_flux_up = output_nc.add_variable("sw_gpt_flux_up" , {"lev", "y", "x"}); - auto nc_sw_gpt_flux_dn = output_nc.add_variable("sw_gpt_flux_dn" , {"lev", "y", "x"}); - auto nc_sw_gpt_flux_dn_dir = output_nc.add_variable("sw_gpt_flux_dn_dir", {"lev", "y", "x"}); - auto nc_sw_gpt_flux_net = output_nc.add_variable("sw_gpt_flux_net" , {"lev", "y", "x"}); + nc_rt_flux_sfc_dif.add_attribute("long_name","Downwelling diffuse shortwave surface fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_sfc_dif.add_attribute("units","W m-2"); - nc_sw_gpt_flux_up .insert(sw_gpt_flux_up_cpu .v(), {0, 0, 0}); - nc_sw_gpt_flux_dn .insert(sw_gpt_flux_dn_cpu .v(), {0, 0, 0}); - nc_sw_gpt_flux_dn_dir.insert(sw_gpt_flux_dn_dir_cpu.v(), {0, 0, 0}); - nc_sw_gpt_flux_net .insert(sw_gpt_flux_net_cpu .v(), {0, 0, 0}); + nc_rt_flux_sfc_up.add_attribute("long_name","Upwelling shortwave surface fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_sfc_up.add_attribute("units","W m-2"); - nc_sw_gpt_flux_up.add_attribute("long_name","Upwelling shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_up.add_attribute("units","W m-2"); + nc_rt_flux_abs_dir.add_attribute("long_name","Absorbed direct shortwave fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_abs_dir.add_attribute("units","W m-3"); - nc_sw_gpt_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_dn.add_attribute("units","W m-2"); - - nc_sw_gpt_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_dn_dir.add_attribute("units","W m-2"); - - nc_sw_gpt_flux_net.add_attribute("long_name","Net shortwave fluxes for g-point "+std::to_string(single_gpt)+" (TwoStream solver)"); - nc_sw_gpt_flux_net.add_attribute("units","W m-2"); - } - } + nc_rt_flux_abs_dif.add_attribute("long_name","Absorbed diffuse shortwave fluxes (Monte Carlo ray tracer)"); + nc_rt_flux_abs_dif.add_attribute("units","W m-3"); } } @@ -1374,7 +1320,7 @@ void solve_radiation(int argc, char** argv) grid_cells, grid_d, kn_grid, - photons_per_pixel, + bw_photons_per_pixel, gas_concs_gpu, p_lay_gpu, p_lev_gpu, t_lay_gpu, t_lev_gpu, @@ -1448,7 +1394,7 @@ void solve_radiation(int argc, char** argv) grid_cells, grid_d, kn_grid, - photons_per_pixel, + bw_photons_per_pixel, gas_concs_gpu, p_lay_gpu, p_lev_gpu, t_lay_gpu, t_lev_gpu,