Skip to content

Commit 103d7a4

Browse files
MennoVeermanMennoVeerman
authored
Merge pull request #51 from magpowell/tilted_column3
Tilted column3
2 parents af2aed8 + 23327ec commit 103d7a4

File tree

7 files changed

+1613
-8
lines changed

7 files changed

+1613
-8
lines changed

‎include/Array.h‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -172,6 +172,22 @@ class Array
172172
offsets = {};
173173
}
174174

175+
inline void expand_dims(const std::array<int, N>& dims)
176+
{
177+
if (this->ncells == 0)
178+
throw std::runtime_error("Only arrays of size > 0 can be expanded");
179+
180+
this->dims = dims;
181+
ncells = product<N>(dims);
182+
183+
if (data.size() > ncells)
184+
throw std::runtime_error("new dimensions too small for data");
185+
186+
data.resize(ncells);
187+
strides = calc_strides<N>(dims);
188+
offsets = {};
189+
}
190+
175191
inline std::vector<T>& v() { return data; }
176192
inline const std::vector<T>& v() const { return data; }
177193

‎include_test/Radiation_solver_rt.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -108,6 +108,7 @@ class Radiation_solver_shortwave
108108
const bool switch_single_gpt,
109109
const bool switch_delta_cloud,
110110
const bool switch_delta_aerosol,
111+
const bool switch_attenuate_tica,
111112
const int single_gpt,
112113
const Int ray_count,
113114
const Vector<int> grid_cells,

‎include_tilt/tilt_utils.h‎

