Skip to content

Commit 99c83d9

Browse files
MennoVeermanMennoVeerman
authored
Merge pull request #56 from magpowell/tica_fixes
initial commit: fix aerosols , remove unused function, etc.
2 parents 103d7a4 + 96621d4 commit 99c83d9

File tree

3 files changed

+37
-177
lines changed

3 files changed

+37
-177
lines changed

‎include_tilt/tilt_utils.h‎

Lines changed: 0 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -45,14 +45,6 @@ void post_process_output(const std::vector<ColumnResult>& col_results,
4545
const bool switch_liq_cloud_optics,
4646
const bool switch_ice_cloud_optics);
4747

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-
5648
void restore_bkg_profile(const int n_x, const int n_y,
5749
const int n_full,
5850
const int n_tilt,

‎src_test/test_rte_rrtmgp_rt.cu‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -469,6 +469,18 @@ void solve_radiation(int argc, char** argv)
469469
"aermr01", "aermr02", "aermr03", "aermr04", "aermr05", "aermr06", "aermr07",
470470
"aermr08", "aermr09", "aermr10","aermr11"
471471
};
472+
473+
for (const auto& aerosol_name : aerosol_names) {
474+
if (!aerosol_concs.exists(aerosol_name)) {
475+
continue;
476+
}
477+
const Array<Float,2>& gas = aerosol_concs.get_vmr(aerosol_name);
478+
if (gas.size() > 1) {
479+
if (gas.get_dims()[0] == 1){
480+
aerosol_concs.set_vmr(aerosol_name, aerosol_concs.get_vmr(aerosol_name).subset({ {{1,n_col}, {1, n_lay}}} ));
481+
}
482+
}
483+
}
472484

473485
Array<Float,1> xh;
474486
Array<Float,1> yh;

‎src_tilt/tilt_utils.cpp‎

