diff --git a/config/snellius.cmake b/config/snellius.cmake index d79b4bc..4cf4f4e 100644 --- a/config/snellius.cmake +++ b/config/snellius.cmake @@ -68,8 +68,9 @@ set(LIBS ${NETCDF_LIB_C} ${HDF5_LIB} ${SZIP_LIB}) if(USECUDA) set(CUDA_PROPAGATE_HOST_FLAGS OFF) + set(CMAKE_CUDA_ARCHITECTURES 80) set(LIBS ${LIBS} -rdynamic $ENV{EBROOTCUDA}/lib64/libcufft.so) - set(USER_CUDA_FLAGS "-arch=sm_80 -std=c++14 --expt-relaxed-constexpr") + set(USER_CUDA_FLAGS "-arch=sm_80 -std=c++14 -O3 --expt-relaxed-constexpr") set(USER_CUDA_FLAGS_RELEASE "-DNDEBUG") add_definitions(-DRTE_RRTMGP_GPU_MEMPOOL_CUDA) endif() diff --git a/data/mie_lut_broadband.nc b/data/mie_lut_broadband.nc index 45e90ed..ab349d8 100644 Binary files a/data/mie_lut_broadband.nc and b/data/mie_lut_broadband.nc differ diff --git a/include_rt/Raytracer.h b/include_rt/Raytracer.h index 11e5981..9e93a7d 100644 --- a/include_rt/Raytracer.h +++ b/include_rt/Raytracer.h @@ -23,6 +23,7 @@ class Raytracer void trace_rays( const int qrng_gpt_offset, + const bool switch_independent_column, const Int photons_per_pixel, const Raytracer_definitions::Vector grid_cells, const Raytracer_definitions::Vector grid_d, diff --git a/include_rt/Raytracer_bw.h b/include_rt/Raytracer_bw.h index 0184c04..7f5f7d0 100644 --- a/include_rt/Raytracer_bw.h +++ b/include_rt/Raytracer_bw.h @@ -96,6 +96,18 @@ class Raytracer_bw const Float total_source, Array_gpu& XYZ); + void accumulate_clouds( + const Vector& grid_d, + const Vector& grid_cells, + const Array_gpu& lwp, + const Array_gpu& iwp, + const Array_gpu& tau_cloud, + const Camera& camera, + Array_gpu& liwp_cam, + Array_gpu& tauc_cam, + Array_gpu& dist_cam, + Array_gpu& zen_cam); + private: }; diff --git a/include_rt/raytracer_functions.h b/include_rt/raytracer_functions.h index 71b0d07..8e614d5 100644 --- a/include_rt/raytracer_functions.h +++ b/include_rt/raytracer_functions.h @@ -20,6 +20,8 @@ namespace Raytracer_functions template static inline __host__ __device__ Vector operator-(const Vector v1, const Vector v2) { return Vector{v1.x-v2.x, v1.y-v2.y, v1.z-v2.z}; } template static inline __host__ __device__ + Vector operator-(const Vector v, const Float s) { return Vector{v.x-s, v.y-s, v.z-s}; } + template static inline __host__ __device__ Vector operator+(const Vector v, const Float s) { return Vector{s+v.x, s+v.y, s+v.z}; } template static inline __host__ __device__ Vector operator+(const Vector v1, const Vector v2) { return Vector{v1.x+v2.x, v1.y+v2.y, v1.z+v2.z}; } diff --git a/include_rt_kernels/raytracer_kernels.h b/include_rt_kernels/raytracer_kernels.h index 0328c90..25bd436 100644 --- a/include_rt_kernels/raytracer_kernels.h +++ b/include_rt_kernels/raytracer_kernels.h @@ -20,6 +20,7 @@ constexpr Float k_null_gas_min = Float(1.e-3); __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_rt_kernels/raytracer_kernels_bw.h b/include_rt_kernels/raytracer_kernels_bw.h index f1c8a20..d69471c 100644 --- a/include_rt_kernels/raytracer_kernels_bw.h +++ b/include_rt_kernels/raytracer_kernels_bw.h @@ -25,6 +25,11 @@ constexpr int bw_kernel_grid = 256; #endif constexpr Float k_null_gas_min = Float(1.e-3); +// sun has a half angle of .266 degrees +constexpr Float cos_half_angle = Float(0.9999891776066407); // cos(half_angle); +constexpr Float sun_solid_angle = Float(6.799910294339209e-05); // 2.*M_PI*(1-cos_half_angle); +constexpr Float sun_solid_angle_reciprocal = Float(14706.07635563193); + struct Grid_knull { @@ -116,4 +121,19 @@ void ray_tracer_kernel_bw( const Float* __restrict__ mie_phase, const Float* __restrict__ mie_phase_ang, const int mie_table_size); + +__global__ +void accumulate_clouds_kernel( + const Float* __restrict__ lwp, + const Float* __restrict__ iwp, + const Float* __restrict__ tau_cloud, + const Vector grid_d, + const Vector grid_size, + const Vector grid_cells, + Float* __restrict__ liwp_cam, + Float* __restrict__ tauc_cam, + Float* __restrict__ dist_cam, + Float* __restrict__ zen_cam, + const Camera camera); + #endif diff --git a/include_test/Radiation_solver_bw.h b/include_test/Radiation_solver_bw.h index 3c0ff14..3072a1b 100644 --- a/include_test/Radiation_solver_bw.h +++ b/include_test/Radiation_solver_bw.h @@ -145,8 +145,10 @@ class Radiation_solver_shortwave Array& sw_bnd_flux_dn_dir, Array& sw_bnd_flux_net) const; void load_mie_tables( - const std::string& file_name_mie, - const bool switch_broadband); + const std::string& file_name_mie_bb, + const std::string& file_name_mie_im, + const bool switch_broadband, + const bool switch_image); #ifdef __CUDACC__ void solve_gpu( @@ -157,6 +159,8 @@ class Radiation_solver_shortwave const bool switch_lu_albedo, const bool switch_delta_cloud, const bool switch_delta_aerosol, + const bool switch_cloud_cam, + const bool switch_raytracing, const Vector& grid_cells, const Vector& grid_d, const Vector& kn_grid, @@ -175,7 +179,11 @@ class Radiation_solver_shortwave const Array_gpu& rh, const Aerosol_concs_gpu& aerosol_concs, const Camera& camera, - Array_gpu& XYZ); + Array_gpu& XYZ, + Array_gpu& liwp_cam, + Array_gpu& tauc_cam, + Array_gpu& dist_cam, + Array_gpu& zen_cam); #endif #ifdef __CUDACC__ @@ -186,6 +194,8 @@ class Radiation_solver_shortwave const bool switch_lu_albedo, const bool switch_delta_cloud, const bool switch_delta_aerosol, + const bool switch_cloud_cam, + const bool switch_raytracing, const Vector& grid_cells, const Vector& grid_d, const Vector& kn_grid, @@ -204,7 +214,11 @@ class Radiation_solver_shortwave const Array_gpu& rh, const Aerosol_concs_gpu& aerosol_concs, const Camera& camera, - Array_gpu& radiance); + Array_gpu& radiance, + Array_gpu& liwp_cam, + Array_gpu& tauc_cam, + Array_gpu& dist_cam, + Array_gpu& zen_cam); int get_n_gpt_gpu() const { return this->kdist_gpu->get_ngpt(); }; int get_n_bnd_gpu() const { return this->kdist_gpu->get_nband(); }; @@ -232,11 +246,16 @@ class Radiation_solver_shortwave std::unique_ptr cloud_optical_props; std::unique_ptr aerosol_optical_props; - Array_gpu mie_phase_angs; + Array_gpu mie_phase_angs_bb; + Array_gpu mie_cdfs_bb; + Array_gpu mie_angs_bb; + Array_gpu mie_phase_bb; + Array_gpu mie_phase_angs; Array_gpu mie_cdfs; Array_gpu mie_angs; Array_gpu mie_phase; + #endif }; #endif diff --git a/include_test/Radiation_solver_rt.h b/include_test/Radiation_solver_rt.h index 81d7878..bb56d27 100644 --- a/include_test/Radiation_solver_rt.h +++ b/include_test/Radiation_solver_rt.h @@ -100,13 +100,13 @@ class Radiation_solver_shortwave void solve_gpu( const bool switch_fluxes, const bool switch_raytracing, + const bool switch_independent_column, const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, const bool switch_single_gpt, const bool switch_delta_cloud, const bool switch_delta_aerosol, - const bool switch_clear_sky_tod, const int single_gpt, const Int ray_count, const Vector grid_cells, diff --git a/python/set_virtual_camera.py b/python/set_virtual_camera.py index 0f365fc..849eddc 100644 --- a/python/set_virtual_camera.py +++ b/python/set_virtual_camera.py @@ -61,11 +61,12 @@ cam["roll"][:] = 0 cam["fisheye"][:]= 1 cam["f_zoom"][:]= 1 + cam["fov"][:] = 80 cam["px"][:] = 0 cam["py"][:] = 0 cam["pz"][:] = 0 - cam["nx"][:] = 1024 - cam["ny"][:] = 1024 + cam["nx"][:] = 256 + cam["ny"][:] = 256 # example imagery settings if args['image']: @@ -73,12 +74,13 @@ cam["pitch"][:] = 0 cam["roll"][:] = 0 cam["fisheye"][:]= 0 + cam["f_zoom"][:]= 1 cam["fov"][:] = 80 cam["px"][:] = 0. cam["py"][:] = 0. cam["pz"][:] = 500. - cam["nx"][:] = 1024 - cam["ny"][:] = 1024 + cam["nx"][:] = 256 + cam["ny"][:] = 256 for var in camera_variables: if not args[var] is None: diff --git a/src_cuda_rt/Gas_optics_rrtmgp_rt.cu b/src_cuda_rt/Gas_optics_rrtmgp_rt.cu index 1e90b18..3cd02de 100644 --- a/src_cuda_rt/Gas_optics_rrtmgp_rt.cu +++ b/src_cuda_rt/Gas_optics_rrtmgp_rt.cu @@ -831,33 +831,38 @@ void Gas_optics_rrtmgp_rt::find_relevant_gases_gpt( const int minor_start_lower = first_last_minor_lower({1, igpt}); const int minor_end_lower = first_last_minor_lower({2, igpt}); - for (int imnr=minor_start_lower+1; imnr<=minor_end_lower+1; ++imnr) + if (minor_start_lower >= 0) { - const int minor_gas = idx_minor_lower({imnr}); - add_if_not_in_vector(gases, minor_gas); - - if ( (minor_scales_with_density_lower({imnr})) && (idx_minor_scaling_lower({imnr}) > 0)) + for (int imnr=minor_start_lower+1; imnr<=minor_end_lower+1; ++imnr) { - const int scale_gas = idx_minor_scaling_lower({imnr}); - add_if_not_in_vector(gases, scale_gas); + const int minor_gas = idx_minor_lower({imnr}); + add_if_not_in_vector(gases, minor_gas); + + if ( (minor_scales_with_density_lower({imnr})) && (idx_minor_scaling_lower({imnr}) > 0)) + { + const int scale_gas = idx_minor_scaling_lower({imnr}); + add_if_not_in_vector(gases, scale_gas); + } } } const int minor_start_upper = first_last_minor_upper({1, igpt}); const int minor_end_upper = first_last_minor_upper({2, igpt}); - for (int imnr=minor_start_upper+1; imnr<=minor_end_upper+1; ++imnr) + if (minor_start_upper >= 0) { - const int minor_gas = idx_minor_upper({imnr}); - add_if_not_in_vector(gases, minor_gas); - - if ( (minor_scales_with_density_upper({imnr})) && (idx_minor_scaling_upper({imnr}) > 0)) + for (int imnr=minor_start_upper+1; imnr<=minor_end_upper+1; ++imnr) { - const int scale_gas = idx_minor_scaling_upper({imnr}); - add_if_not_in_vector(gases, scale_gas); + const int minor_gas = idx_minor_upper({imnr}); + add_if_not_in_vector(gases, minor_gas); + + if ( (minor_scales_with_density_upper({imnr})) && (idx_minor_scaling_upper({imnr}) > 0)) + { + const int scale_gas = idx_minor_scaling_upper({imnr}); + add_if_not_in_vector(gases, scale_gas); + } } } - } __global__ diff --git a/src_cuda_rt/Raytracer.cu b/src_cuda_rt/Raytracer.cu index 38aca01..10eb3d1 100644 --- a/src_cuda_rt/Raytracer.cu +++ b/src_cuda_rt/Raytracer.cu @@ -86,8 +86,8 @@ namespace const int icol_x = blockIdx.x*blockDim.x + threadIdx.x; const int icol_y = blockIdx.y*blockDim.y + threadIdx.y; const int iz = blockIdx.z*blockDim.z + threadIdx.z; - - if ( (icol_x < grid_cells.x) && (icol_y < grid_cells.y) && (iz < grid_cells.z) ) + + if ( (icol_x < grid_cells.x) && (icol_y < grid_cells.y) && (iz < (grid_cells.z - 1)) ) { const int idx = icol_x + icol_y*grid_cells.x + iz*grid_cells.y*grid_cells.x; const Float kext_tot = tau_tot[idx] / grid_d.z; @@ -106,6 +106,61 @@ namespace } } + __global__ + void bundles_optical_props_tod( + const Vector grid_cells, const Vector grid_d, const int n_lev, + const Float* __restrict__ tau_tot, const Float* __restrict__ ssa_tot, + const Float* __restrict__ tau_cld, const Float* __restrict__ ssa_cld, const Float* __restrict__ asy_cld, + const Float* __restrict__ tau_aer, const Float* __restrict__ ssa_aer, const Float* __restrict__ asy_aer, + Float* __restrict__ k_ext, Optics_scat* __restrict__ scat_asy) + { + const int icol_x = blockIdx.x*blockDim.x + threadIdx.x; + const int icol_y = blockIdx.y*blockDim.y + threadIdx.y; + + const int z_tod = grid_cells.z - 1; + + if ( (icol_x < grid_cells.x) && (icol_y < grid_cells.y)) + { + Float tau_tot_sum = Float(0.); + Float tausca_tot_sum = Float(0.); + + Float tausca_cld_sum = Float(0.); + Float tauscag_cld_sum = Float(0.); + + Float tausca_aer_sum = Float(0.); + Float tauscag_aer_sum = Float(0.); + + for (int iz=z_tod; iz grid_cells, const Vector grid_d, @@ -197,6 +253,8 @@ void Raytracer::trace_rays( Array_gpu& flux_abs_dir, Array_gpu& flux_abs_dif) { + const int n_lay = tau_total.dim(2); + // set of block and grid dimensions used in data processing kernels - requires some proper tuning later const int block_col_x = 8; const int block_col_y = 8; @@ -215,12 +273,21 @@ void Raytracer::trace_rays( Array_gpu k_ext({grid_cells.x, grid_cells.y, grid_cells.z}); Array_gpu scat_asy({grid_cells.x, grid_cells.y, grid_cells.z}); + // first on the whole grid expect the extra layer bundles_optical_props<<>>( grid_cells, grid_d, tau_total.ptr(), ssa_total.ptr(), tau_cloud.ptr(), ssa_cloud.ptr(), asy_cloud.ptr(), tau_aeros.ptr(), ssa_aeros.ptr(), asy_aeros.ptr(), k_ext.ptr(), scat_asy.ptr()); + + // second, integrate from TOD to TOA + bundles_optical_props_tod<<>>( + grid_cells, grid_d, n_lay, + tau_total.ptr(), ssa_total.ptr(), + tau_cloud.ptr(), ssa_cloud.ptr(), asy_cloud.ptr(), + tau_aeros.ptr(), ssa_aeros.ptr(), asy_aeros.ptr(), + k_ext.ptr(), scat_asy.ptr()); // create k_null_grid const int block_kn_x = 8; @@ -303,6 +370,7 @@ void Raytracer::trace_rays( 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, diff --git a/src_cuda_rt/Raytracer_bw.cu b/src_cuda_rt/Raytracer_bw.cu index 49c7471..6b8813f 100644 --- a/src_cuda_rt/Raytracer_bw.cu +++ b/src_cuda_rt/Raytracer_bw.cu @@ -655,9 +655,59 @@ void Raytracer_bw::trace_rays_bb( count_to_flux_2d<<>>( camera, photons_per_pixel, toa_src, - Float(1.), + sun_solid_angle_reciprocal, camera_count.ptr(), flux_camera.ptr()); } +void Raytracer_bw::accumulate_clouds( + const Vector& grid_d, + const Vector& grid_cells, + const Array_gpu& lwp, + const Array_gpu& iwp, + const Array_gpu& tau_cloud, + const Camera& camera, + Array_gpu& liwp_cam, + Array_gpu& tauc_cam, + Array_gpu& dist_cam, + Array_gpu& zen_cam) +{ + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, liwp_cam.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, tauc_cam.ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, dist_cam.ptr()); + + const int n_pix = camera.nx * camera.ny; + const int n_block = std::min(n_pix, 512); + const int n_grid = std::ceil(Float(n_pix)/n_block); + + dim3 grid(n_grid); + dim3 block(n_block); + + // domain sizes + const Vector grid_size = grid_d * grid_cells; + + accumulate_clouds_kernel<<>>( + lwp.ptr(), + iwp.ptr(), + tau_cloud.ptr(), + grid_d, + grid_size, + grid_cells, + liwp_cam.ptr(), + tauc_cam.ptr(), + dist_cam.ptr(), + zen_cam.ptr(), + camera); + +} + + + + + + + + + + diff --git a/src_kernels_cuda_rt/raytracer_kernels.cu b/src_kernels_cuda_rt/raytracer_kernels.cu index 4960ac0..a4e3ca3 100644 --- a/src_kernels_cuda_rt/raytracer_kernels.cu +++ b/src_kernels_cuda_rt/raytracer_kernels.cu @@ -119,6 +119,7 @@ namespace __global__ void ray_tracer_kernel( + const Bool independent_column, const Int photons_to_shoot, const Int qrng_grid_x, const Int qrng_grid_y, @@ -199,7 +200,7 @@ void ray_tracer_kernel( const Float sx = abs((photon.direction.x > 0) ? ((i_n+1) * kn_grid_d.x - photon.position.x)/photon.direction.x : (i_n*kn_grid_d.x - photon.position.x)/photon.direction.x); const Float sy = abs((photon.direction.y > 0) ? ((j_n+1) * kn_grid_d.y - photon.position.y)/photon.direction.y : (j_n*kn_grid_d.y - photon.position.y)/photon.direction.y); const Float sz = abs((photon.direction.z > 0) ? ((k_n+1) * kn_grid_d.z - photon.position.z)/photon.direction.z : (k_n*kn_grid_d.z - photon.position.z)/photon.direction.z); - d_max = min(sx, min(sy, sz)); + d_max = independent_column ? sz : min(sx, min(sy, sz)); const int ijk_n = i_n + j_n*kn_grid.x + k_n*kn_grid.x*kn_grid.y; k_ext_null = k_null_grid[ijk_n]; } @@ -213,12 +214,16 @@ void ray_tracer_kernel( if (dn >= d_max) { - const Float dx = photon.direction.x * (s_min + d_max); - const Float dy = photon.direction.y * (s_min + d_max); - const Float dz = photon.direction.z * (s_min + d_max); + if (!independent_column) + { + const Float dx = photon.direction.x * (s_min + d_max); + const Float dy = photon.direction.y * (s_min + d_max); - photon.position.x += dx; - photon.position.y += dy; + photon.position.x += dx; + photon.position.y += dy; + } + + const Float dz = photon.direction.z * (s_min + d_max); photon.position.z += dz; // surface hit @@ -309,20 +314,22 @@ void ray_tracer_kernel( // regular cell crossing: adjust tau and apply periodic BC else { - photon.position.x += photon.direction.x>0 ? s_min : -s_min; - photon.position.y += photon.direction.y>0 ? s_min : -s_min; photon.position.z += photon.direction.z>0 ? s_min : -s_min; - - // Cyclic boundary condition in x. - photon.position.x = fmod(photon.position.x, grid_size.x); - if (photon.position.x < Float(0.)) - photon.position.x += grid_size.x; - - // Cyclic boundary condition in y. - photon.position.y = fmod(photon.position.y, grid_size.y); - if (photon.position.y < Float(0.)) - photon.position.y += grid_size.y; - + if (!independent_column) + { + photon.position.x += photon.direction.x>0 ? s_min : -s_min; + photon.position.y += photon.direction.y>0 ? s_min : -s_min; + + // Cyclic boundary condition in x. + photon.position.x = fmod(photon.position.x, grid_size.x); + if (photon.position.x < Float(0.)) + photon.position.x += grid_size.x; + + // Cyclic boundary condition in y. + photon.position.y = fmod(photon.position.y, grid_size.y); + if (photon.position.y < Float(0.)) + photon.position.y += grid_size.y; + } tau -= d_max * k_ext_null; d_max = Float(0.); transition = true; @@ -330,13 +337,17 @@ void ray_tracer_kernel( } else { - Float dx = photon.direction.x * dn; - Float dy = photon.direction.y * dn; Float dz = photon.direction.z * dn; - - photon.position.x = (dx > 0) ? min(photon.position.x + dx, (i_n+1) * kn_grid_d.x - s_min) : max(photon.position.x + dx, (i_n) * kn_grid_d.x + s_min); - photon.position.y = (dy > 0) ? min(photon.position.y + dy, (j_n+1) * kn_grid_d.y - s_min) : max(photon.position.y + dy, (j_n) * kn_grid_d.y + s_min); photon.position.z = (dz > 0) ? min(photon.position.z + dz, (k_n+1) * kn_grid_d.z - s_min) : max(photon.position.z + dz, (k_n) * kn_grid_d.z + s_min); + + if (!independent_column) + { + Float dx = photon.direction.x * dn; + Float dy = photon.direction.y * dn; + + photon.position.x = (dx > 0) ? min(photon.position.x + dx, (i_n+1) * kn_grid_d.x - s_min) : max(photon.position.x + dx, (i_n) * kn_grid_d.x + s_min); + photon.position.y = (dy > 0) ? min(photon.position.y + dy, (j_n+1) * kn_grid_d.y - s_min) : max(photon.position.y + dy, (j_n) * kn_grid_d.y + s_min); + } // Calculate the 3D index. const int i = float_to_int(photon.position.x, grid_d.x, grid_cells.x); diff --git a/src_kernels_cuda_rt/raytracer_kernels_bw.cu b/src_kernels_cuda_rt/raytracer_kernels_bw.cu index 47bb77b..abc2747 100644 --- a/src_kernels_cuda_rt/raytracer_kernels_bw.cu +++ b/src_kernels_cuda_rt/raytracer_kernels_bw.cu @@ -12,10 +12,6 @@ namespace constexpr Float w_thres = 0.5; - // sun has a half angle of .266 degrees - constexpr Float cos_half_angle = Float(0.9999891776066407); // cos(half_angle); - constexpr Float solid_angle = Float(6.799910294339209e-05); // 2.*M_PI*(1-cos_half_angle); - enum class Phase_kind {Lambertian, Specular, Rayleigh, HG, Mie}; template __device__ @@ -246,7 +242,7 @@ namespace inline Float probability_from_sun( Photon photon, const Vector& sun_direction, - const Float solid_angle, const Float g, + const Float sun_solid_angle, const Float g, const Float* __restrict__ mie_phase_ang, const Float* __restrict__ mie_phase, const Float r_eff, @@ -257,20 +253,20 @@ namespace const Float cos_angle = dot(photon.direction, sun_direction); if (kind == Phase_kind::HG) { - return henyey_phase(g, cos_angle) * solid_angle; + return henyey_phase(g, cos_angle) * sun_solid_angle; } 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) * solid_angle; - return mie_interpolate_phase_table(mie_phase_ang, mie_phase, acos(cos_angle), r_eff, mie_table_size) * solid_angle; + // 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; } else if (kind == Phase_kind::Rayleigh) { - return rayleigh_phase(cos_angle) * solid_angle; + return rayleigh_phase(cos_angle) * sun_solid_angle; } else if (kind == Phase_kind::Lambertian) { - return lambertian_phase() * solid_angle; + return lambertian_phase() * sun_solid_angle; } else if (kind == Phase_kind::Specular) { @@ -466,7 +462,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, 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, Float(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, @@ -554,7 +550,7 @@ void ray_tracer_kernel_bw( : Phase_kind::Lambertian; // SUN SCATTERING GOES HERE - const Float p_sun = probability_from_sun(photon, sun_direction, 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, Float(0.), 0, surface_normal, surface_kind); const Float trans_sun = transmission_direct_sun(photon,n,rng,sun_direction, k_null_grid,k_ext, @@ -693,7 +689,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, solid_angle, g, mie_phase_ang_shared, mie_phase, r_eff[ijk], mie_table_size, + const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, g, mie_phase_ang_shared, mie_phase, r_eff[ijk], mie_table_size, surface_normal, kind); const Float trans_sun = transmission_direct_sun(photon,n,rng,sun_direction, k_null_grid,k_ext, @@ -738,4 +734,117 @@ void ray_tracer_kernel_bw( } } } + __global__ + void accumulate_clouds_kernel( + const Float* __restrict__ lwp, + const Float* __restrict__ iwp, + const Float* __restrict__ tau_cloud, + const Vector grid_d, + const Vector grid_size, + const Vector grid_cells, + Float* __restrict__ liwp_cam, + Float* __restrict__ tauc_cam, + Float* __restrict__ dist_cam, + Float* __restrict__ zen_cam, + const Camera camera) + { + const int pix = blockDim.x * blockIdx.x + threadIdx.x; + const Float s_eps = max(max(grid_size.z, grid_size.x), grid_size.y) * Float_epsilon; + Vector direction; + Vector position; + + if (pix < camera.nx * camera.ny) + { + Float liwp_sum = 0; + Float tauc_sum = 0; + Float dist = 0; + bool reached_cloud = false; + + //const int pix = n * pixels_per_thread + ipix; + const Float i = (Float(pix % camera.nx) + Float(0.5))/ Float(camera.nx); + const Float j = (Float(pix / camera.nx) + Float(0.5))/ Float(camera.ny); + + position = camera.position + s_eps; + + if (camera.fisheye) + { + const Float photon_zenith = i * Float(.5) * M_PI / camera.f_zoom; + const Float photon_azimuth = j * Float(2.) * M_PI; + const Vector dir_tmp = {sin(photon_zenith) * sin(photon_azimuth), sin(photon_zenith) * cos(photon_azimuth), cos(photon_zenith)}; + + direction = {dot(camera.mx, dir_tmp), dot(camera.my, dir_tmp), dot(camera.mz, dir_tmp) * Float(-1)}; + } + else + { + direction = normalize(camera.cam_width * (Float(2.)*i-Float(1.0)) + camera.cam_height * (Float(2.)*j-Float(1.0)) + camera.cam_depth); + } + + // first bring photon to top of dynamical domain + if ((position.z >= (grid_size.z - s_eps)) && (direction.z < Float(0.))) + { + const Float s = abs((position.z - grid_size.z)/direction.z); + position = position + direction * s - s_eps; + + // Cyclic boundary condition in x. + position.x = fmod(position.x, grid_size.x); + if (position.x < Float(0.)) + position.x += grid_size.x; + + // Cyclic boundary condition in y. + position.y = fmod(position.y, grid_size.y); + if (position.y < Float(0.)) + position.y += grid_size.y; + } + + while ((position.z <= grid_size.z - s_eps) && (position.z >= s_eps)) + { + const int i = float_to_int(position.x, grid_d.x, grid_cells.x); + const int j = float_to_int(position.y, grid_d.y, grid_cells.y); + const int k = float_to_int(position.z, grid_d.z, grid_cells.z); + const int ijk = i + j*grid_cells.x + k*grid_cells.x*grid_cells.y; + + const Float sx = abs((direction.x > 0) ? ((i+1) * grid_d.x - position.x)/direction.x : (i*grid_d.x - position.x)/direction.x); + const Float sy = abs((direction.y > 0) ? ((j+1) * grid_d.y - position.y)/direction.y : (j*grid_d.y - position.y)/direction.y); + const Float sz = abs((direction.z > 0) ? ((k+1) * grid_d.z - position.z)/direction.z : (k*grid_d.z - position.z)/direction.z); + const Float s_min = min(sx, min(sy, sz)); + + liwp_sum += s_min * (lwp[ijk] + iwp[ijk]); + tauc_sum += s_min * tau_cloud[ijk]; + if (!reached_cloud) + { + dist += s_min; + reached_cloud = tau_cloud[ijk] > 0; + } + + position = position + direction * s_min; + + position.x += direction.x >= 0 ? s_eps : -s_eps; + position.y += direction.y >= 0 ? s_eps : -s_eps; + position.z += direction.z >= 0 ? s_eps : -s_eps; + + // Cyclic boundary condition in x. + position.x = fmod(position.x, grid_size.x); + if (position.x < Float(0.)) + position.x += grid_size.x; + + // Cyclic boundary condition in y. + position.y = fmod(position.y, grid_size.y); + if (position.y < Float(0.)) + position.y += grid_size.y; + + } + + // divide out initial layer thicknes, equivalent to first converting lwp (g/m2) to lwc (g/m3) or optical depth to k_ext(1/m) + liwp_cam[pix] = liwp_sum / grid_d.z; + tauc_cam[pix] = tauc_sum / grid_d.z; + dist_cam[pix] = reached_cloud ? dist : Float(-1) ; + zen_cam[pix] = acos(direction.z); + } + + + + + } + + diff --git a/src_test/Radiation_solver_bw.cu b/src_test/Radiation_solver_bw.cu index 00d4cd4..6d10d4d 100644 --- a/src_test/Radiation_solver_bw.cu +++ b/src_test/Radiation_solver_bw.cu @@ -733,31 +733,39 @@ Radiation_solver_shortwave::Radiation_solver_shortwave( } void Radiation_solver_shortwave::load_mie_tables( - const std::string& file_name_mie, - const bool switch_broadband) + const std::string& file_name_mie_bb, + const std::string& file_name_mie_im, + const bool switch_broadband, + const bool switch_image) { - Netcdf_file mie_nc(file_name_mie, Netcdf_mode::Read); - const int n_bnd_sw = this->get_n_bnd_gpu(); - const int n_re = mie_nc.get_dimension_size("r_eff"); - const int n_mie = mie_nc.get_dimension_size("n_ang"); - if (switch_broadband) { - Array mie_cdf(mie_nc.get_variable("phase_cdf", {n_bnd_sw, 1, n_mie}), {n_mie, 1, n_bnd_sw}); - Array mie_ang(mie_nc.get_variable("phase_cdf_angle", {n_bnd_sw, 1, n_re, n_mie}), {n_mie, n_re, 1, n_bnd_sw}); - - Array mie_phase(mie_nc.get_variable("phase", {n_bnd_sw, 1, n_re, n_mie}), {n_mie, n_re, 1, n_bnd_sw}); + Netcdf_file mie_nc(file_name_mie_bb, Netcdf_mode::Read); + const int n_bnd_sw = this->get_n_bnd_gpu(); + const int n_re = mie_nc.get_dimension_size("r_eff"); + const int n_mie = mie_nc.get_dimension_size("n_ang"); + + Array mie_cdf(mie_nc.get_variable("phase_cdf", {n_bnd_sw, n_mie}), {n_mie, 1, n_bnd_sw}); + Array mie_ang(mie_nc.get_variable("phase_cdf_angle", {n_bnd_sw, n_re, n_mie}), {n_mie, n_re, 1, n_bnd_sw}); + + Array mie_phase(mie_nc.get_variable("phase", {n_bnd_sw, n_re, n_mie}), {n_mie, n_re, 1, n_bnd_sw}); Array mie_phase_ang(mie_nc.get_variable("phase_angle", {n_mie}), {n_mie}); - this->mie_cdfs = mie_cdf; - this->mie_angs = mie_ang; + this->mie_cdfs_bb = mie_cdf; + this->mie_angs_bb = mie_ang; - this->mie_phase = mie_phase; - this->mie_phase_angs = mie_phase_ang; + this->mie_phase_bb = mie_phase; + this->mie_phase_angs_bb = mie_phase_ang; } - else + + if (switch_image) { + Netcdf_file mie_nc(file_name_mie_im, Netcdf_mode::Read); + const int n_bnd_sw = this->get_n_bnd_gpu(); + const int n_re = mie_nc.get_dimension_size("r_eff"); + const int n_mie = mie_nc.get_dimension_size("n_ang"); const int n_sub = mie_nc.get_dimension_size("sub_band"); + Array mie_cdf(mie_nc.get_variable("phase_cdf", {n_bnd_sw, n_sub, n_mie}), {n_mie, n_sub, n_bnd_sw}); Array mie_ang(mie_nc.get_variable("phase_cdf_angle", {n_bnd_sw, n_sub, n_re, n_mie}), {n_mie, n_re, n_sub, n_bnd_sw}); @@ -783,6 +791,8 @@ void Radiation_solver_shortwave::solve_gpu( const bool switch_lu_albedo, const bool switch_delta_cloud, const bool switch_delta_aerosol, + const bool switch_cloud_cam, + const bool switch_raytracing, const Vector& grid_cells, const Vector& grid_d, const Vector& kn_grid, @@ -801,7 +811,11 @@ void Radiation_solver_shortwave::solve_gpu( const Array_gpu& rh, const Aerosol_concs_gpu& aerosol_concs, const Camera& camera, - Array_gpu& XYZ) + Array_gpu& XYZ, + Array_gpu& liwp_cam, + Array_gpu& tauc_cam, + Array_gpu& dist_cam, + Array_gpu& zen_cam) { const int n_col = p_lay.dim(1); @@ -844,214 +858,240 @@ void Radiation_solver_shortwave::solve_gpu( const Array& band_limits_gpt(this->kdist_gpu->get_band_lims_gpoint()); - int previous_band = 0; - for (int igpt=1; igpt<=n_gpt; ++igpt) + if (switch_raytracing) { - int band = 0; - for (int ibnd=1; ibnd<=n_bnd; ++ibnd) + int previous_band = 0; + for (int igpt=1; igpt<=n_gpt; ++igpt) { - if (igpt <= band_limits_gpt({2, ibnd})) + int band = 0; + for (int ibnd=1; ibnd<=n_bnd; ++ibnd) { - band = ibnd; - break; + if (igpt <= band_limits_gpt({2, ibnd})) + { + band = ibnd; + break; + } } - } - 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; - - const Float solar_source_band = kdist_gpu->band_source(band_limits_gpt({1,band}), band_limits_gpt({2,band})); - - constexpr int n_col_block = 1<<14; // 2^14 - - Array_gpu toa_src_temp({n_col_block}); - auto gas_optics_subset = [&]( - const int col_s, const int n_col_subset) - { - // Run the gas_optics on a subset. - kdist_gpu->gas_optics( - igpt, - col_s, - n_col_subset, - n_col, - p_lay, - p_lev, - t_lay, - gas_concs, - optical_props, - toa_src_temp, - col_dry); - }; - - const int n_blocks = n_col / n_col_block; - const int n_col_residual = n_col % n_col_block; - - std::unique_ptr optical_props_block = - std::make_unique(n_col_block, n_lay, *kdist_gpu); - - for (int n=0; n 0) - { - const int col_s = n_blocks*n_col_block; - gas_optics_subset(col_s, n_col_residual); - } - - if (tune_step) return; - - toa_src.fill(toa_src_temp({1}) * tsi_scaling({1})); - if (switch_cloud_optics) - { - if (band > previous_band) + 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; + + const Float solar_source_band = kdist_gpu->band_source(band_limits_gpt({1,band}), band_limits_gpt({2,band})); + + constexpr int n_col_block = 1<<14; // 2^14 + + Array_gpu toa_src_temp({n_col_block}); + auto gas_optics_subset = [&]( + const int col_s, const int n_col_subset) { - cloud_optics_gpu->cloud_optics( - band, - lwp, - iwp, - rel, - rei, - *cloud_optical_props); - - if (switch_delta_cloud) - cloud_optical_props->delta_scale(); + // Run the gas_optics on a subset. + kdist_gpu->gas_optics( + igpt, + col_s, + n_col_subset, + n_col, + p_lay, + p_lev, + t_lay, + gas_concs, + optical_props, + toa_src_temp, + col_dry); + }; + + const int n_blocks = n_col / n_col_block; + const int n_col_residual = n_col % n_col_block; + + std::unique_ptr optical_props_block = + std::make_unique(n_col_block, n_lay, *kdist_gpu); + + for (int n=0; n(*optical_props), - dynamic_cast(*cloud_optical_props)); - } - else - { - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_tau().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_ssa().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_g().ptr()); - } - - if (switch_aerosol_optics) - { - if (band > previous_band) + + if (n_col_residual > 0) { - Aerosol_concs_gpu aerosol_concs_subset(aerosol_concs, 1, n_col); - aerosol_optics_gpu->aerosol_optics( - band, - aerosol_concs_subset, - rh, p_lev, - *aerosol_optical_props); - - if (switch_delta_aerosol) - aerosol_optical_props->delta_scale(); + const int col_s = n_blocks*n_col_block; + gas_optics_subset(col_s, n_col_residual); } - - // Add the aerosol optical props to the gas optical properties. - add_to( - dynamic_cast(*optical_props), - dynamic_cast(*aerosol_optical_props)); - } - else - { - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_tau().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_ssa().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_g().ptr()); - } - - - /* rrtmgp's bands are quite broad, we divide each spectral band in three equally broad spectral intervals - and run each g-point for each spectral interval, using the mean rayleigh scattering coefficient of each spectral interval - in stead of RRTMGP's rayleigh scattering coefficients. - The contribution of each spectral interval to the spectral band is based on the integrated (<>) Planck source function: - / , with a sun temperature of 5778 K. This is not entirely accurate because - the sun is not a black body radiatior, but the approximations comes close enough. - - */ - - // number of intervals - const Array& band_limits_wn(this->kdist_gpu->get_band_lims_wavenumber()); - const int nwv = n_sub; - - const Float wv1 = 1. / band_limits_wn({2,band}) * Float(1.e7); - const Float wv2 = 1. / band_limits_wn({1,band}) * Float(1.e7); - const Float dwv = (wv2-wv1)/Float(nwv); - - // - const Float total_planck = Planck_integrator(wv1,wv2); - - for (int iwv=0; iwv xyz_factor({3}); - xyz_factor({1}) = xyz_irradiance(wv1_sub,wv2_sub,&get_x); - xyz_factor({2}) = xyz_irradiance(wv1_sub,wv2_sub,&get_y); - xyz_factor({3}) = xyz_irradiance(wv1_sub,wv2_sub,&get_z); - Array_gpu xyz_factor_gpu(xyz_factor); - - if (switch_lu_albedo) + + if (tune_step) return; + + toa_src.fill(toa_src_temp({1}) * tsi_scaling({1})); + if (switch_cloud_optics) { - if (albedo.size() == 0) albedo.set_dims({1, n_col}); - spectral_albedo(n_col, wv1, wv2, land_use_map, albedo); + if (band > previous_band) + { + cloud_optics_gpu->cloud_optics( + band, + lwp, + iwp, + rel, + rei, + *cloud_optical_props); + + if (switch_delta_cloud) + cloud_optical_props->delta_scale(); + } + + // Add the cloud optical props to the gas optical properties. + add_to( + dynamic_cast(*optical_props), + dynamic_cast(*cloud_optical_props)); } - - const Float zenith_angle = std::acos(mu0({1})); - const Float azimuth_angle = azi({1}); - - if (switch_cloud_mie) + else { - mie_cdfs_sub = mie_cdfs.subset({{ {1, n_mie}, {iwv+1,iwv+1}, {band, band} }}); - mie_angs_sub = mie_angs.subset({{ {1, n_mie}, {1, n_re}, {iwv+1,iwv+1}, {band, band} }}); - mie_phase_sub = mie_phase.subset({{ {1, n_mie}, {1, n_re}, {iwv+1,iwv+1}, {band, band} }}); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_tau().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_ssa().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_g().ptr()); } - - raytracer.trace_rays( - igpt, - photons_per_pixel, n_lay, - grid_cells, grid_d, kn_grid, - z_lev, - mie_cdfs_sub, - mie_angs_sub, - mie_phase_sub, - mie_phase_angs, - rel, - 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(), - albedo, - land_use_map, - zenith_angle, - azimuth_angle, - toa_src({1}), - toa_factor, - rayleigh, - col_dry, - gas_concs.get_vmr("h2o"), - camera, - flux_camera); - - raytracer.add_xyz_camera( - camera, - xyz_factor_gpu, - flux_camera, - XYZ); - + + if (switch_aerosol_optics) + { + if (band > previous_band) + { + Aerosol_concs_gpu aerosol_concs_subset(aerosol_concs, 1, n_col); + aerosol_optics_gpu->aerosol_optics( + band, + aerosol_concs_subset, + rh, p_lev, + *aerosol_optical_props); + + if (switch_delta_aerosol) + aerosol_optical_props->delta_scale(); + } + + // Add the aerosol optical props to the gas optical properties. + add_to( + dynamic_cast(*optical_props), + dynamic_cast(*aerosol_optical_props)); + } + else + { + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_tau().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_ssa().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_g().ptr()); + } + + + /* rrtmgp's bands are quite broad, we divide each spectral band in three equally broad spectral intervals + and run each g-point for each spectral interval, using the mean rayleigh scattering coefficient of each spectral interval + in stead of RRTMGP's rayleigh scattering coefficients. + The contribution of each spectral interval to the spectral band is based on the integrated (<>) Planck source function: + / , with a sun temperature of 5778 K. This is not entirely accurate because + the sun is not a black body radiatior, but the approximations comes close enough. + + */ + + // number of intervals + const Array& band_limits_wn(this->kdist_gpu->get_band_lims_wavenumber()); + const int nwv = n_sub; + + const Float wv1 = 1. / band_limits_wn({2,band}) * Float(1.e7); + const Float wv2 = 1. / band_limits_wn({1,band}) * Float(1.e7); + const Float dwv = (wv2-wv1)/Float(nwv); + + // + const Float total_planck = Planck_integrator(wv1,wv2); + + for (int iwv=0; iwv xyz_factor({3}); + xyz_factor({1}) = xyz_irradiance(wv1_sub,wv2_sub,&get_x); + xyz_factor({2}) = xyz_irradiance(wv1_sub,wv2_sub,&get_y); + xyz_factor({3}) = xyz_irradiance(wv1_sub,wv2_sub,&get_z); + Array_gpu xyz_factor_gpu(xyz_factor); + + if (switch_lu_albedo) + { + if (albedo.size() == 0) albedo.set_dims({1, n_col}); + spectral_albedo(n_col, wv1, wv2, land_use_map, albedo); + } + + const Float zenith_angle = std::acos(mu0({1})); + const Float azimuth_angle = azi({1}); + + if (switch_cloud_mie) + { + mie_cdfs_sub = mie_cdfs.subset({{ {1, n_mie}, {iwv+1,iwv+1}, {band, band} }}); + mie_angs_sub = mie_angs.subset({{ {1, n_mie}, {1, n_re}, {iwv+1,iwv+1}, {band, band} }}); + mie_phase_sub = mie_phase.subset({{ {1, n_mie}, {1, n_re}, {iwv+1,iwv+1}, {band, band} }}); + } + + raytracer.trace_rays( + igpt, + photons_per_pixel, n_lay, + grid_cells, grid_d, kn_grid, + z_lev, + mie_cdfs_sub, + mie_angs_sub, + mie_phase_sub, + mie_phase_angs, + rel, + 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(), + albedo, + land_use_map, + zenith_angle, + azimuth_angle, + toa_src({1}), + toa_factor, + rayleigh, + col_dry, + gas_concs.get_vmr("h2o"), + camera, + flux_camera); + + raytracer.add_xyz_camera( + camera, + xyz_factor_gpu, + flux_camera, + XYZ); + + } + previous_band = band; } - previous_band = band; + } + + if (switch_cloud_cam) + { + cloud_optics_gpu->cloud_optics( + 11, // 441-615 nm band + lwp, + iwp, + rel, + rei, + *cloud_optical_props); + + raytracer.accumulate_clouds( + grid_d, + grid_cells, + lwp, + iwp, + dynamic_cast(*cloud_optical_props).get_tau(), + camera, + liwp_cam, + tauc_cam, + dist_cam, + zen_cam); } } @@ -1062,6 +1102,8 @@ void Radiation_solver_shortwave::solve_gpu_bb( const bool switch_lu_albedo, const bool switch_delta_cloud, const bool switch_delta_aerosol, + const bool switch_cloud_cam, + const bool switch_raytracing, const Vector& grid_cells, const Vector& grid_d, const Vector& kn_grid, @@ -1080,8 +1122,11 @@ void Radiation_solver_shortwave::solve_gpu_bb( const Array_gpu& rh, const Aerosol_concs_gpu& aerosol_concs, const Camera& camera, - Array_gpu& radiance) - + Array_gpu& radiance, + Array_gpu& liwp_cam, + Array_gpu& tauc_cam, + Array_gpu& dist_cam, + Array_gpu& zen_cam) { const int n_col = p_lay.dim(1); const int n_lay = p_lay.dim(2); @@ -1089,9 +1134,9 @@ void Radiation_solver_shortwave::solve_gpu_bb( const int n_gpt = this->kdist_gpu->get_ngpt(); const int n_bnd = this->kdist_gpu->get_nband(); - const int n_mie = (switch_cloud_mie) ? this->mie_angs.dim(1) : 0; - const int n_re = (switch_cloud_mie) ? this->mie_angs.dim(2) : 0; - + const int n_mie = (switch_cloud_mie) ? this->mie_angs_bb.dim(1) : 0; + const int n_re = (switch_cloud_mie) ? this->mie_angs_bb.dim(2) : 0; + const int cam_nx = radiance.dim(1); const int cam_ny = radiance.dim(2); const Bool top_at_1 = p_lay({1, 1}) < p_lay({1, n_lay}); @@ -1121,167 +1166,193 @@ void Radiation_solver_shortwave::solve_gpu_bb( const Array& band_limits_gpt(this->kdist_gpu->get_band_lims_gpoint()); int previous_band = 0; - for (int igpt=1; igpt<=n_gpt; ++igpt) + if (switch_raytracing) { - int band = 0; - for (int ibnd=1; ibnd<=n_bnd; ++ibnd) + for (int igpt=1; igpt<=n_gpt; ++igpt) { - if (igpt <= band_limits_gpt({2, ibnd})) + int band = 0; + for (int ibnd=1; ibnd<=n_bnd; ++ibnd) { - band = ibnd; - break; + if (igpt <= band_limits_gpt({2, ibnd})) + { + band = ibnd; + break; + } } - } - constexpr int n_col_block = 1<<14; // 2^14 + constexpr int n_col_block = 1<<14; // 2^14 - Array_gpu toa_src_temp({n_col_block}); + Array_gpu toa_src_temp({n_col_block}); - auto gas_optics_subset = [&]( - const int col_s, const int n_col_subset) - { - // Run the gas_optics on a subset. - kdist_gpu->gas_optics( - igpt, - col_s, - n_col_subset, - n_col, - p_lay, - p_lev, - t_lay, - gas_concs, - optical_props, - toa_src_temp, - col_dry); - }; + auto gas_optics_subset = [&]( + const int col_s, const int n_col_subset) + { + // Run the gas_optics on a subset. + kdist_gpu->gas_optics( + igpt, + col_s, + n_col_subset, + n_col, + p_lay, + p_lev, + t_lay, + gas_concs, + optical_props, + toa_src_temp, + col_dry); + }; + + const int n_blocks = n_col / n_col_block; + const int n_col_residual = n_col % n_col_block; + + for (int n=0; n 0) + { + const int col_s = n_blocks*n_col_block; + gas_optics_subset(col_s, n_col_residual); + } - for (int n=0; n 0) - { - const int col_s = n_blocks*n_col_block; - gas_optics_subset(col_s, n_col_residual); - } + if (switch_cloud_optics) + { + if (band > previous_band) + { + cloud_optics_gpu->cloud_optics( + band, + lwp, + iwp, + rel, + rei, + *cloud_optical_props); + + if (switch_delta_cloud) + cloud_optical_props->delta_scale(); + } + // Add the cloud optical props to the gas optical properties. + add_to( + dynamic_cast(*optical_props), + dynamic_cast(*cloud_optical_props)); + } + else + { + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_tau().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_ssa().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_g().ptr()); + } - toa_src.fill(toa_src_temp({1}) * tsi_scaling({1})); + if (switch_aerosol_optics) + { + if (band > previous_band) + { + Aerosol_concs_gpu aerosol_concs_subset(aerosol_concs, 1, n_col); + aerosol_optics_gpu->aerosol_optics( + band, + aerosol_concs_subset, + rh, p_lev, + *aerosol_optical_props); + + if (switch_delta_aerosol) + aerosol_optical_props->delta_scale(); + } + + // Add the aerosol optical props to the gas optical properties. + add_to( + dynamic_cast(*optical_props), + dynamic_cast(*aerosol_optical_props)); + } + else + { + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_tau().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_ssa().ptr()); + Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_g().ptr()); + } - if (switch_cloud_optics) - { - if (band > previous_band) + const Float zenith_angle = std::acos(mu0({1})); + const Float azimuth_angle = azi({1}); + + Array_gpu albedo; + if (switch_lu_albedo) { - cloud_optics_gpu->cloud_optics( - band, - lwp, - iwp, - rel, - rei, - *cloud_optical_props); + albedo.set_dims({1, n_col}); - if (switch_delta_cloud) - cloud_optical_props->delta_scale(); + const Array& band_limits_wn(this->kdist_gpu->get_band_lims_wavenumber()); + const Float wv1 = 1. / band_limits_wn({2,band}) * Float(1.e7); + const Float wv2 = 1. / band_limits_wn({1,band}) * Float(1.e7); + spectral_albedo(n_col, wv1, wv2, land_use_map, albedo); } - // Add the cloud optical props to the gas optical properties. - add_to( - dynamic_cast(*optical_props), - dynamic_cast(*cloud_optical_props)); - } - else - { - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_tau().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_ssa().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_g().ptr()); - } - - if (switch_aerosol_optics) - { - if (band > previous_band) + else { - Aerosol_concs_gpu aerosol_concs_subset(aerosol_concs, 1, n_col); - aerosol_optics_gpu->aerosol_optics( - band, - aerosol_concs_subset, - rh, p_lev, - *aerosol_optical_props); - - if (switch_delta_aerosol) - aerosol_optical_props->delta_scale(); + albedo = sfc_alb.subset({{ {band, band}, {1, n_col}}}); } - // Add the aerosol optical props to the gas optical properties. - add_to( - dynamic_cast(*optical_props), - dynamic_cast(*aerosol_optical_props)); - } - else - { - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_tau().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_ssa().ptr()); - Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_g().ptr()); - } - - const Float zenith_angle = std::acos(mu0({1})); - const Float azimuth_angle = azi({1}); + if (switch_cloud_mie) + { + mie_cdfs_sub = mie_cdfs_bb.subset({{ {1, n_mie}, {1,1}, {band, band} }}); + mie_angs_sub = mie_angs_bb.subset({{ {1, n_mie}, {1, n_re}, {1,1}, {band, band} }}); + mie_phase_sub = mie_phase_bb.subset({{ {1, n_mie}, {1, n_re}, {1,1}, {band, band} }}); + } - Array_gpu albedo; - if (switch_lu_albedo) - { - albedo.set_dims({1, n_col}); + raytracer.trace_rays_bb( + igpt, + photons_per_pixel, n_lay, + grid_cells, grid_d, kn_grid, + z_lev, + mie_cdfs_sub, + mie_angs_sub, + mie_phase_sub, + mie_phase_angs_bb, + rel, + 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(), + albedo, + land_use_map, + zenith_angle, + azimuth_angle, + toa_src({1}), + camera, + flux_camera); - const Array& band_limits_wn(this->kdist_gpu->get_band_lims_wavenumber()); - const Float wv1 = 1. / band_limits_wn({2,band}) * Float(1.e7); - const Float wv2 = 1. / band_limits_wn({1,band}) * Float(1.e7); - spectral_albedo(n_col, wv1, wv2, land_use_map, albedo); - } - else - { - albedo = sfc_alb.subset({{ {band, band}, {1, n_col}}}); - } + raytracer.add_camera( + camera, + flux_camera, + radiance); - if (switch_cloud_mie) - { - mie_cdfs_sub = mie_cdfs.subset({{ {1, n_mie}, {1,1}, {band, band} }}); - mie_angs_sub = mie_angs.subset({{ {1, n_mie}, {1, n_re}, {1,1}, {band, band} }}); - mie_phase_sub = mie_phase.subset({{ {1, n_mie}, {1, n_re}, {1,1}, {band, band} }}); + previous_band = band; } + } - raytracer.trace_rays_bb( - igpt, - photons_per_pixel, n_lay, - grid_cells, grid_d, kn_grid, - z_lev, - mie_cdfs_sub, - mie_angs_sub, - mie_phase_sub, - mie_phase_angs, + if (switch_cloud_cam) + { + cloud_optics_gpu->cloud_optics( + 11, // 441-615 nm band + lwp, + iwp, rel, - 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(), - albedo, - land_use_map, - zenith_angle, - azimuth_angle, - toa_src({1}), + rei, + *cloud_optical_props); + + raytracer.accumulate_clouds( + grid_d, + grid_cells, + lwp, + iwp, + dynamic_cast(*cloud_optical_props).get_tau(), camera, - flux_camera); - - raytracer.add_camera( - camera, - flux_camera, - radiance); - - previous_band = band; + liwp_cam, + tauc_cam, + dist_cam, + zen_cam); } } diff --git a/src_test/Radiation_solver_rt.cu b/src_test/Radiation_solver_rt.cu index 259b498..90b4c8b 100644 --- a/src_test/Radiation_solver_rt.cu +++ b/src_test/Radiation_solver_rt.cu @@ -42,48 +42,6 @@ namespace { - __global__ - void compute_tod_flux_kernel( - const int ncol, const int nlay, const int col_per_thread, const Float* __restrict__ flux_dn, const Float* __restrict__ flux_dn_dir, Float* __restrict__ tod_dir_diff) - { - const int i = blockIdx.x*blockDim.x + threadIdx.x; - Float flx_dir = 0; - Float flx_tot = 0; - for (int icol = i*col_per_thread; icol < (i+1)*col_per_thread; ++icol) - { - if ( ( icol < ncol) ) - { - const int idx = icol + nlay*ncol; - flx_dir += flux_dn_dir[idx]; - flx_tot += flux_dn[idx]; - } - } - atomicAdd(&tod_dir_diff[0], flx_dir); - atomicAdd(&tod_dir_diff[1], flx_tot - flx_dir); - } - - void compute_tod_flux( - const int ncol, const int nlay, const Array_gpu& flux_dn, const Array_gpu& flux_dn_dir, Array& tod_dir_diff) - { - const int col_per_thread = 32; - const int nthread = int(ncol/col_per_thread) + 1; - const int block_col = 16; - const int grid_col = nthread/block_col + (nthread%block_col > 0); - - dim3 grid_gpu(grid_col, 1); - dim3 block_gpu(block_col, 1); - - Array_gpu tod_dir_diff_g({2}); - tod_dir_diff_g.fill(Float(0.)); - - compute_tod_flux_kernel<<>>( - ncol, nlay, col_per_thread, flux_dn.ptr(), flux_dn_dir.ptr(), tod_dir_diff_g.ptr()); - Array tod_dir_diff_c(tod_dir_diff_g); - - tod_dir_diff({1}) = tod_dir_diff_c({1}) / Float(ncol); - tod_dir_diff({2}) = tod_dir_diff_c({2}) / Float(ncol); - } - std::vector get_variable_string( const std::string& var_name, std::vector i_count, @@ -604,13 +562,13 @@ void Radiation_solver_shortwave::load_mie_tables( void Radiation_solver_shortwave::solve_gpu( const bool switch_fluxes, const bool switch_raytracing, + const bool switch_independent_column, const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, const bool switch_single_gpt, const bool switch_delta_cloud, const bool switch_delta_aerosol, - const bool switch_clear_sky_tod, const int single_gpt, const Int ray_count, const Vector grid_cells, @@ -770,29 +728,6 @@ void Radiation_solver_shortwave::solve_gpu( Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_g().ptr()); } - Array tod_dir_diff({2}); - - if (switch_fluxes && switch_raytracing && switch_clear_sky_tod) - { - // If switch_clear_sky_tod, compute TOD irradiance for ray tracer with clear-sky atmosphere, otherwise compute TOD irradiance after cloud optics - std::unique_ptr fluxes = - std::make_unique(grid_cells.x, grid_cells.y, n_lev); - - 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()); - - compute_tod_flux(n_col, grid_cells.z, (*fluxes).get_flux_dn(), (*fluxes).get_flux_dn_dir(), tod_dir_diff); - } - if (switch_cloud_optics) { if (band > previous_band) @@ -855,9 +790,6 @@ void Radiation_solver_shortwave::solve_gpu( Float zenith_angle = std::acos(mu0({1})); Float azimuth_angle = azi({1}); // sun approximately from south - if (!switch_clear_sky_tod) - compute_tod_flux(n_col, grid_cells.z, (*fluxes).get_flux_dn(), (*fluxes).get_flux_dn_dir(), tod_dir_diff); - if (switch_cloud_mie) { mie_cdfs_sub = mie_cdfs.subset({{ {1, n_mie}, {band, band} }}); @@ -866,6 +798,7 @@ void Radiation_solver_shortwave::solve_gpu( raytracer.trace_rays( igpt, + switch_independent_column, ray_count, grid_cells, grid_d, @@ -884,8 +817,8 @@ void Radiation_solver_shortwave::solve_gpu( sfc_alb_dir.subset({{ {band, band}, {1, n_col}}}), zenith_angle, azimuth_angle, - tod_dir_diff({1}), - tod_dir_diff({2}), + toa_src({1}) * mu0({1}), + Float(0.), (*fluxes).get_flux_tod_dn(), (*fluxes).get_flux_tod_up(), (*fluxes).get_flux_sfc_dir(), diff --git a/src_test/test_rt_lite.cu b/src_test/test_rt_lite.cu index df912bf..66691a4 100644 --- a/src_test/test_rt_lite.cu +++ b/src_test/test_rt_lite.cu @@ -126,8 +126,9 @@ void solve_radiation(int argc, char** argv) ////// FLOW CONTROL SWITCHES ////// // Parse the command line options. std::map> command_line_switches { - {"raytracing" , { true, "Use raytracing for flux computation. '--raytracing 256': use 256 rays per pixel" }}, - {"profiling" , { false, "Perform additional profiling run." }} }; + {"raytracing" , { true, "Use raytracing for flux computation. '--raytracing 256': use 256 rays per pixel" }}, + {"independent-column", { false, "run raytracer in independent column mode"}}, + {"profiling" , { false, "Perform additional profiling run." }} }; std::map> command_line_ints { {"raytracing", {32, "Number of rays initialised at TOD per pixel per quadraute."}}} ; @@ -135,9 +136,10 @@ void solve_radiation(int argc, char** argv) if (parse_command_line_options(command_line_switches, command_line_ints, argc, argv)) return; - const bool switch_raytracing = command_line_switches.at("raytracing" ).first; - const bool switch_profiling = command_line_switches.at("profiling" ).first; - + const bool switch_raytracing = command_line_switches.at("raytracing" ).first; + const bool switch_independent_column = command_line_switches.at("independent-column").first; + const bool switch_profiling = command_line_switches.at("profiling" ).first; + // Print the options to the screen. print_command_line_options(command_line_switches, command_line_ints); @@ -257,6 +259,7 @@ void solve_radiation(int argc, char** argv) raytracer.trace_rays( 0, + switch_independent_column, photons_per_pixel, grid_cells, grid_d, diff --git a/src_test/test_rte_rrtmgp_bw.cu b/src_test/test_rte_rrtmgp_bw.cu index 24cc880..4929e12 100644 --- a/src_test/test_rte_rrtmgp_bw.cu +++ b/src_test/test_rte_rrtmgp_bw.cu @@ -238,17 +238,19 @@ void solve_radiation(int argc, char** argv) {"shortwave" , { true, "Enable computation of shortwave radiation." }}, {"longwave" , { false, "Enable computation of longwave radiation." }}, {"fluxes" , { true, "Enable computation of fluxes." }}, - {"raytracing" , { false, "Use raytracing for flux computation." }}, + {"raytracing" , { true, "Use raytracing for flux computation." }}, {"cloud-optics" , { false, "Enable cloud optics." }}, {"cloud-mie" , { false, "mie cloud droplet scattering." }}, {"aerosol-optics" , { false, "Enable aerosol optics." }}, {"output-optical" , { false, "Enable output of optical properties." }}, {"output-bnd-fluxes", { false, "Enable output of band fluxes." }}, {"lu-albedo" , { false, "Compute spectral albedo from land use map" }}, - {"broadband" , { false, "Compute broadband fluxes" }}, + {"image" , { true, "Compute XYZ values to generate RGB images" }}, + {"broadband" , { false, "Compute broadband radiances" }}, {"profiling" , { false, "Perform additional profiling run." }}, {"delta-cloud" , { false, "delta-scaling of cloud optical properties" }}, - {"delta-aerosol" , { false, "delta-scaling of aerosol optical properties" }}}; + {"delta-aerosol" , { false, "delta-scaling of aerosol optical properties" }}, + {"cloud-cam" , { false, "accumulate cloud water & ice paths for each camera pixel" }}}; int photons_per_pixel = 1; if (parse_command_line_options(command_line_options, photons_per_pixel, argc, argv)) @@ -264,10 +266,13 @@ void solve_radiation(int argc, char** argv) const bool switch_output_optical = command_line_options.at("output-optical" ).first; const bool switch_output_bnd_fluxes = command_line_options.at("output-bnd-fluxes").first; const bool switch_lu_albedo = command_line_options.at("lu-albedo" ).first; + const bool switch_image = command_line_options.at("image" ).first; const bool switch_broadband = command_line_options.at("broadband" ).first; const bool switch_profiling = command_line_options.at("profiling" ).first; const bool switch_delta_cloud = command_line_options.at("delta-cloud" ).first; const bool switch_delta_aerosol = command_line_options.at("delta-aerosol" ).first; + const bool switch_cloud_cam = command_line_options.at("cloud-cam" ).first; + const bool switch_raytracing = command_line_options.at("raytracing" ).first; if (switch_longwave) { @@ -378,14 +383,14 @@ void solve_radiation(int argc, char** argv) Array rel; Array rei; - if (switch_cloud_optics) + if (switch_cloud_optics || switch_cloud_cam) { lwp.set_dims({n_col, n_lay}); lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); iwp.set_dims({n_col, n_lay}); iwp = std::move(input_nc.get_variable("iwp", {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})); @@ -397,6 +402,7 @@ void solve_radiation(int argc, char** argv) rel.set_dims({n_col, n_lay}); rel.fill(Float(0.)); } + Array rh; Aerosol_concs aerosol_concs; @@ -418,6 +424,8 @@ void solve_radiation(int argc, char** argv) read_and_set_aer("aermr11", n_col_x, n_col_y, n_lay, input_nc, aerosol_concs); } + + ////// CREATE THE OUTPUT FILE ////// // Create the general dimensions and arrays. Status::print_message("Preparing NetCDF output file."); @@ -676,17 +684,29 @@ void solve_radiation(int argc, char** argv) if (switch_broadband) { radiance.set_dims({camera.nx, camera.ny}); - - if (switch_cloud_mie) - rad_sw.load_mie_tables("mie_lut_broadband.nc", switch_broadband); } - else + if (switch_image) { XYZ.set_dims({camera.nx, camera.ny, 3}); - - if (switch_cloud_mie) - rad_sw.load_mie_tables("mie_lut_visualisation.nc", switch_broadband); } + + if (switch_cloud_mie) + rad_sw.load_mie_tables("mie_lut_broadband.nc", "mie_lut_visualisation.nc", switch_broadband, switch_image); + + + Array_gpu liwp_cam; + Array_gpu tauc_cam; + Array_gpu dist_cam; + 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}); + } + // Solve the radiation. Status::print_message("Solving the shortwave radiation."); @@ -727,6 +747,8 @@ void solve_radiation(int argc, char** argv) switch_lu_albedo, switch_delta_cloud, switch_delta_aerosol, + switch_cloud_cam, + switch_raytracing, grid_cells, grid_d, kn_grid, @@ -745,7 +767,11 @@ void solve_radiation(int argc, char** argv) rh_gpu, aerosol_concs, camera, - radiance); + radiance, + liwp_cam, + tauc_cam, + dist_cam, + zen_cam); cudaEventRecord(stop, 0); cudaEventSynchronize(stop); @@ -755,7 +781,7 @@ void solve_radiation(int argc, char** argv) cudaEventDestroy(start); cudaEventDestroy(stop); - Status::print_message("Duration shortwave solver: " + std::to_string(duration) + " (ms)"); + Status::print_message("Duration shortwave solver (broadband version): " + std::to_string(duration) + " (ms)"); }; auto run_solver = [&](const bool tune_step) @@ -796,6 +822,8 @@ void solve_radiation(int argc, char** argv) switch_lu_albedo, switch_delta_cloud, switch_delta_aerosol, + switch_cloud_cam, + switch_raytracing, grid_cells, grid_d, kn_grid, @@ -814,7 +842,11 @@ void solve_radiation(int argc, char** argv) rh_gpu, aerosol_concs, camera, - XYZ); + XYZ, + liwp_cam, + tauc_cam, + dist_cam, + zen_cam); cudaEventRecord(stop, 0); cudaEventSynchronize(stop); @@ -824,7 +856,7 @@ void solve_radiation(int argc, char** argv) cudaEventDestroy(start); cudaEventDestroy(stop); - Status::print_message("Duration shortwave solver: " + std::to_string(duration) + " (ms)"); + Status::print_message("Duration shortwave solver (image version): " + std::to_string(duration) + " (ms)"); }; if (switch_broadband) @@ -839,7 +871,7 @@ void solve_radiation(int argc, char** argv) cudaProfilerStop(); } } - else + if (switch_image) { // tune step run_solver(true); @@ -859,35 +891,59 @@ void solve_radiation(int argc, char** argv) // Store the output. Status::print_message("Storing the shortwave output."); - if (switch_broadband) + if (switch_raytracing) { - Array radiance_cpu(radiance); output_nc.add_dimension("gpt_sw", n_gpt_sw); output_nc.add_dimension("band_sw", n_bnd_sw); auto nc_sw_band_lims_wvn = output_nc.add_variable("sw_band_lims_wvn", {"band_sw", "pair"}); nc_sw_band_lims_wvn.insert(rad_sw.get_band_lims_wavenumber_gpu().v(), {0, 0}); + + if (switch_broadband) + { + Array radiance_cpu(radiance); - auto nc_var = output_nc.add_variable("radiance", {"y","x"}); - nc_var.insert(radiance_cpu.v(), {0, 0}); - nc_var.add_attribute("long_name", "shortwave radiance"); - nc_var.add_attribute("units", "W m-2 sr-1"); - } - else - { - Array xyz_cpu(XYZ); - output_nc.add_dimension("gpt_sw", n_gpt_sw); - output_nc.add_dimension("band_sw", n_bnd_sw); - output_nc.add_dimension("n",3); + auto nc_var = output_nc.add_variable("radiance", {"y","x"}); + nc_var.insert(radiance_cpu.v(), {0, 0}); + nc_var.add_attribute("long_name", "shortwave radiance"); + nc_var.add_attribute("units", "W m-2 sr-1"); + } + if (switch_image) + { + Array xyz_cpu(XYZ); + output_nc.add_dimension("n",3); - auto nc_sw_band_lims_wvn = output_nc.add_variable("sw_band_lims_wvn", {"band_sw", "pair"}); - nc_sw_band_lims_wvn.insert(rad_sw.get_band_lims_wavenumber_gpu().v(), {0, 0}); + auto nc_xyz = output_nc.add_variable("XYZ", {"n","y","x"}); + nc_xyz.insert(xyz_cpu.v(), {0, 0, 0}); + + nc_xyz.add_attribute("long_name", "X Y Z tristimulus values"); + } + } - auto nc_xyz = output_nc.add_variable("XYZ", {"n","y","x"}); - nc_xyz.insert(xyz_cpu.v(), {0, 0, 0}); + 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", {"y","x"}); + nc_var_liwp.insert(liwp_cam_cpu.v(), {0, 0}); + nc_var_liwp.add_attribute("long_name", "accumulated liquid+ice water path"); - nc_xyz.add_attribute("long_name", "X Y Z tristimulus values"); + auto nc_var_tauc = output_nc.add_variable("tau_cld_cam", {"y","x"}); + 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", {"y","x"}); + 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", {"y","x"}); + 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}); diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index 005e082..cfa2062 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -224,18 +224,18 @@ void solve_radiation(int argc, char** argv) ////// FLOW CONTROL SWITCHES ////// // Parse the command line options. std::map> command_line_switches { - {"shortwave" , { true, "Enable computation of shortwave radiation."}}, - {"longwave" , { false, "Enable computation of longwave radiation." }}, - {"fluxes" , { true, "Enable computation of fluxes." }}, - {"raytracing" , { true, "Use raytracing for flux computation. '--raytracing 256': use 256 rays per pixel" }}, - {"cloud-optics" , { false, "Enable cloud optics." }}, - {"cloud-mie" , { false, "Use Mie tables for cloud scattering in ray tracer" }}, - {"aerosol-optics" , { false, "Enable aerosol optics." }}, - {"single-gpt" , { false, "Output optical properties and fluxes for a single g-point. '--single-gpt 100': output 100th g-point" }}, - {"profiling" , { false, "Perform additional profiling run." }}, - {"delta-cloud" , { false, "delta-scaling of cloud optical properties" }}, - {"delta-aerosol" , { false, "delta-scaling of aerosol optical properties" }}, - {"clear-sky-tod" , { false, "compute ray tracing TOD irradiance from clear-sky (no clouds) 2-stream computations" }} }; + {"shortwave" , { true, "Enable computation of shortwave radiation."}}, + {"longwave" , { false, "Enable computation of longwave radiation." }}, + {"fluxes" , { true, "Enable computation of fluxes." }}, + {"raytracing" , { true, "Use raytracing for flux computation. '--raytracing 256': use 256 rays per pixel" }}, + {"independent-column", { false, "run raytracer in independent column mode"}}, + {"cloud-optics" , { false, "Enable cloud optics." }}, + {"cloud-mie" , { false, "Use Mie tables for cloud scattering in ray tracer" }}, + {"aerosol-optics" , { false, "Enable aerosol optics." }}, + {"single-gpt" , { false, "Output optical properties and fluxes for a single g-point. '--single-gpt 100': output 100th g-point" }}, + {"profiling" , { false, "Perform additional profiling run." }}, + {"delta-cloud" , { false, "delta-scaling of cloud optical properties" }}, + {"delta-aerosol" , { false, "delta-scaling of aerosol optical properties" }} }; std::map> command_line_ints { {"raytracing", {32, "Number of rays initialised at TOD per pixel per quadraute."}}, @@ -244,18 +244,18 @@ void solve_radiation(int argc, char** argv) if (parse_command_line_options(command_line_switches, command_line_ints, argc, argv)) return; - const bool switch_shortwave = command_line_switches.at("shortwave" ).first; - const bool switch_longwave = command_line_switches.at("longwave" ).first; - const bool switch_fluxes = command_line_switches.at("fluxes" ).first; - const bool switch_raytracing = command_line_switches.at("raytracing" ).first; - const bool switch_cloud_optics = command_line_switches.at("cloud-optics" ).first; - const bool switch_cloud_mie = command_line_switches.at("cloud-mie" ).first; - const bool switch_aerosol_optics = command_line_switches.at("aerosol-optics" ).first; - const bool switch_single_gpt = command_line_switches.at("single-gpt" ).first; - const bool switch_profiling = command_line_switches.at("profiling" ).first; - const bool switch_delta_cloud = command_line_switches.at("delta-cloud" ).first; - const bool switch_delta_aerosol = command_line_switches.at("delta-aerosol" ).first; - const bool switch_clear_sky_tod = command_line_switches.at("clear-sky-tod" ).first; + const bool switch_shortwave = command_line_switches.at("shortwave" ).first; + const bool switch_longwave = command_line_switches.at("longwave" ).first; + const bool switch_fluxes = command_line_switches.at("fluxes" ).first; + const bool switch_raytracing = command_line_switches.at("raytracing" ).first; + const bool switch_independent_column= command_line_switches.at("independent-column").first; + const bool switch_cloud_optics = command_line_switches.at("cloud-optics" ).first; + const bool switch_cloud_mie = command_line_switches.at("cloud-mie" ).first; + const bool switch_aerosol_optics = command_line_switches.at("aerosol-optics" ).first; + const bool switch_single_gpt = command_line_switches.at("single-gpt" ).first; + const bool switch_profiling = command_line_switches.at("profiling" ).first; + const bool switch_delta_cloud = command_line_switches.at("delta-cloud" ).first; + const bool switch_delta_aerosol = command_line_switches.at("delta-aerosol" ).first; Int photons_per_pixel = Int(command_line_ints.at("raytracing").first); @@ -288,14 +288,18 @@ void solve_radiation(int argc, char** argv) const int n_col = n_col_x * n_col_y; const int n_lay = input_nc.get_dimension_size("lay"); const int n_lev = input_nc.get_dimension_size("lev"); - const int n_z = input_nc.get_dimension_size("z"); + + const int n_z_in = input_nc.get_dimension_size("z"); + + // number of vertical levels in the raytrace grid. We add 1 layer on top, in which we add integrated optical properties between TOD and TOA + const int n_z = n_z_in + 1; Array grid_x(input_nc.get_variable("x", {n_col_x}), {n_col_x}); Array grid_xh(input_nc.get_variable("xh", {n_col_x+1}), {n_col_x+1}); 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}), {n_z}); - + Array grid_z(input_nc.get_variable("z", {n_z_in}), {n_z_in}); + const Vector grid_cells = {n_col_x, n_col_y, n_z}; 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"), @@ -394,7 +398,7 @@ void solve_radiation(int argc, char** argv) output_nc.add_dimension("lev", n_lev); output_nc.add_dimension("pair", 2); - output_nc.add_dimension("z", n_z); + output_nc.add_dimension("z", n_z_in); auto nc_x = output_nc.add_variable("x", {"x"}); auto nc_y = output_nc.add_variable("y", {"y"}); auto nc_z = output_nc.add_variable("z", {"z"}); @@ -760,13 +764,13 @@ void solve_radiation(int argc, char** argv) rad_sw.solve_gpu( switch_fluxes, switch_raytracing, + switch_independent_column, switch_cloud_optics, switch_cloud_mie, switch_aerosol_optics, switch_single_gpt, switch_delta_cloud, switch_delta_aerosol, - switch_clear_sky_tod, single_gpt, photons_per_pixel, grid_cells, @@ -931,9 +935,10 @@ void solve_radiation(int argc, char** argv) nc_rt_flux_sfc_dir.insert(rt_flux_sfc_dir_cpu.v(), {0,0}); nc_rt_flux_sfc_dif.insert(rt_flux_sfc_dif_cpu.v(), {0,0}); nc_rt_flux_sfc_up .insert(rt_flux_sfc_up_cpu .v(), {0,0}); + nc_rt_flux_abs_dir.insert(rt_flux_abs_dir_cpu.v(), {0,0,0}); nc_rt_flux_abs_dif.insert(rt_flux_abs_dif_cpu.v(), {0,0,0}); - + 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");