Lines changed: 132 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,132 @@
1+
#include <cmath>
2+
#include <chrono>
3+
#include <iomanip>
4+
#include <iostream>
5+
#include <numeric>
6+
#include <random>
7+
#include <vector>
8+
9+
#include "Array.h"
10+
#include "Source_functions.h"
11+
12+
struct ColumnResult {
13+
Array<Float,1> lwp;
14+
Array<Float,1> iwp;
15+
Array<Float,1> rel;
16+
Array<Float,1> dei;
17+
};
18+
19+
struct ijk
20+
{
21+
int i;
22+
int j;
23+
int k;
24+
};
25+
26+
inline int sign(Float value)
27+
{
28+
return (Float(0.) < value) - (value < Float(0.));
29+
}
30+
31+
void tilted_path(std::vector<Float>& xh, std::vector<Float>& yh,
32+
std::vector<Float>& zh, std::vector<Float>& z,
33+
Float sza, Float azi,
34+
Float x_start, Float y_start,
35+
std::vector<ijk>& tilted_path,
36+
std::vector<Float>& zh_tilted);
37+
38+
void post_process_output(const std::vector<ColumnResult>& col_results,
39+
const int n_col_x, const int n_col_y,
40+
const int n_z, const int n_zh,
41+
Array<Float,2>* lwp_out,
42+
Array<Float,2>* iwp_out,
43+
Array<Float,2>* rel_out,
44+
Array<Float,2>* dei_out,
45+
const bool switch_liq_cloud_optics,
46+
const bool switch_ice_cloud_optics);
47+
48+
bool prepare_netcdf(Netcdf_handle& input_nc, std::string file_name, int n_lay, int n_lev, int n_col_x, int n_col_y,
49+
int n_zh, int n_z,
50+
Float sza, std::vector<Float> zh, std::vector<Float> z,
51+
Array<Float,2> p_lay, Array<Float,2> t_lay, Array<Float,2> p_lev, Array<Float,2> t_lev,
52+
Array<Float,2> lwp, Array<Float,2> iwp, Array<Float,2> rel, Array<Float,2> dei,
53+
Gas_concs& gas_concs, std::vector<std::string> gas_names,
54+
bool switch_cloud_optics, bool switch_liq_cloud_optics, bool switch_ice_cloud_optics);
55+
56+
void restore_bkg_profile(const int n_x, const int n_y,
57+
const int n_full,
58+
const int n_tilt,
59+
const int bkg_start,
60+
std::vector<Float>& var,
61+
std::vector<Float>& var_w_bkg);
62+
63+
void restore_bkg_profile_bundle(const int n_col_x, const int n_col_y,
64+
const int n_lay, const int n_lev,
65+
const int n_lay_tot, const int n_lev_tot,
66+
const int n_z_in, const int n_zh_in,
67+
const int bkg_start_z, const int bkg_start_zh,
68+
Array<Float,2>* p_lay_copy, Array<Float,2>* t_lay_copy, Array<Float,2>* p_lev_copy, Array<Float,2>* t_lev_copy,
69+
Array<Float,2>* lwp_copy, Array<Float,2>* iwp_copy, Array<Float,2>* rel_copy, Array<Float,2>* dei_copy, Array<Float,2>* rh_copy,
70+
Gas_concs& gas_concs_copy, Aerosol_concs& aerosol_concs_copy,
71+
Array<Float,2>* p_lay, Array<Float,2>* t_lay, Array<Float,2>* p_lev, Array<Float,2>* t_lev,
72+
Array<Float,2>* lwp, Array<Float,2>* iwp, Array<Float,2>* rel, Array<Float,2>* dei, Array<Float,2>* rh,
73+
Gas_concs& gas_concs, Aerosol_concs& aerosol_concs,
74+
std::vector<std::string> gas_names, std::vector<std::string> aerosol_names,
75+
bool switch_liq_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics
76+
);
77+
78+
void compress_columns_weighted_avg(const int n_x, const int n_y,
79+
const int n_out,
80+
const int n_tilt,
81+
const int compress_lay_start_idx,
82+
std::vector<Float>& var, std::vector<Float>& var_weighting);
83+
84+
void compress_columns_p_or_t(const int n_x, const int n_y,
85+
const int n_out_lay, const int n_tilt,
86+
const int compress_lay_start_idx,
87+
std::vector<Float>& var_lev, std::vector<Float>& var_lay);
88+
89+
void tilt_fields(const int n_z_in, const int n_zh_in, const int n_col_x, const int n_col_y,
90+
const int n_z_tilt, const int n_zh_tilt, const int n_col,
91+
const Array<Float,1> zh, const Array<Float,1> z,
92+
const Array<Float,1> zh_tilt, const Array<ijk,1> path,
93+
Array<Float,2>* p_lay_copy, Array<Float,2>* t_lay_copy, Array<Float,2>* p_lev_copy, Array<Float,2>* t_lev_copy,
94+
Array<Float,2>* rh_copy,
95+
Gas_concs& gas_concs_copy, const std::vector<std::string> gas_names,
96+
Aerosol_concs& aerosol_concs_copy, const std::vector<std::string> aerosol_names, const bool switch_aerosol_optics
97+
);
98+
99+
void compress_fields(const int compress_lay_start_idx, const int n_col_x, const int n_col_y,
100+
const int n_z_in, const int n_zh_in, const int n_z_tilt,
101+
Array<Float,2>* p_lay_copy, Array<Float,2>* t_lay_copy, Array<Float,2>* p_lev_copy, Array<Float,2>* t_lev_copy,
102+
Array<Float,2>* rh_copy,
103+
Gas_concs& gas_concs_copy, std::vector<std::string> gas_names,
104+
Aerosol_concs& aerosol_concs_copy, std::vector<std::string> aerosol_names, const bool switch_aerosol_optics);
105+
106+
void create_tilted_columns(const int n_x, const int n_y, const int n_lay_in, const int n_lev_in,
107+
const std::vector<Float>& zh_tilted, const std::vector<ijk>& tilted_path,
108+
std::vector<Float>& var);
109+
110+
void interpolate(const int n_x, const int n_y, const int n_lay_in, const int n_lev_in,
111+
const std::vector<Float>& zh_in, const std::vector<Float>& zf_in,
112+
const std::vector<Float>& play_in, const std::vector<Float>& plev_in,
113+
const Float zp, const ijk offset,
114+
Float* p_out);
115+
116+
void create_tilted_columns_levlay(const int n_x, const int n_y, const int n_lay_in, const int n_lev_in,
117+
const std::vector<Float>& zh_in, const std::vector<Float>& z_in,
118+
const std::vector<Float>& zh_tilted, const std::vector<ijk>& tilted_path,
119+
std::vector<Float>& var_lay, std::vector<Float>& var_lev);
120+
121+
void tica_tilt(const Float sza, const Float azi,
122+
const int n_col_x, const int n_col_y, const int n_col,
123+
const int n_lay, const int n_lev, const int n_z_in, const int n_zh_in ,
124+
Array<Float,1> xh, Array<Float,1> yh, Array<Float,1> zh, Array<Float,1> z,
125+
Array<Float,2> p_lay, Array<Float,2> t_lay, Array<Float,2> p_lev, Array<Float,2> t_lev,
126+
Array<Float,2> lwp, Array<Float,2> iwp, Array<Float,2> rel, Array<Float,2> dei, Array<Float,2> rh,
127+
Gas_concs gas_concs, Aerosol_concs aerosol_concs,
128+
Array<Float,2>& p_lay_out, Array<Float,2>& t_lay_out, Array<Float,2>& p_lev_out, Array<Float,2>& t_lev_out,
129+
Array<Float,2>& lwp_out, Array<Float,2>& iwp_out, Array<Float,2>& rel_out, Array<Float,2>& dei_out, Array<Float,2>& rh_out,
130+
Gas_concs& gas_concs_out, Aerosol_concs aerosol_concs_out,
131+
std::vector<std::string> gas_names, std::vector<std::string> aerosol_names,
132+
bool switch_cloud_optics, bool switch_liq_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics);