Lines changed: 25 additions & 169 deletions
Original file line numberDiff line numberDiff line change
@@ -138,14 +138,18 @@ void tilted_path(std::vector<Float>& xh, std::vector<Float>& yh,
138138
// Record the path segment
139139
dz_tilted[z_idx] += dz0;
140140

141-
// Check z boundary crossing
142-
if (std::abs(l - lz) < epsilon || zp >= zh[k+1]) {
143-
k += 1;
141+
// Record boundary crossing
142+
if ((std::abs(l - lx) < epsilon) || (std::abs(l - ly) < epsilon) || ( (std::abs(l - lz) < epsilon || zp >= zh[k+1]))){
144143
// Create a new path segment after crossing boundary
145144
tilted_path.push_back({i, j, k});
146145
dz_tilted.push_back(0.0);
147146
z_idx += 1;
148147
}
148+
149+
// Check z boundary crossing
150+
if ((std::abs(l - lz) < epsilon || zp >= zh[k+1])) {
151+
k += 1;
152+
}
149153

150154
// Check y boundary crossing
151155
if (std::abs(l - ly) < epsilon) {
@@ -171,149 +175,6 @@ void tilted_path(std::vector<Float>& xh, std::vector<Float>& yh,
171175
}
172176
}
173177

174-
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,
175-
int n_zh, int n_z,
176-
Float sza, std::vector<Float> zh, std::vector<Float> z,
177-
Array<Float,2> p_lay, Array<Float,2> t_lay, Array<Float,2> p_lev, Array<Float,2> t_lev,
178-
Array<Float,2> lwp, Array<Float,2> iwp, Array<Float,2> rel, Array<Float,2> dei,
179-
Gas_concs& gas_concs, std::vector<std::string> gas_names,
180-
bool switch_cloud_optics, bool switch_liq_cloud_optics, bool switch_ice_cloud_optics) {
181-
const int n_col = n_col_x * n_col_y;
182-
183-
// Create the new NetCDF file
184-
Netcdf_file output_nc(file_name, Netcdf_mode::Create);
185-
// Copy dimensions from the input file
186-
const int n_bnd_sw = 14;
187-
const int n_bnd_lw = 16;
188-
output_nc.add_dimension("band_sw", n_bnd_sw);
189-
output_nc.add_dimension("band_lw", n_bnd_lw);
190-
191-
output_nc.add_dimension("lay", n_lay);
192-
output_nc.add_dimension("lev", n_lev);
193-
194-
output_nc.add_dimension("x", n_col_x);
195-
output_nc.add_dimension("y", n_col_y);
196-
output_nc.add_dimension("z", n_z);
197-
198-
output_nc.add_dimension("xh", n_col_x+1);
199-
output_nc.add_dimension("yh", n_col_y+1);
200-
output_nc.add_dimension("zh", n_zh);
201-
202-
Array<Float,1> grid_x(input_nc.get_variable<Float>("x", {n_col_x}), {n_col_x});
203-
Array<Float,1> grid_xh(input_nc.get_variable<Float>("xh", {n_col_x+1}), {n_col_x+1});
204-
Array<Float,1> grid_y(input_nc.get_variable<Float>("y", {n_col_y}), {n_col_y});
205-
Array<Float,1> grid_yh(input_nc.get_variable<Float>("yh", {n_col_y+1}), {n_col_y+1});
206-
207-
// Create and write the grid coordinate variables
208-
auto nc_x = output_nc.add_variable<Float>("x", {"x"});
209-
auto nc_y = output_nc.add_variable<Float>("y", {"y"});
210-
auto nc_xh = output_nc.add_variable<Float>("xh", {"xh"});
211-
auto nc_yh = output_nc.add_variable<Float>("yh", {"yh"});
212-
213-
nc_x.insert(grid_x.v(), {0});
214-
nc_y.insert(grid_y.v(), {0});
215-
nc_xh.insert(grid_xh.v(), {0});
216-
nc_yh.insert(grid_yh.v(), {0});
217-
218-
// Create and write surface properties
219-
Array<Float,2> sfc_alb_dir(input_nc.get_variable<Float>("sfc_alb_dir", {n_col_y, n_col_x, n_bnd_sw}), {n_bnd_sw, n_col});
220-
Array<Float,2> sfc_alb_dif(input_nc.get_variable<Float>("sfc_alb_dif", {n_col_y, n_col_x, n_bnd_sw}), {n_bnd_sw, n_col});
221-
Array<Float,2> emis_sfc(input_nc.get_variable<Float>("emis_sfc", {n_col_y, n_col_x, n_bnd_lw}), {n_bnd_lw, n_col});
222-
Array<Float,1> t_sfc(input_nc.get_variable<Float>("t_sfc", {n_col_y, n_col_x}), {n_col});
223-
224-
auto nc_sfc_alb_dir = output_nc.add_variable<Float>("sfc_alb_dir", {"y", "x", "band_sw"});
225-
auto nc_sfc_alb_dif = output_nc.add_variable<Float>("sfc_alb_dif", {"y", "x", "band_sw"});
226-
auto nc_emis_sfc = output_nc.add_variable<Float>("emis_sfc", {"y", "x", "band_lw"});
227-
auto nc_t_sfc = output_nc.add_variable<Float>("t_sfc", {"y", "x"});
228-
229-
nc_sfc_alb_dir.insert(sfc_alb_dir.v(), {0, 0, 0});
230-
nc_sfc_alb_dif.insert(sfc_alb_dif.v(), {0, 0, 0});
231-
nc_emis_sfc.insert(emis_sfc.v(), {0, 0, 0});
232-
nc_t_sfc.insert(t_sfc.v(), {0, 0});
233-
234-
// Write tsi scaling
235-
auto nc_tsi_scaling = output_nc.add_variable<Float>("tsi_scaling", {});
236-
nc_tsi_scaling.insert(std::cos(sza), {0});
237-
238-
// Write the grid coordinates
239-
auto nc_z = output_nc.add_variable<Float>("z", {"z"});
240-
auto nc_zh = output_nc.add_variable<Float>("zh", {"zh"});
241-
nc_zh.insert(zh, {0});
242-
nc_z.insert(z, {0});
243-
244-
// Write the atmospheric fields
245-
auto nc_play = output_nc.add_variable<Float>("p_lay", {"lay", "y", "x"});
246-
auto nc_plev = output_nc.add_variable<Float>("p_lev", {"lev", "y", "x"});
247-
248-
nc_play.insert(p_lay.v(), {0, 0, 0});
249-
nc_plev.insert(p_lev.v(), {0, 0, 0});
250-
251-
auto nc_tlay = output_nc.add_variable<Float>("t_lay", {"lay", "y", "x"});
252-
auto nc_tlev = output_nc.add_variable<Float>("t_lev", {"lev", "y", "x"});
253-
254-
nc_tlay.insert(t_lay.v(), {0, 0, 0});
255-
nc_tlev.insert(t_lev.v(), {0, 0, 0});
256-
257-
// Write the cloud optical properties if applicable
258-
if (switch_cloud_optics) {
259-
if (switch_liq_cloud_optics)
260-
{
261-
auto nc_lwp = output_nc.add_variable<Float>("lwp", {"lay", "y", "x"});
262-
nc_lwp.insert(lwp.v(), {0, 0, 0});
263-
auto nc_rel = output_nc.add_variable<Float>("rel", {"lay", "y", "x"});
264-
nc_rel.insert(rel.v(), {0, 0, 0});
265-
}
266-
267-
if (switch_ice_cloud_optics)
268-
{
269-
auto nc_iwp = output_nc.add_variable<Float>("iwp", {"lay", "y", "x"});
270-
auto nc_dei = output_nc.add_variable<Float>("dei", {"lay", "y", "x"});
271-
nc_iwp.insert(iwp.v(), {0, 0, 0});
272-
nc_dei.insert(dei.v(), {0, 0, 0});
273-
}
274-
}
275-
276-
// Write the gas concentrations
277-
for (const auto& gas_name : gas_names) {
278-
if (!gas_concs.exists(gas_name)) {
279-
continue;
280-
}
281-
const Array<Float,2>& gas = gas_concs.get_vmr(gas_name);
282-
std::string var_name = "vmr_" + gas_name;
283-
if (gas.size() == 1) {
284-
auto nc_gas = output_nc.add_variable<Float>(var_name);
285-
nc_gas.insert(gas.v(), {});
286-
} else {
287-
auto nc_gas = output_nc.add_variable<Float>(var_name, {"lay", "y", "x"});
288-
const std::vector<Float>& flat_data = gas.v();
289-
nc_gas.insert(flat_data, {0, 0, 0});
290-
}
291-
}
292-
auto nc_mu = output_nc.add_variable<Float>("mu0", {"y", "x"});
293-
auto nc_azi = output_nc.add_variable<Float>("azi", {"y", "x"});
294-
295-
Array<Float, 2> mu_array(std::array<int, 2>{n_col_y, n_col_x});
296-
Array<Float, 2> azi_array(std::array<int, 2>{n_col_y, n_col_x});
297-
298-
mu_array.fill(1.0);
299-
azi_array.fill(0.0);
300-
301-
nc_mu.insert(mu_array.v(), {0, 0});
302-
nc_azi.insert(azi_array.v(), {0, 0});
303-
304-
auto nc_ng_x = output_nc.add_variable<Float>("ngrid_x", {});
305-
auto nc_ng_y = output_nc.add_variable<Float>("ngrid_y", {});
306-
auto nc_ng_z = output_nc.add_variable<Float>("ngrid_z", {});
307-
308-
nc_ng_x.insert(48, {0});
309-
nc_ng_y.insert(48, {0});
310-
nc_ng_z.insert(32, {0});
311-
312-
output_nc.sync();
313-
return true;
314-
}
315-
316-
317178
void restore_bkg_profile(const int n_x, const int n_y,
318179
const int n_full,
319180
const int n_tilt,
@@ -412,7 +273,6 @@ void restore_bkg_profile_bundle(const int n_col_x, const int n_col_y,
412273
}
413274
if (switch_aerosol_optics)
414275
{
415-
// TODO RH
416276
restore_bkg_profile(n_col_x, n_col_y, n_lay, n_z_in, bkg_start_z, rh_copy->v(), rh->v());
417277
rh_copy->expand_dims({n_col, n_lay_tot});
418278

@@ -662,19 +522,17 @@ void tilt_fields(const int n_z_in, const int n_zh_in, const int n_col_x, const i
662522
continue;
663523
}
664524
const Array<Float,2>& gas = gas_concs_copy.get_vmr(gas_name);
665-
666525
if (gas.size() > 1) {
667-
Array<Float,2> gas_tmp(gas);
668-
create_tilted_columns(n_col_x, n_col_y, n_z_in, n_zh_in, zh_tilt.v(), path.v(), gas_tmp.v());
669-
gas_tmp.expand_dims({n_col, n_z_tilt});
670-
gas_concs_copy.set_vmr(gas_name, gas_tmp);
671-
}
672-
else if (gas.size() == 1) {
673-
// Do nothing for single profiles
526+
if (gas.get_dims()[0] > 1) { // checking: do we have 3D field?
527+
Array<Float,2> gas_tmp(gas);
528+
create_tilted_columns(n_col_x, n_col_y, n_z_in, n_zh_in, zh_tilt.v(), path.v(), gas_tmp.v());
529+
gas_tmp.expand_dims({n_col, n_z_tilt});
530+
gas_concs_copy.set_vmr(gas_name, gas_tmp);
531+
}
532+
else {
533+
throw std::runtime_error("No tilted column implementation for single column profiles.");
534+
}
674535
}
675-
else {
676-
throw std::runtime_error("No tilted column implementation for single profiles.");
677-
}
678536
}
679537

680538
if (switch_aerosol_optics)
@@ -689,17 +547,16 @@ void tilt_fields(const int n_z_in, const int n_zh_in, const int n_col_x, const i
689547
const Array<Float,2>& gas = aerosol_concs_copy.get_vmr(aerosol_name);
690548

691549
if (gas.size() > 1) {
692-
Array<Float,2> gas_tmp(gas);
693-
create_tilted_columns(n_col_x, n_col_y, n_z_in, n_zh_in, zh_tilt.v(), path.v(), gas_tmp.v());
694-
gas_tmp.expand_dims({n_col, n_z_tilt});
695-
aerosol_concs_copy.set_vmr(aerosol_name, gas_tmp);
696-
}
697-
else if (gas.size() == 1) {
698-
// Do nothing for single profiles
550+
if (gas.get_dims()[0] > 1) { // checking: do we have 3D field?
551+
Array<Float,2> gas_tmp(gas);
552+
create_tilted_columns(n_col_x, n_col_y, n_z_in, n_zh_in, zh_tilt.v(), path.v(), gas_tmp.v());
553+
gas_tmp.expand_dims({n_col, n_z_tilt});
554+
aerosol_concs_copy.set_vmr(aerosol_name, gas_tmp);
555+
}
556+
else {
557+
throw std::runtime_error("No tilted column implementation for single column profiles.");
558+
}
699559
}
700-
else {
701-
throw std::runtime_error("No tilted column implementation for single profiles.");
702-
}
703560
}
704561
}
705562

@@ -742,7 +599,6 @@ void compress_fields(const int compress_lay_start_idx, const int n_col_x, const
742599

743600
if (switch_aerosol_optics)
744601
{
745-
// TODO RH
746602
std::vector<Float> ones(p_lay_copy->v().size(), 1.0);
747603
compress_columns_weighted_avg(n_col_x, n_col_y,
748604
n_z_in, n_z_tilt,

0 commit comments

Comments
 (0)