‎src_test/CMakeLists.txt‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
#
22
# This file is part of a C++ interface to the Radiative Transfer for Energetics (RTE)
33
#
4-
include_directories(${INCLUDE_DIRS} "../include" "../include_test" "../include_rt" "../include_rt_kernels/" "../include_kernels_cuda")
4+
include_directories(${INCLUDE_DIRS} "../include" "../include_test" "../include_rt" "../include_rt_kernels/" "../include_kernels_cuda" "../include_tilt")
55

66
# retrieve the git hash from the current commit
77
find_package(Git)
@@ -28,7 +28,7 @@ if(USECUDA)
2828
add_executable(test_rte_rrtmgp_gpu Radiation_solver.cu test_rte_rrtmgp.cu)
2929
target_link_libraries(test_rte_rrtmgp_gpu rte_rrtmgp rte_rrtmgp_cuda ${LIBS} m)
3030

31-
add_executable(test_rte_rrtmgp_rt_gpu Radiation_solver_rt.cu test_rte_rrtmgp_rt.cu)
31+
add_executable(test_rte_rrtmgp_rt_gpu Radiation_solver_rt.cu test_rte_rrtmgp_rt.cu ../src_tilt/tilt_utils.cpp)
3232
target_link_libraries(test_rte_rrtmgp_rt_gpu rte_rrtmgp rte_rrtmgp_cuda rte_rrtmgp_cuda_rt curand ${LIBS} m)
3333

3434
add_executable(test_rte_rrtmgp_bw_gpu Radiation_solver_bw.cu test_rte_rrtmgp_bw.cu)
@@ -40,3 +40,4 @@ endif()
4040

4141
add_executable(test_rte_rrtmgp Radiation_solver.cpp test_rte_rrtmgp.cpp)
4242
target_link_libraries(test_rte_rrtmgp rte_rrtmgp ${LIBS} m)
43+

‎src_test/Radiation_solver_rt.cu‎

Lines changed: 38 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,29 @@
4242

4343
namespace
4444
{
45+
__global__
46+
void scale_tau_kernel(Float* tau, const int ncol, const int nlay, Float scale_factor) {
47+
const int icol = blockIdx.x*blockDim.x + threadIdx.x;
48+
const int ilay = blockIdx.y*blockDim.y + threadIdx.y;
49+
50+
if ( (icol < ncol) && (ilay < nlay) )
51+
{
52+
const int idx = icol + ilay*ncol;
53+
tau[idx] = tau[idx] * scale_factor;
54+
}
55+
}
56+
57+
void scale_tau(Float* tau, const int ncol, const int nlay, Float scale_factor) {
58+
const int block_col = 64;
59+
const int block_lay = 1;
60+
const int grid_col = ncol/block_col + (ncol%block_col > 0);
61+
const int grid_lay = nlay/block_lay + (nlay%block_lay > 0);
62+
63+
dim3 grid_gpu(grid_col, grid_lay);
64+
dim3 block_gpu(block_col, block_lay);
65+
scale_tau_kernel<<<grid_gpu, block_gpu>>>(tau, ncol, nlay, scale_factor);
66+
}
67+
4568
std::vector<std::string> get_variable_string(
4669
const std::string& var_name,
4770
std::vector<int> i_count,
@@ -576,6 +599,7 @@ void Radiation_solver_shortwave::solve_gpu(
576599
const bool switch_single_gpt,
577600
const bool switch_delta_cloud,
578601
const bool switch_delta_aerosol,
602+
const bool switch_attenuate_tica,
579603
const int single_gpt,
580604
const Int ray_count,
581605
const Vector<int> grid_cells,
@@ -657,6 +681,7 @@ void Radiation_solver_shortwave::solve_gpu(
657681

658682
const Array<int, 2>& band_limits_gpt(this->kdist_gpu->get_band_lims_gpoint());
659683
int previous_band = 0;
684+
Float attenuate_scale_factor = 1/tsi_scaling({1});
660685

661686
for (int igpt=1; igpt<=n_gpt; ++igpt)
662687
{
@@ -712,6 +737,11 @@ void Radiation_solver_shortwave::solve_gpu(
712737

713738
toa_src.fill(toa_src_temp({1}) * tsi_scaling({1}));
714739

740+
if (switch_attenuate_tica)
741+
{
742+
scale_tau(dynamic_cast<Optical_props_2str_rt&>(*optical_props).get_tau().ptr(), n_col, n_lay, attenuate_scale_factor);
743+
}
744+
715745
if (switch_aerosol_optics)
716746
{
717747
if (band > previous_band)
@@ -724,6 +754,10 @@ void Radiation_solver_shortwave::solve_gpu(
724754
*aerosol_optical_props);
725755
if (switch_delta_aerosol)
726756
aerosol_optical_props->delta_scale();
757+
if (switch_attenuate_tica)
758+
{
759+
scale_tau(dynamic_cast<Optical_props_2str_rt&>(*aerosol_optical_props).get_tau().ptr(), n_col, n_lay, attenuate_scale_factor);
760+
}
727761
}
728762

729763
// Add the cloud optical props to the gas optical properties.
@@ -752,6 +786,10 @@ void Radiation_solver_shortwave::solve_gpu(
752786

753787
if (switch_delta_cloud)
754788
cloud_optical_props->delta_scale();
789+
if (switch_attenuate_tica)
790+
{
791+
scale_tau(dynamic_cast<Optical_props_2str_rt&>(*cloud_optical_props).get_tau().ptr(), n_col, n_lay, attenuate_scale_factor);
792+
}
755793
}
756794
// Add the cloud optical props to the gas optical properties.
757795
add_to(
@@ -765,7 +803,6 @@ void Radiation_solver_shortwave::solve_gpu(
765803
Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_g().ptr());
766804
}
767805

768-
769806
// Store the optical properties, if desired
770807
if (switch_single_gpt && igpt == single_gpt)
771808
{

0 commit comments

Comments
 (0)