From d4a2087336d09f116399eba53e3912cab36e5050 Mon Sep 17 00:00:00 2001 From: pzmij Date: Tue, 15 Feb 2022 11:38:02 +0100 Subject: [PATCH 1/5] befor pull --- drawbicyc/drawbicyc.cpp | 40 +++++++++---------- drawbicyc/include/plot_prof.hpp | 64 +++++++++++++++---------------- drawbicyc/include/plot_series.hpp | 64 ++++++++++++++++++++++++++++++- drawbicyc/include/plots.hpp | 7 ++++ 4 files changed, 122 insertions(+), 53 deletions(-) diff --git a/drawbicyc/drawbicyc.cpp b/drawbicyc/drawbicyc.cpp index 9fb9313..cada582 100644 --- a/drawbicyc/drawbicyc.cpp +++ b/drawbicyc/drawbicyc.cpp @@ -1,5 +1,5 @@ #include "include/plot_series.hpp" -#include "include/plot_prof.hpp" +///#include "include/plot_prof.hpp" #include "include/plot_fields.hpp" #include "include/plot_qv_qc_2_6_10_min.hpp" #include "include/plot_lgrngn_spec.hpp" @@ -19,7 +19,7 @@ int main(int argc, char** argv) ("qv_qc_2_6_10_min", po::value()->default_value(false) , "plot comparison of qv and qc fields at 2, 6 and 10 min?") ("dir", po::value()->required() , "directory containing out_lgrngn") ("micro", po::value()->required(), "one of: blk_1m, blk_2m, lgrngn") - ("type", po::value()->required(), "one of: dycoms, moist_thermal, rico, Lasher_Trapp, gccn_ccn_conc")//, base_prflux_vs_clhght") + ("type", po::value()->required(), "one of: dycoms, moist_thermal, rico, Lasher_Trapp,ICMW2020_cc, gccn_ccn_conc")//, base_prflux_vs_clhght") ; po::variables_map vm; @@ -33,8 +33,8 @@ int main(int argc, char** argv) // handling the "type" option std::string type = vm["type"].as(); - if(type != "dycoms" && type != "moist_thermal" && type != "rico" && type != "pi_chamber" && type != "Lasher_Trapp" && type != "pi_chamber_icmw" && type != "gccn_ccn_conc") - throw std::runtime_error("Unrecognized 'type' option, only dycoms, rico, moist_thermal, pi_chamber, pi_chamber_icmw, Lasher_Trapp, gccn_ccn_conc available now");//, base_prflux_vs_clhght available now"); + if(type != "dycoms" && type != "moist_thermal" && type != "rico" && type != "pi_chamber" && type != "Lasher_Trapp" && type != "ICMW2020_cc" && type != "pi_chamber_icmw" && type != "gccn_ccn_conc") + throw std::runtime_error("Unrecognized 'type' option, only dycoms, rico, moist_thermal, pi_chamber, pi_chamber_icmw, Lasher_Trapp, ICMW2020_cc, gccn_ccn_conc available now");//, base_prflux_vs_clhght available now"); // should profiles be normalized by inversion height const bool normalize_prof = type == "dycoms"; @@ -56,32 +56,32 @@ int main(int argc, char** argv) H5::DataSet h5d = h5f.openDataSet("G"); H5::DataSpace h5s = h5d.getSpace(); int NDims = h5s.getSimpleExtentNdims(); - + std::cout<<"HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH"<(h5, micro), plots, type); -// if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), plots, type, normalize_prof); -// if(flag_fields) plot_fields(PlotterMicro_t<2>(h5, micro), plots, type); -// if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro)); + if(flag_series) plot_series(PlotterMicro_t<2>(h5, micro), plots, type); + //if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), plots, type, normalize_prof); + //if(flag_fields) plot_fields(PlotterMicro_t<2>(h5, micro), plots, type); + //if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro)); } else if(NDims == 3) { if(flag_series) plot_series(PlotterMicro_t<3>(h5, micro), plots, type); - if(flag_profiles) plot_profiles(PlotterMicro_t<3>(h5, micro), plots, type, normalize_prof); + //if(flag_profiles) plot_profiles(PlotterMicro_t<3>(h5, micro), plots, type, normalize_prof); // if(flag_fields) plot_fields(PlotterMicro_t<3>(h5, micro), plots, type); // if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro)); // if(flag_lgrngn_spec) { diff --git a/drawbicyc/include/plot_prof.hpp b/drawbicyc/include/plot_prof.hpp index 2e7d086..a8a03ea 100644 --- a/drawbicyc/include/plot_prof.hpp +++ b/drawbicyc/include/plot_prof.hpp @@ -890,13 +890,13 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res = res_tmp; res_prof_hlpr = plotter.horizontal_mean(res); // average in x } - else if (plt == "S_drop") // supersat. weighted by nc - { - typename Plotter_t::arr_t nc(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - nc *= rhod * plotter.dv; - typename Plotter_t::arr_t S(plotter.h5load_timestep("RH", at * n["outfreq"])-1); - res_prof_hlpr = plotter.horizontal_weighted_mean(S, nc); - } + //else if (plt == "S_drop") // supersat. weighted by nc + //{ + //typename Plotter_t::arr_t nc(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); + //nc *= rhod * plotter.dv; + //typename Plotter_t::arr_t S(plotter.h5load_timestep("RH", at * n["outfreq"])-1); + //res_prof_hlpr = plotter.horizontal_weighted_mean(S, nc); + //} else if (plt == "Sigma2_S") // variance of supersaturation { typename Plotter_t::arr_t S(plotter.h5load_timestep("RH", at * n["outfreq"])-1); @@ -904,18 +904,18 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res = S * S; res_prof_hlpr = plotter.horizontal_mean(res); // average in x } - else if (plt == "Sigma2_S_drop") // variance of supersaturation weighted by nc - { - typename Plotter_t::arr_t nc(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - nc *= rhod * plotter.dv; - typename Plotter_t::arr_t S(plotter.h5load_timestep("RH", at * n["outfreq"])-1); - - if(blitz::sum(nc) > 0) - { - plotter.subtract_horizontal_weighted_mean(S, nc); - res = S * S; - res_prof_hlpr = plotter.horizontal_weighted_mean(res, nc); - } + //else if (plt == "Sigma2_S_drop") // variance of supersaturation weighted by nc + // { +// typename Plotter_t::arr_t nc(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); +// nc *= rhod * plotter.dv; +// typename Plotter_t::arr_t S(plotter.h5load_timestep("RH", at * n["outfreq"])-1); +// + // if(blitz::sum(nc) > 0) + // { + // plotter.subtract_horizontal_weighted_mean(S, nc); + // res = S * S; + // res_prof_hlpr = plotter.horizontal_weighted_mean(res, nc); + // } else { res_prof_hlpr = 0; @@ -976,23 +976,23 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res_pos_out_done = true; } - if (plt != "base_prflux_vs_clhght") - res_prof_sum /= last_timestep - first_timestep + 1; - else - res_prof_sum = where(occur_no > 0, res_prof_sum/occur_no, 0); + //if (plt != "base_prflux_vs_clhght") + // res_prof_sum /= last_timestep - first_timestep + 1; + //else + // res_prof_sum = where(occur_no > 0, res_prof_sum/occur_no, 0); // set labels for the gnuplot plot - gnuplot_profs_set_labels(gp, plt, normalize); + // gnuplot_profs_set_labels(gp, plt, normalize); // do the plotting - gp << "plot '-' with line\n"; - gp.send1d(boost::make_tuple(res_prof_sum, res_pos)); - - if (plt == "base_prflux_vs_clhght") - { - oprof_file << plt << " number of occurances" << endl; - oprof_file << occur_no; - } + //gp << "plot '-' with line\n"; + //gp.send1d(boost::make_tuple(res_prof_sum, res_pos)); + +// if (plt == "base_prflux_vs_clhght") +// { +// oprof_file << plt << " number of occurances" << endl; +// oprof_file << occur_no; + // } oprof_file << plt << endl; oprof_file << res_prof_sum ; diff --git a/drawbicyc/include/plot_series.hpp b/drawbicyc/include/plot_series.hpp index d2ac874..2d82299 100644 --- a/drawbicyc/include/plot_series.hpp +++ b/drawbicyc/include/plot_series.hpp @@ -794,7 +794,69 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) } } catch(...) {if(at==first_timestep) data_found[plt]=0;} - } + } + // else if (plt == "cwp") + // { + // // cloud water path + // try + // { +// { +// typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); +// snap *= rhod *1e3; // water per cubic metre (should be wet density ...) & g/kg +// res_prof(at) = blitz::mean(snap); +// } +// } +// catch(...) {;} +// } +// else if (plt == "lwm") +// { +// //liquid water mass +// try +// { +// { + // auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]) * rhod; +// typename Plotter_t::arr_t snap(tmp); +// snap += plotter.h5load_rr_timestep(at * n["outfreq"]) * rhod; +// snap *= plotter.CellVol; +// snap(plotter.hrzntl_slice(0)) = snap(plotter.hrzntl_slice(0))/2; +// snap(plotter.hrzntl_slice(-1)) = snap(plotter.hrzntl_slice(-1))/2; +// res_prof(at) = blitz::sum( snap); + // } + // } + // catch(...) {;} + // } + //else if (plt == "cwm") + //{ + //cloud water mass +// try +// { +// { + // auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]); +// typename Plotter_t::arr_t snap(tmp); + // snap *= rhod * plotter.CellVol; + // snap(plotter.hrzntl_slice(0)) = snap(plotter.hrzntl_slice(0))/2; + // snap(plotter.hrzntl_slice(-1)) = snap(plotter.hrzntl_slice(-1))/2; + // res_prof(at) = blitz::sum(snap); + // } +// } +// catch(...) {;} + // } + // else if (plt == "rwm") + // { +// //liquid water mass +// try +// { +// { +// auto tmp = plotter.h5load_rr_timestep(at * n["outfreq"]); + // typename Plotter_t::arr_t snap(tmp); + // snap *= rhod * plotter.CellVol; + // snap(plotter.hrzntl_slice(0)) = snap(plotter.hrzntl_slice(0))/2; + // snap(plotter.hrzntl_slice(-1)) = snap(plotter.hrzntl_slice(-1))/2; + // res_prof(at) = blitz::sum(snap); +// } +// } +// catch(...) {;} + // } else if (plt == "surf_flux_latent") { try diff --git a/drawbicyc/include/plots.hpp b/drawbicyc/include/plots.hpp index 619b592..f0baabd 100644 --- a/drawbicyc/include/plots.hpp +++ b/drawbicyc/include/plots.hpp @@ -7,6 +7,8 @@ #include "cases/PiChamber/plots_icmw.hpp" #include "cases/Lasher_Trapp/plots.hpp" #include "cases/common_plots/GCCN_CCN_conc/plots.hpp" +#include "cases/ICMW2020_cumulus_congestus/plots.hpp" + const std::vector series_sgs({ // "tot_tke" @@ -58,6 +60,11 @@ class Plots series.insert(series.end(), series_Lasher_Trapp.begin(), series_Lasher_Trapp.end()); fields.insert(fields.end(), fields_Lasher_Trapp.begin(), fields_Lasher_Trapp.end()); } + else if(type == "ICMW2020_cc") { + profs.insert(profs.end(), profs_ICMW2020.begin(), profs_ICMW2020.end()); + series.insert(series.end(), series_ICMW2020.begin(), series_ICMW2020.end()); + fields.insert(fields.end(), fields_ICMW2020.begin(), fields_ICMW2020.end()); + } else if(type == "pi_chamber_icmw") { profs.insert(profs.end(), profs_PiChamberICMW.begin(), profs_PiChamberICMW.end()); series.insert(series.end(), series_PiChamberICMW.begin(), series_PiChamberICMW.end()); From e7900f4847a57ae2fee8935e2de7a43b22bb6648 Mon Sep 17 00:00:00 2001 From: pzmij Date: Wed, 16 Feb 2022 01:23:30 +0100 Subject: [PATCH 2/5] Cumu_Cong update --- .../Cumulus_Congestus_series_comparison.py | 74 +++++++++++ cases/Cumulus_Congestus/plot_ranges.py | 124 ++++++++++++++++++ drawbicyc/average.cpp | 2 +- drawbicyc/drawbicyc.cpp | 11 +- .../include/cases/Cumulus_Congestus/plots.hpp | 93 +++++++++++++ drawbicyc/include/plot_prof.hpp | 62 ++++----- 6 files changed, 328 insertions(+), 38 deletions(-) create mode 100644 cases/Cumulus_Congestus/Cumulus_Congestus_series_comparison.py create mode 100644 cases/Cumulus_Congestus/plot_ranges.py create mode 100644 drawbicyc/include/cases/Cumulus_Congestus/plots.hpp diff --git a/cases/Cumulus_Congestus/Cumulus_Congestus_series_comparison.py b/cases/Cumulus_Congestus/Cumulus_Congestus_series_comparison.py new file mode 100644 index 0000000..40ac10b --- /dev/null +++ b/cases/Cumulus_Congestus/Cumulus_Congestus_series_comparison.py @@ -0,0 +1,74 @@ +from matplotlib import rc +import matplotlib.pyplot as plt +from matplotlib.ticker import AutoMinorLocator, MultipleLocator + +import os +import sys +sys.path.append(os.path.dirname(os.path.abspath(__file__)) + "/../Matplotlib_common/") + +from plot_ranges import xscaledict, yscaledict, xlimdict_series, ylimdict_series +from plot_series import * + +# activate latex text rendering +rc('text', usetex=True) + +Cumulus_Congestus_vars = ["RH_max", "cloud_top_height", "surf_precip", "acc_precip", "acc_vol_precip"] + +# init the plot +nplotx = 2 +nploty= 3 +fig, axarr = plt.subplots(nplotx,nploty) + +plot_series( Cumulus_Congestus_vars, 0, nplotx, nploty, axarr, xscaledict, yscaledict, xlimdict_series, ylimdict_series, xlabel='Time [h]') + +# legend font size +plt.rcParams.update({'font.size': 8}) + +# hide axes on empty plots +if len( Cumulus_Congestus_vars) % nploty == 0: + nemptyplots = 0 +else: + nemptyplots = nploty - len( Cumulus_Congestus_vars) % nploty + emptyplots = np.arange(nploty - nemptyplots, nploty) +for empty in emptyplots: + axarr[nplotx-1, empty].axis('off') + +# hide hrzntl tic labels +x_empty_label = np.arange(0, nplotx-1) +y_empty_label = np.arange(nploty) +for x in x_empty_label: + for y in y_empty_label: + axarr[x,y].set_xticklabels([]) + +#axes = plt.gca() +#axes.tick_params(direction='in') +x_arr = np.arange(nplotx) +y_arr = np.arange(nploty) +for x in x_arr: + for y in y_arr: + #tics inside + axarr[x,y].tick_params(direction='in', which='both', top=1, right=1) + #minor tics + axarr[x,y].xaxis.set_minor_locator(AutoMinorLocator()) + axarr[x,y].yaxis.set_minor_locator(AutoMinorLocator()) + #labels and tics font size + for item in ([axarr[x,y].xaxis.label, axarr[x,y].yaxis.label] + axarr[x,y].get_xticklabels() + axarr[x,y].get_yticklabels()): + item.set_fontsize(8) + # subplot numbering + if y < nploty - nemptyplots or x < (nplotx - 1): + axarr[x,y].text(0.2, 0.9, labeldict[y + x*nploty], fontsize=8, transform=axarr[x,y].transAxes) + +#single legend for the whole figure +handles, labels = axarr[0,0].get_legend_handles_labels() + +lgd = fig.legend(handles, labels, handlelength=4, loc='lower center', bbox_to_anchor=(0.45,0)) + +#figure size +fig.set_size_inches(7.874, 5. + (len(labels) - 2) * 0.2)# 5.214)#20.75,13.74) +#distances between subplots and from bottom of the plot +fig.subplots_adjust(bottom=0.18 + (len(labels) - 2) * 0.03, hspace=0.1, wspace=0.4) + +#fig.tight_layout(pad=0, w_pad=0, h_pad=0) + +#plt.show() +fig.savefig(argv[len(sys.argv)-1], bbox_inches='tight', dpi=300)#, bbox_extra_artists=(lgd,)) diff --git a/cases/Cumulus_Congestus/plot_ranges.py b/cases/Cumulus_Congestus/plot_ranges.py new file mode 100644 index 0000000..f4e2d51 --- /dev/null +++ b/cases/Cumulus_Congestus/plot_ranges.py @@ -0,0 +1,124 @@ +xscaledict = { + "thl" : "linear", + "00rtot" : "linear", + "rliq" : "linear", + "prflux" : "linear", + "cl_nc" : "linear", + "clfrac" : "linear", + "wvar" : "linear", + "w3rd" : "linear", + "sat_RH" : "linear", + "rad_flx" : "linear", + "lwp" : "linear", + "RH_max" : "linear", + "cloud_top_height" : "linear", + "total_droplets_number" : "linear", + "acc_precip_vol" : "linear", + "rwp" : "linear", + "er" : "linear", + "wvarmax" : "linear", + "surf_precip" : "linear", + "acc_precip" : "linear", + "cloud_base" : "linear", + "gccn_rw_cl" : "linear", + "non_gccn_rw_cl" : "linear", + "base_prflux_vs_clhght" : "log", + "cl_gccn_conc" : "linear" + } + +yscaledict = { + "thl" : "linear", + "00rtot" : "linear", + "rliq" : "linear", + "prflux" : "linear", + "cl_nc" : "linear", + "clfrac" : "linear", + "wvar" : "linear", + "w3rd" : "linear", + "sat_RH" : "linear", + "rad_flx" : "linear", + "lwp" : "linear", + "RH_max" : "linear", + "cloud_top_height" : "linear", + "total_droplets_number" : "linear", + "acc_precip_vol" : "linear", + "rwp" : "linear", + "er" : "linear", + "wvarmax" : "linear", + "surf_precip" : "linear", + "acc_precip" : "linear", + "cloud_base" : "linear", + "gccn_rw_cl" : "linear", + "non_gccn_rw_cl" : "linear", + "base_prflux_vs_clhght" : "linear", + "cl_gccn_conc" : "log" + } + +xlimdict_profs = { + "thl" : None, + "00rtot" : None, + "rliq" : None, + "prflux" : (0,20), + "cl_nc" : (0,90), + "clfrac" : None, + "wvar" : None, + "w3rd" : None, + "sat_RH" : None, + #"RH_max" : None, + #"cloud_top_height" : None, + "rad_flx" : None, + "gccn_rw_cl" : (0,90), + "non_gccn_rw_cl" : (0,12), + "base_prflux_vs_clhght" : (1,10000) + } +ylimdict_profs = { + "thl" : (0,3000), + "00rtot" : (0,3000), + "rliq" : (0,3000), + "prflux" : (0,3000), + "cl_nc" : (0,3000), + "clfrac" : (0,3000), + "wvar" : (0,3000), + "w3rd" : (0,3000), + "sat_RH" : (0,3000), + "RH_max" : (0,3000), + "cloud_top_height" : (0,3000), + "rad_flx" : (0,3000), + "gccn_rw_cl" : (0,3000), + "non_gccn_rw_cl" : (0,3000), + "base_prflux_vs_clhght" : (0,2500) + } + +xlimdict_series = { + "clfrac" : None, + "cl_nc" : None, + "lwp" : None, + "RH_max" : None, + "cloud_top_height" : None, + "total_droplets_number" : None, + "rwp" : None, + "er" : None, + "wvarmax" : None, + "surf_precip" : None, + "acc_precip" : None, + "acc_precip_vol" : None, + "cl_gccn_conc" : None, + "cloud_base" : None + } + +ylimdict_series = { + "clfrac" : None, + "cl_nc" : None, + "lwp" : None, + "RH_max" : None, + "cloud_top_height" : None, + "total_droplets_number" : None, + "acc_precip_vol" : None, + "rwp" : None, + "er" : None, + "wvarmax" : None, + "surf_precip" : None, + "acc_precip" : None,#(0,0.07), + "cl_gccn_conc" : (1e-6, 1), + "cloud_base" : None +} diff --git a/drawbicyc/average.cpp b/drawbicyc/average.cpp index 4fccc4c..5dca5f2 100644 --- a/drawbicyc/average.cpp +++ b/drawbicyc/average.cpp @@ -193,7 +193,7 @@ int main(int argc, char* argv[]) const string profiles_suffix = string("_profiles_")+argv[2]+string("_")+argv[3]; string plot_type; // determine type of plots based on the name of the first file - const string types[] = {"rico", "dycoms", "moist_thermal", "Cumulus_Congestus"}; + const string types[] = {"rico", "dycoms", "moist_thermal", "cumulus_congestus"}; for(auto type : types) { if(hasEnding(string(argv[4]), type)) diff --git a/drawbicyc/drawbicyc.cpp b/drawbicyc/drawbicyc.cpp index 6954c1a..c87ee7d 100644 --- a/drawbicyc/drawbicyc.cpp +++ b/drawbicyc/drawbicyc.cpp @@ -56,8 +56,7 @@ int main(int argc, char** argv) H5::DataSet h5d = h5f.openDataSet("G"); H5::DataSpace h5s = h5d.getSpace(); int NDims = h5s.getSimpleExtentNdims(); - std::cout<<"HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH"<(h5, micro), plots, type); //if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), plots, type, normalize_prof); - if(flag_fields) plot_fields(PlotterMicro_t<2>(h5, micro), plots, type); + //if(flag_fields) plot_fields(PlotterMicro_t<2>(h5, micro), plots, type); //if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro)); } else if(NDims == 3) { if(flag_series) plot_series(PlotterMicro_t<3>(h5, micro), plots, type); //if(flag_profiles) plot_profiles(PlotterMicro_t<3>(h5, micro), plots, type, normalize_prof); -// if(flag_fields) plot_fields(PlotterMicro_t<3>(h5, micro), plots, type); -// if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro)); + //if(flag_fields) plot_fields(PlotterMicro_t<3>(h5, micro), plots, type); + //if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro)); // if(flag_lgrngn_spec) { // plot_lgrngn_spec_positions(PlotterMicro_t<3>(h5, "lgrngn")); // plot_lgrngn_spec(PlotterMicro_t<3>(h5, "lgrngn")); diff --git a/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp b/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp new file mode 100644 index 0000000..bff92bd --- /dev/null +++ b/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp @@ -0,0 +1,93 @@ +#pragma once + + +const std::vector series_Cumulus_Congestus({ + //"clfrac", + "lwp", + "rwp", + "cwp", + "lwm", + "cwm", + "rwm", + //"surf_precip", + "acc_precip", + "acc_vol_precip", + //"cl_nc", + //"cloud_base", + //"RH_max", + "cloud_top_height", + "total_droplets_number", + //"surf_flux_latent", + //"surf_flux_sensible", + //"sd_conc_avg", + //"mass_dry", + //"cl_gccn_conc", + //"gccn_conc", + //"cl_non_gccn_conc", + //"non_gccn_conc", + //"cl_gccn_to_non_gccn_conc_ratio" + //, "cl_gccn_meanr" + //,"cl_avg_cloud_rad" + // "sd_conc_std_dev", + // "tot_water" +}); + +std::vector profs_Cumulus_Congestus({ + "00rtot", + "rliq", + //"thl", + //"wvar", + //"prflux", + //"clfrac", + //"sd_conc", + "cl_nc", + //"cl_nc_up", + //"w", + //"u", + //"v", + //"base_prflux_vs_clhght", + //"non_gccn_rw_cl", + //"gccn_rw_cl," + //, "N_c", + "actrw_reff_cl", + "ratio_mean_volume_r_to_eff_r_cubed", + "cloud_water_std", + "rain_water_std", + "cloud_std_dev" + //,"vel_div" + //, "nc_up" + //,"sat_RH_up" + //, "act_conc_up" + //, "nc_down" +}); + +std::vector fields_Cumulus_Congestus({ + //"rl", + "mass_conc_rain_drops", + "mass_conc_cloud_droplets", + "num_conc_droplets", + "LWC", + //"mass_conc_cloud", + "ratio_mean_volume_r_to_eff_r_cubed", + "std_dev_droplet_size_dist", + //"rr", + "num_conc_rain_drops", + "effective_radius", + //"na", + //"th", + //"rv" + //"u", + //"w" + //"sd_conc", + //"r_dry", + //"RH", + //"supersat", + //"nc", + //"rr", + //"nr", + //"ef", "na", + //"th", + //"rv", + //"u", + //"w" +}); diff --git a/drawbicyc/include/plot_prof.hpp b/drawbicyc/include/plot_prof.hpp index 3fae9f1..deeed69 100644 --- a/drawbicyc/include/plot_prof.hpp +++ b/drawbicyc/include/plot_prof.hpp @@ -978,13 +978,13 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res = res_tmp; res_prof_hlpr = plotter.horizontal_mean(res); // average in x } - //else if (plt == "S_drop") // supersat. weighted by nc - //{ - //typename Plotter_t::arr_t nc(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - //nc *= rhod * plotter.dv; - //typename Plotter_t::arr_t S(plotter.h5load_timestep("RH", at * n["outfreq"])-1); - //res_prof_hlpr = plotter.horizontal_weighted_mean(S, nc); - //} + else if (plt == "S_drop") // supersat. weighted by nc + { + typename Plotter_t::arr_t nc(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); + nc *= rhod * plotter.dv; + typename Plotter_t::arr_t S(plotter.h5load_timestep("RH", at * n["outfreq"])-1); + res_prof_hlpr = plotter.horizontal_weighted_mean(S, nc); + } else if (plt == "Sigma2_S") // variance of supersaturation { typename Plotter_t::arr_t S(plotter.h5load_timestep("RH", at * n["outfreq"])-1); @@ -992,18 +992,18 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res = S * S; res_prof_hlpr = plotter.horizontal_mean(res); // average in x } - //else if (plt == "Sigma2_S_drop") // variance of supersaturation weighted by nc - // { -// typename Plotter_t::arr_t nc(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); -// nc *= rhod * plotter.dv; -// typename Plotter_t::arr_t S(plotter.h5load_timestep("RH", at * n["outfreq"])-1); -// - // if(blitz::sum(nc) > 0) - // { - // plotter.subtract_horizontal_weighted_mean(S, nc); - // res = S * S; - // res_prof_hlpr = plotter.horizontal_weighted_mean(res, nc); - // } + else if (plt == "Sigma2_S_drop") // variance of supersaturation weighted by nc + { + typename Plotter_t::arr_t nc(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); + nc *= rhod * plotter.dv; + typename Plotter_t::arr_t S(plotter.h5load_timestep("RH", at * n["outfreq"])-1); + + if(blitz::sum(nc) > 0) + { + plotter.subtract_horizontal_weighted_mean(S, nc); + res = S * S; + res_prof_hlpr = plotter.horizontal_weighted_mean(res, nc); + } else { res_prof_hlpr = 0; @@ -1064,23 +1064,23 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res_pos_out_done = true; } - //if (plt != "base_prflux_vs_clhght") - // res_prof_sum /= last_timestep - first_timestep + 1; - //else - // res_prof_sum = where(occur_no > 0, res_prof_sum/occur_no, 0); + if (plt != "base_prflux_vs_clhght") + res_prof_sum /= last_timestep - first_timestep + 1; + else + res_prof_sum = where(occur_no > 0, res_prof_sum/occur_no, 0); // set labels for the gnuplot plot - // gnuplot_profs_set_labels(gp, plt, normalize); + gnuplot_profs_set_labels(gp, plt, normalize); // do the plotting - //gp << "plot '-' with line\n"; - //gp.send1d(boost::make_tuple(res_prof_sum, res_pos)); + gp << "plot '-' with line\n"; + gp.send1d(boost::make_tuple(res_prof_sum, res_pos)); -// if (plt == "base_prflux_vs_clhght") -// { -// oprof_file << plt << " number of occurances" << endl; -// oprof_file << occur_no; - // } + if (plt == "base_prflux_vs_clhght") + { + oprof_file << plt << " number of occurances" << endl; + oprof_file << occur_no; + } oprof_file << plt << endl; oprof_file << res_prof_sum ; From 5d8c0cc65dfc52fe3b12aa5f9a98c4a2bf5851ab Mon Sep 17 00:00:00 2001 From: ryba183 Date: Sun, 20 Mar 2022 20:50:30 +0100 Subject: [PATCH 3/5] update CC --- drawbicyc/drawbicyc.cpp | 4 +- .../include/cases/Cumulus_Congestus/plots.hpp | 15 +- .../include/gnuplot_profs_set_labels.hpp | 4 + drawbicyc/include/plot_prof.hpp | 138 +++++++++++------- 4 files changed, 99 insertions(+), 62 deletions(-) diff --git a/drawbicyc/drawbicyc.cpp b/drawbicyc/drawbicyc.cpp index c87ee7d..f1d9bfb 100644 --- a/drawbicyc/drawbicyc.cpp +++ b/drawbicyc/drawbicyc.cpp @@ -1,5 +1,5 @@ #include "include/plot_series.hpp" -//#include "include/plot_prof.hpp" +#include "include/plot_prof.hpp" #include "include/plot_fields.hpp" #include "include/plot_qv_qc_2_6_10_min.hpp" #include "include/plot_lgrngn_spec.hpp" @@ -73,7 +73,7 @@ int main(int argc, char** argv) { if(flag_series) plot_series(PlotterMicro_t<2>(h5, micro), plots, type); - //if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), plots, type, normalize_prof); + if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), plots, type, normalize_prof); //if(flag_fields) plot_fields(PlotterMicro_t<2>(h5, micro), plots, type); //if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro)); } diff --git a/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp b/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp index bff92bd..0431945 100644 --- a/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp +++ b/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp @@ -37,7 +37,7 @@ std::vector profs_Cumulus_Congestus({ "rliq", //"thl", //"wvar", - //"prflux", + "prflux", //"clfrac", //"sd_conc", "cl_nc", @@ -45,15 +45,16 @@ std::vector profs_Cumulus_Congestus({ //"w", //"u", //"v", - //"base_prflux_vs_clhght", + "coal_tele", + "base_prflux_vs_clhght", //"non_gccn_rw_cl", //"gccn_rw_cl," //, "N_c", - "actrw_reff_cl", - "ratio_mean_volume_r_to_eff_r_cubed", - "cloud_water_std", - "rain_water_std", - "cloud_std_dev" + // "actrw_reff_cl", + // "ratio_mean_volume_r_to_eff_r_cubed", + // "cloud_water_std", + //"rain_water_std", + // "cloud_std_dev" //,"vel_div" //, "nc_up" //,"sat_RH_up" diff --git a/drawbicyc/include/gnuplot_profs_set_labels.hpp b/drawbicyc/include/gnuplot_profs_set_labels.hpp index f2cfd8e..b6448b4 100644 --- a/drawbicyc/include/gnuplot_profs_set_labels.hpp +++ b/drawbicyc/include/gnuplot_profs_set_labels.hpp @@ -162,6 +162,10 @@ void gnuplot_profs_set_labels(Gnuplot &gp, std::string plt, const bool normalize { gp << "set title 'precipitation flux [W/m^2]'\n"; } + else if (plt == "coal_tele") + { + gp << "set title 'coal tele mass flux [W/m^2]'\n"; + } else if (plt == "rad_flx") { gp << "set title 'radiative flux [W/m2]'\n"; diff --git a/drawbicyc/include/plot_prof.hpp b/drawbicyc/include/plot_prof.hpp index deeed69..14cce55 100644 --- a/drawbicyc/include/plot_prof.hpp +++ b/drawbicyc/include/plot_prof.hpp @@ -42,7 +42,8 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool std::cerr << int(n["dt"] * n["outfreq"]+0.5) << std::endl; int first_timestep = vm["prof_start"].as() / int(n["dt"] * n["outfreq"]+0.5); int last_timestep = vm["prof_end"].as() / int(n["dt"] * n["outfreq"]+0.5); - + //int start_time = vm["prof_start"].as(); + //int end_time = vm["prof_end"].as(); // some ugly constants const double p_1000 = 100000.; const double L = 2.5e6; @@ -68,17 +69,38 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool blitz::secondIndex j; typename Plotter_t::arr_t res(rhod.shape()); typename Plotter_t::arr_t res_tmp(rhod.shape()); + //typename Plotter_t::arr_t res_tmp1(rhod.shape()); typename Plotter_t::arr_t res_tmp2(rhod.shape()); + //typename Plotter_t::arr_t res_tmp3(rhod.shape()); + //typename Plotter_t::arr_t res_tmp4(rhod.shape()); + //typename Plotter_t::arr_t res_sum(rhod.shape()); + //typename Plotter_t::arr_t res_mean(rhod.shape()); + //typename Plotter_t::arr_t rc_mask(rhod.shape()); blitz::Array res_prof_sum(n["z"]); // profile interpolate to the uniform grid summed over timesteps blitz::Array res_prof(n["z"]); // profile interpolate to the uniform grid blitz::Array res_prof_hlpr(n["z"]); // actual profile blitz::Array prof_tmp(n["z"]); + //blitz::Array res_prof_num(n["z"]); blitz::Array occur_no(n["z"]); // number of occurances - for unusual profiles like base_prflux_vs_clhght blitz::Range all = blitz::Range::all(); res_prof_sum = 0; occur_no = 0; - + + if (plt == "coal_tele") + { + //coal tele mass flux in W/m^2 + //res = plotter.h5load_timestep("coal_tele_mass_flux",first_timestep * n["outfreq"]) * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * first_timestep * n["outfreq"] * n["dt"]); + res = plotter.h5load_timestep("coal_tele_mass_flux",last_timestep * n["outfreq"]) * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * last_timestep * n["outfreq"] * n["dt"]); + res -= plotter.h5load_timestep("coal_tele_mass_flux",first_timestep * n["outfreq"]) * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * first_timestep * n["outfreq"] * n["dt"]); + // - plotter.h5load_timestep("coal_tele_mass_flux",first_timestep * n["outfreq"]) * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * first_timestep * n["outfreq"] * n["dt"]); + res_prof_hlpr = plotter.horizontal_mean(res); + //auto tmp2 = plotter.h5load_timestep("coal_tele_mass_flux", last_timestep * n["outfreq"]); + //typename Plotter_t::arr_t snap(tmp); + //typename Plotter_t::arr_t snap2(tmp2); + //res_prof_sum = (res_prof * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * last_timestep * n["outfreq"] * n["dt"])) - (res_prof_hlpr * 4./3 * 3.14 * 1e3 * 2264.76e3 / (plotter.CellVol * first_timestep * n["outfreq"] * n["dt"])); + //res = tmp; + } for (int at = first_timestep; at <= last_timestep; ++at) // TODO: mark what time does it actually mean! { std::cout << at * n["outfreq"] << std::endl; @@ -97,33 +119,44 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool // liquid water content res += plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud res += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; // rain - typename Plotter_t::arr_t rc_mask(plotter.h5load_rc_timestep(at * n["outfreq"])); - rc_mask = iscloudy_rc_rico(rc_mask); - prof_tmp = plotter.horizontal_sum(rc_mask); - res_prof_hlpr = where(prof_tmp > 0, plotter.horizontal_sum(res) / prof_tmp, 0); + typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); + res_tmp = iscloudy_rc_rico(snap); + res_prof_sum = plotter.horizontal_sum(res_tmp); + res_prof_hlpr = where(prof_tmp > 0, plotter.horizontal_sum(res) / res_prof_sum, 0); } - if (plt == "cloud_water_std") +/* if (plt == "cloud_water_std") { // cloud_water_std - res = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud - rc_mask = iscloudy_rc_rico(res); - res_prof = res * rc_mask; - res_prof_num = plotter.horizontal_sum(rc_mask); - res_mean = plotter.horizontal_sum(res_prof) / res_prof_num; - res_sum = plotter.horizontal_sum((res_prof - res_mean) * (res_prof - res_mean)) / (res_prof_num -1); - res_prof_hlpr = where(rc_mask >0 , sqrt(res_sum) , 0); + //auto res = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud + //typename Plotter_t::arr_t snap(res); + res = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud + res_prof_sum = plotter.horizontal_sum(res); + res *= iscloudy_rc_rico(res); + res_tmp = plotter.horizontal_sum(res_tmp) / res_prof_sum; + res_tmp2 = plotter.horizontal_sum((res - res_tmp) * (res - res_tmp)) / (res_prof_sum -1); + res_prof_hlpr = where(res >0, sqrt(res_tmp2), 0); + + //rc_mask = iscloudy_rc_rico(res); + //res_tmp = res * rc_mask; + //res_prof_num = plotter.horizontal_sum(rc_mask); + //res_mean = plotter.horizontal_sum(res_tmp) / res_prof_num; + //res_sum = plotter.horizontal_sum((res_tmp - res_mean) * (res_tmp - res_mean)) / (res_prof_num -1); + //res_prof_hlpr = where(rc_mask >0 , sqrt(res_sum) , 0); } - if (plt == "rain_water_std") - { +*/ + //if (plt == "rain_water_std") + //{ // rain_water_std - res = plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; // cloud - rc_mask = iscloudy_rc_rico(res); - res_prof = res * rc_mask; - res_prof_num = plotter.horizontal_sum(rc_mask); - res_mean = plotter.horizontal_sum(res_prof) / res_prof_num; - res_sum = plotter.horizontal_sum((res_prof - res_mean) * (res_prof - res_mean)) / (res_prof_num -1); - res_prof_hlpr = where(rc_mask >0 , sqrt(res_sum) , 0); - } + // auto tmp = plotter.h5load_rr_timestep(at * n["outfreq"]) + + // res = plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; // cloud + // rc_mask = iscloudy_rc_rico(res); + // res_tmp = res * rc_mask; + // res_prof_num = plotter.horizontal_sum(rc_mask); + // res_mean = plotter.horizontal_sum(res_tmp) / res_prof_num; + // res_sum = plotter.horizontal_sum((res_prof - res_mean) * (res_prof - res_mean)) / (res_prof_num -1); + // res_prof_hlpr = where(rc_mask >0 , sqrt(res_sum) , 0); + //} if (plt == "gccn_conc") { res = plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"]) * rhod / 1e6; @@ -463,35 +496,35 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool prof_tmp = plotter.horizontal_sum(res_tmp2); // number of downdraft cells on a given level res_prof_hlpr = where(prof_tmp > 0 , plotter.horizontal_sum(res_tmp) / prof_tmp, 0); } - if (plt == "ratio_mean_volume_r_to_eff_r_cubed") - { + //if (plt == "ratio_mean_volume_r_to_eff_r_cubed") + //{ // ratio mean (r volume / r effective) ^3 - { - typename Plotter_t::arr_t snap(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"])); - res_tmp = snap; - } - { - auto tmp2 = plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp2); - res_tmp2 = snap; - } - { - auto tmp3 = plotter.h5load_timestep("actrw_rw_mom3", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp3); - res_tmp3 = snap; - } - { - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - res_tmp4 = iscloudy_rc_rico(snap); - } + //{ + // typename Plotter_t::arr_t snap(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"])); + // res_tmp = snap; + // } + // { + // auto tmp2 = plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]); + // typename Plotter_t::arr_t snap(tmp2); + // res_tmp2 = snap; + // } + // { + // auto tmp3 = plotter.h5load_timestep("actrw_rw_mom3", at * n["outfreq"]); + // typename Plotter_t::arr_t snap(tmp3); + // res_tmp3 = snap; + // } + // { + // typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); + // res_tmp4 = iscloudy_rc_rico(snap); + // } - res_tmp = where(res_tmp > 0, res_tmp2 * res_tmp2 * res_tmp2 / res_tmp, 0. ); - res_tmp = where(res_tmp3 > 0, res_tmp / res_tmp3 / res_tmp3, 0.); - res_tmp *= res_tmp4; - prof_tmp = plotter.horizontal_sum(res_tmp4); // number of downdraft cells on a given level - res_prof_hlpr = where(prof_tmp > 0 , plotter.horizontal_sum(res_tmp) / prof_tmp, 0); - } - if (plt == "cloud_std_dev") + // res_tmp = where(res_tmp > 0, res_tmp2 * res_tmp2 * res_tmp2 / res_tmp, 0. ); + // res_tmp = where(res_tmp3 > 0, res_tmp / res_tmp3 / res_tmp3, 0.); + // res_tmp *= res_tmp4; + // prof_tmp = plotter.horizontal_sum(res_tmp4); // number of downdraft cells on a given level + // res_prof_hlpr = where(prof_tmp > 0 , plotter.horizontal_sum(res_tmp) / prof_tmp, 0); + // } + /* if (plt == "cloud_std_dev") { { auto tmp1 = plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]) *1e6; @@ -519,6 +552,7 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool prof_tmp = plotter.horizontal_sum(res_tmp3); // number of downdraft cells on a given level res_prof_hlpr = where(prof_tmp > 0 , plotter.horizontal_sum(res_tmp) / prof_tmp, 0); } + */ if (plt == "nc_up") { // updraft only @@ -1063,7 +1097,6 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool oprof_file << res_pos; res_pos_out_done = true; } - if (plt != "base_prflux_vs_clhght") res_prof_sum /= last_timestep - first_timestep + 1; else @@ -1075,7 +1108,6 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool // do the plotting gp << "plot '-' with line\n"; gp.send1d(boost::make_tuple(res_prof_sum, res_pos)); - if (plt == "base_prflux_vs_clhght") { oprof_file << plt << " number of occurances" << endl; From 567c0f41c4b73d7745d4df7cf9f66793498988b4 Mon Sep 17 00:00:00 2001 From: ryba183 Date: Thu, 1 Sep 2022 14:42:58 +0200 Subject: [PATCH 4/5] 01.09.2022 working --- UWLCM_plotters/include/PlotterCommon.hpp | 10 ++--- drawbicyc/drawbicyc.cpp | 17 ++++--- .../include/cases/Cumulus_Congestus/plots.hpp | 3 +- .../include/gnuplot_series_set_labels.hpp | 44 ++++++++++++++++++- drawbicyc/include/plot_series.hpp | 2 +- 5 files changed, 60 insertions(+), 16 deletions(-) diff --git a/UWLCM_plotters/include/PlotterCommon.hpp b/UWLCM_plotters/include/PlotterCommon.hpp index 1c75c83..cc34d00 100644 --- a/UWLCM_plotters/include/PlotterCommon.hpp +++ b/UWLCM_plotters/include/PlotterCommon.hpp @@ -30,10 +30,10 @@ class PlotterCommon { if(h5f.getFileName() != file) { - notice_macro("about to close current file: " << h5f.getFileName()) + notice_macro("Teraz about to close current file: " << h5f.getFileName()) h5f.close(); - notice_macro("about to open file: " << file) + notice_macro("Teraaz 1about to open file: " << file) h5f.openFile(file, H5F_ACC_RDONLY); } @@ -46,10 +46,10 @@ class PlotterCommon { if(h5f.getFileName() != file) { - notice_macro("about to close current file: " << h5f.getFileName()) + notice_macro("Teraz 2 about to close current file: " << h5f.getFileName()) h5f.close(); - notice_macro("about to open file: " << file) + notice_macro("Teraz 2about to open file: " << file) h5f.openFile(file, H5F_ACC_RDONLY); } @@ -85,7 +85,7 @@ class PlotterCommon file(file) { // init h5f - notice_macro("about to open file: " << file << "/const.h5") + notice_macro("Teraz 3about to open file: " << file << "/const.h5") h5f.openFile(file + "/const.h5", H5F_ACC_RDONLY); // init dt and outfreq diff --git a/drawbicyc/drawbicyc.cpp b/drawbicyc/drawbicyc.cpp index f1d9bfb..f46735e 100644 --- a/drawbicyc/drawbicyc.cpp +++ b/drawbicyc/drawbicyc.cpp @@ -37,7 +37,7 @@ int main(int argc, char** argv) throw std::runtime_error("Unrecognized 'type' option, only dycoms, rico, moist_thermal, pi_chamber, pi_chamber_icmw, cumulus_congestus, gccn_ccn_conc available now");//, base_prflux_vs_clhght available now"); // should profiles be normalized by inversion height - const bool normalize_prof = type == "dycoms"; + //const bool normalize_prof = type == "dycoms"; // parse dir name std::string @@ -56,8 +56,10 @@ int main(int argc, char** argv) H5::DataSet h5d = h5f.openDataSet("G"); H5::DataSpace h5s = h5d.getSpace(); int NDims = h5s.getSimpleExtentNdims(); - // detecting if subgrid model was on + std::cout << NDims <(h5, micro), plots, type); - if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), plots, type, normalize_prof); +// if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), plots, type, normalize_prof); //if(flag_fields) plot_fields(PlotterMicro_t<2>(h5, micro), plots, type); //if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro)); } - else if(NDims == 3) - { - if(flag_series) plot_series(PlotterMicro_t<3>(h5, micro), plots, type); +// else if(NDims == 3) +// { +// if(flag_series) plot_series(PlotterMicro_t<3>(h5, micro), plots, type); //if(flag_profiles) plot_profiles(PlotterMicro_t<3>(h5, micro), plots, type, normalize_prof); //if(flag_fields) plot_fields(PlotterMicro_t<3>(h5, micro), plots, type); //if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro)); @@ -87,7 +90,7 @@ int main(int argc, char** argv) // plot_lgrngn_spec_positions(PlotterMicro_t<3>(h5, "lgrngn")); // plot_lgrngn_spec(PlotterMicro_t<3>(h5, "lgrngn")); // } - } +// } else assert(false && "need 2d or 3d input data"); diff --git a/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp b/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp index 0431945..637a626 100644 --- a/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp +++ b/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp @@ -9,7 +9,7 @@ const std::vector series_Cumulus_Congestus({ "lwm", "cwm", "rwm", - //"surf_precip", + "surf_precip", "acc_precip", "acc_vol_precip", //"cl_nc", @@ -17,6 +17,7 @@ const std::vector series_Cumulus_Congestus({ //"RH_max", "cloud_top_height", "total_droplets_number", + "cloud_cover_rico", //"surf_flux_latent", //"surf_flux_sensible", //"sd_conc_avg", diff --git a/drawbicyc/include/gnuplot_series_set_labels.hpp b/drawbicyc/include/gnuplot_series_set_labels.hpp index b258c9d..49757f1 100644 --- a/drawbicyc/include/gnuplot_series_set_labels.hpp +++ b/drawbicyc/include/gnuplot_series_set_labels.hpp @@ -187,6 +187,24 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt) gp << "set xlabel ''\n"; gp << "set ylabel ''\n"; } + else if (plt == "cloud_top_height") + { + gp << "set title 'cloud top height [m]'\n"; + gp << "set xlabel 'time [min]'\n"; + gp << "set ylabel 'height [m]'\n"; + } + else if (plt == "total_droplets_number") + { + gp << "set title 'Total droplets number [#]'\n"; + gp << "set xlabel 'time [min]'\n"; + gp << "set ylabel ''\n"; + } + else if (plt == "acc_vol_precip") + { + gp << "set title 'accumulated volume precipitation [mm^3]'\n"; + gp << "set xlabel 'time [h]'\n"; + gp << "set ylabel ''\n"; + } else if (plt == "mass_dry") gp << "set title 'total dry mass [g]'\n"; else if (plt == "RH_max") @@ -198,13 +216,35 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt) else if (plt == "lwp") { gp << "set title 'liquid water path [g / m^2]'\n"; - gp << "set xlabel ''\n"; + gp << "set xlabel 'time[h]'\n"; + gp << "set ylabel ''\n"; + } + else if (plt == "lwm") + { + gp << "set title 'liquid water mass [kg]'\n"; + gp << "set xlabel 'time[h]'\n"; + gp << "set ylabel ''\n"; + }else if (plt == "cwm") + { + gp << "set title 'cloud water mass [kg]'\n"; + gp << "set xlabel 'time[h]'\n"; + gp << "set ylabel ''\n"; + }else if (plt == "rwm") + { + gp << "set title 'rain water mass [kg]'\n"; + gp << "set xlabel 'time[h]'\n"; gp << "set ylabel ''\n"; } else if (plt == "rwp") { gp << "set title 'rain water path [g / m^2]'\n"; - gp << "set xlabel ''\n"; + gp << "set xlabel 'time [h]'\n"; + gp << "set ylabel ''\n"; + } + else if (plt == "cwp") + { + gp << "set title 'cloud water path [g / m^2]'\n"; + gp << "set xlabel 'time [h]'\n"; gp << "set ylabel ''\n"; } else if (plt == "cloud_base_dycoms") diff --git a/drawbicyc/include/plot_series.hpp b/drawbicyc/include/plot_series.hpp index 7dd583d..66abcfd 100644 --- a/drawbicyc/include/plot_series.hpp +++ b/drawbicyc/include/plot_series.hpp @@ -28,7 +28,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) typename Plotter_t::arr_t rtot(rhod.shape()); typename Plotter_t::arr_t res_tmp(rhod.shape()); - + // for calculating running averages of u and w, needed in TKE calc in Pi chamber LES std::vector prev_u_vec, prev_w_vec; // container for the running sum From 347538c40a418413031a1486882ab8b31f6fc9ee Mon Sep 17 00:00:00 2001 From: ryba183 Date: Mon, 14 Aug 2023 23:17:47 +0200 Subject: [PATCH 5/5] prepare to pc test --- UWLCM_plotters/include/Plotter2d.hpp | 11 ++ UWLCM_plotters/include/PlotterCommon.hpp | 10 +- UWLCM_plotters/include/PlotterMicro.hpp | 12 +- drawbicyc/drawbicyc.cpp | 4 +- .../include/cases/Cumulus_Congestus/plots.hpp | 25 ++- .../include/gnuplot_profs_set_labels.hpp | 26 +++- drawbicyc/include/plot_prof.hpp | 142 +++++++++++++++++- drawbicyc/include/plot_series.hpp | 7 +- 8 files changed, 213 insertions(+), 24 deletions(-) diff --git a/UWLCM_plotters/include/Plotter2d.hpp b/UWLCM_plotters/include/Plotter2d.hpp index 2bed638..63ea185 100644 --- a/UWLCM_plotters/include/Plotter2d.hpp +++ b/UWLCM_plotters/include/Plotter2d.hpp @@ -198,6 +198,17 @@ class Plotter_t<2> : public PlotterCommon // edge cells are smaller dv(blitz::Range(0,0),blitz::Range::all()) /= 2.; dv(blitz::Range::all(),blitz::Range(0,0)) /= 2.; + dv(blitz::Range::all(),blitz::Range(dv.cols() - 1, dv.cols() - 1)) /= 2.; + dv(blitz::Range(dv.rows() - 1, dv.rows() - 1), blitz::Range::all()) /= 2.; + + // Printing the modified dv array + //std::cout << "Modified dv array:\n"; + //for (int i = 0; i < dv.extent(0); ++i) { + // for (int j = 0; j < dv.extent(1); ++j) { + // std::cout << dv(i, j) << " "; + // } + // std::cout << "\n"; + //} // other dataset are of the size x*z, resize tmp tmp.resize(n[0]-1, n[1]-1); diff --git a/UWLCM_plotters/include/PlotterCommon.hpp b/UWLCM_plotters/include/PlotterCommon.hpp index cc34d00..9d33f4b 100644 --- a/UWLCM_plotters/include/PlotterCommon.hpp +++ b/UWLCM_plotters/include/PlotterCommon.hpp @@ -30,10 +30,10 @@ class PlotterCommon { if(h5f.getFileName() != file) { - notice_macro("Teraz about to close current file: " << h5f.getFileName()) + //notice_macro("Teraz about to close current file: " << h5f.getFileName()) h5f.close(); - notice_macro("Teraaz 1about to open file: " << file) + //notice_macro("Teraaz 1about to open file: " << file) h5f.openFile(file, H5F_ACC_RDONLY); } @@ -46,10 +46,10 @@ class PlotterCommon { if(h5f.getFileName() != file) { - notice_macro("Teraz 2 about to close current file: " << h5f.getFileName()) + //notice_macro("Teraz 2 about to close current file: " << h5f.getFileName()) h5f.close(); - notice_macro("Teraz 2about to open file: " << file) + //notice_macro("Teraz 2about to open file: " << file) h5f.openFile(file, H5F_ACC_RDONLY); } @@ -85,7 +85,7 @@ class PlotterCommon file(file) { // init h5f - notice_macro("Teraz 3about to open file: " << file << "/const.h5") + //notice_macro("Teraz 3about to open file: " << file << "/const.h5") h5f.openFile(file + "/const.h5", H5F_ACC_RDONLY); // init dt and outfreq diff --git a/UWLCM_plotters/include/PlotterMicro.hpp b/UWLCM_plotters/include/PlotterMicro.hpp index 68dd71a..1da8f28 100644 --- a/UWLCM_plotters/include/PlotterMicro.hpp +++ b/UWLCM_plotters/include/PlotterMicro.hpp @@ -106,10 +106,16 @@ class PlotterMicro_t : public Plotter_t { if(this->micro == "lgrngn") { + + //std::cout << "CellVol = " << this->CellVol << std::endl; + //std::cout << "dv = " << this->dv << std::endl; + //std::cerr << this->dv; res = this->h5load_timestep("precip_rate", at) - * 4./3 * 3.14 * 1e3 // to get mass - / this->CellVol // averaged over cell volume, TODO: make precip rate return specific moment? wouldnt need the dx and dy - * L_evap; + * 4./3 * 3.14 * 1e3 // to get mass + //5000. + / this->dv // averaged over cell volume, TODO: make precip rate return specific moment? wouldnt need the dx and dy + // this->CellVol // averaged over cell volume, TODO: make precip rate return specific moment? wouldnt need the dx and dy + * L_evap; } else if(this->micro == "blk_1m") try diff --git a/drawbicyc/drawbicyc.cpp b/drawbicyc/drawbicyc.cpp index f46735e..a12306f 100644 --- a/drawbicyc/drawbicyc.cpp +++ b/drawbicyc/drawbicyc.cpp @@ -37,7 +37,7 @@ int main(int argc, char** argv) throw std::runtime_error("Unrecognized 'type' option, only dycoms, rico, moist_thermal, pi_chamber, pi_chamber_icmw, cumulus_congestus, gccn_ccn_conc available now");//, base_prflux_vs_clhght available now"); // should profiles be normalized by inversion height - //const bool normalize_prof = type == "dycoms"; + const bool normalize_prof = type == "dycoms"; // parse dir name std::string @@ -76,7 +76,7 @@ int main(int argc, char** argv) { if(flag_series) plot_series(PlotterMicro_t<2>(h5, micro), plots, type); -// if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), plots, type, normalize_prof); + if(flag_profiles) plot_profiles(PlotterMicro_t<2>(h5, micro), plots, type, normalize_prof); //if(flag_fields) plot_fields(PlotterMicro_t<2>(h5, micro), plots, type); //if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMicro_t<2>(h5, micro)); } diff --git a/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp b/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp index 637a626..544e360 100644 --- a/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp +++ b/drawbicyc/include/cases/Cumulus_Congestus/plots.hpp @@ -17,6 +17,14 @@ const std::vector series_Cumulus_Congestus({ //"RH_max", "cloud_top_height", "total_droplets_number", + "cl_meanr", + "cl_nc", + "cl_nr", + "nc", + "nr", + "com_mom0", + "com_mom1", + "com_mom2", "cloud_cover_rico", //"surf_flux_latent", //"surf_flux_sensible", @@ -34,20 +42,29 @@ const std::vector series_Cumulus_Congestus({ }); std::vector profs_Cumulus_Congestus({ - "00rtot", + "actrw_rw_cl_conc", "rliq", + "actrw_rw_cl", + "actrw_mom1", + "actrw_mom2", + "disp_r", + "actrw_rw_cl_sigma", + "sigma_r", + "mean_r", + //"00rtot", + //"rliq", //"thl", //"wvar", "prflux", //"clfrac", //"sd_conc", - "cl_nc", + //"cl_nc", //"cl_nc_up", //"w", //"u", //"v", - "coal_tele", - "base_prflux_vs_clhght", + //"coal_tele", + //"base_prflux_vs_clhght", //"non_gccn_rw_cl", //"gccn_rw_cl," //, "N_c", diff --git a/drawbicyc/include/gnuplot_profs_set_labels.hpp b/drawbicyc/include/gnuplot_profs_set_labels.hpp index b6448b4..eb83cdf 100644 --- a/drawbicyc/include/gnuplot_profs_set_labels.hpp +++ b/drawbicyc/include/gnuplot_profs_set_labels.hpp @@ -13,7 +13,7 @@ void gnuplot_profs_set_labels(Gnuplot &gp, std::string plt, const bool normalize if (plt == "rliq") { - gp << "set title 'liquid water [g/kg]'\n"; + gp << "set title 'liquid water in cloudy cells [g/kg]'\n"; } if (plt == "gccn_rw") { @@ -218,4 +218,28 @@ void gnuplot_profs_set_labels(Gnuplot &gp, std::string plt, const bool normalize { gp << "set title 'mean cloud droplet radius [um]'\n"; } + else if (plt == "actrw_mom1") + { + gp << "set title 'mean radius of actrw droplets [um]'\n"; + } + else if (plt == "actrw_mom2") + { + gp << "set title 'variance of actrw droplet size distribution [um^2] * 1e3'\n"; + } + else if (plt == "actrw_sd") + { + gp << "set title 'standard deviation of actrw droplet size distribution [um]'\n"; + } + else if (plt == "actrw_rw_cl_conc") + { + gp << "set title 'concentration of actrw droplet in cloudy cells [cm^-3]'\n"; + } + else if (plt == "actrw_rw_cl") + { + gp << "set title 'mean radius of actrw droplet in cloudy cells [um]'\n"; + } + else if (plt == "disp_r") + { + gp << "set title 'relative dispersion of r_{drop} [1]'\n"; + } } diff --git a/drawbicyc/include/plot_prof.hpp b/drawbicyc/include/plot_prof.hpp index 14cce55..389eabd 100644 --- a/drawbicyc/include/plot_prof.hpp +++ b/drawbicyc/include/plot_prof.hpp @@ -121,8 +121,27 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; // rain typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); res_tmp = iscloudy_rc_rico(snap); - res_prof_sum = plotter.horizontal_sum(res_tmp); - res_prof_hlpr = where(prof_tmp > 0, plotter.horizontal_sum(res) / res_prof_sum, 0); + res *= res_tmp; + prof_tmp = plotter.horizontal_sum(res_tmp); + res_prof_hlpr = where(prof_tmp > 0, plotter.horizontal_sum(res) / prof_tmp, 0); + } + if (plt == "actrw_rw_cl_conc") + { + // concentration of activated droplets (r > rc) in cloudy cells + { + typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); + res_tmp = iscloudy_rc_rico(snap); + } + { + auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]); + typename Plotter_t::arr_t snap(tmp); + res_tmp2 = snap; + res_tmp2 *= rhod / 1e6; // per cm^3 + } + // cloudy only + prof_tmp = plotter.horizontal_sum(res_tmp); + res_tmp2 *= res_tmp; + res_prof_hlpr = where(prof_tmp > 0 , plotter.horizontal_sum(res_tmp2) / prof_tmp, 0); } /* if (plt == "cloud_water_std") { @@ -485,7 +504,29 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool { auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); - res_tmp = where(res_tmp > 0 , res_tmp / snap, res_tmp); + res_tmp = where(snap > 0 , res_tmp / snap, 0); + } + { + typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); + res_tmp2 = iscloudy_rc_rico(snap); + res_tmp *= res_tmp2; + } + // mean only over downdraught cells + prof_tmp = plotter.horizontal_sum(res_tmp2); // number of downdraft cells on a given level + res_prof_hlpr = where(prof_tmp > 0 , plotter.horizontal_sum(res_tmp) / prof_tmp, 0); + } + if (plt == "actrw_rw_cl_sigma") + { + // sigma(radius) of activated droplets (r > rc) in cloudy cells + { + auto tmp = plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]) * 1e6; + typename Plotter_t::arr_t snap(tmp); + res_tmp = snap; + } + { + auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]); + typename Plotter_t::arr_t snap(tmp); + res_tmp = where(snap > 0 , res_tmp / snap, 0); } { typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); @@ -898,9 +939,14 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool { // precipitation flux(doesnt include vertical volicty w!) { - res = plotter.h5load_prflux_timestep(at * n["outfreq"]); + //res = plotter.dv; + //std::cerr << plotter.h5load_prflux_timestep(at * n["outfreq"]); + //res = (plotter.dv).shape(); + res = plotter.h5load_prflux_timestep(at * n["outfreq"]); // plotter.dv; res_prof_hlpr = plotter.horizontal_mean(res); // average in x - } + //res_prof_hlpr = plotter.horizontal_sum(res)/122.; // average in x + + } // add vertical velocity to precipitation flux (3rd mom of cloud drops * w) /* { @@ -1051,11 +1097,95 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res = where(res > 0 , res / res_tmp, res); res_prof_hlpr = plotter.horizontal_mean(res); // average in x } + if (plt == "actrw_mom1") + { + // mean radius of actrw droplets + res = plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]) * 1e6; + res_tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]); + res = where(res > 0 , res / res_tmp, res); + res_prof_hlpr = plotter.horizontal_mean(res); // average in x + } + else if (plt == "actrw_mom2") + { + // variance of actrw droplet size distribution + res = plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]) * 1e6 * 1e3; // convert to (um)^2 * 1e3 + res_tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]); + res = where(res > 0 , res / res_tmp, res); + res_prof_hlpr = plotter.horizontal_mean(res); // average in x + } + else if (plt == "disp_r") + { + // relative dispersion (std dev / mean) of droplet radius distribution averaged over cells away from walls + //typename Plotter_t::arr_t m0(plotter.nowall(typename Plotter_t::arr_t(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])), distance_from_walls)); + typename Plotter_t::arr_t m0(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"])); + typename Plotter_t::arr_t m1(plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"])); + typename Plotter_t::arr_t m2(plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"])); + // calculate stddev of radius, store in m2 + m2 = where(m0 > 0, + m2 / m0 - m1 / m0 * m1 / m0 , 0.); + // might be slightly negative due to numerical errors + m2 = where(m2 < 0, 0, m2); + m2 = sqrt(m2); // sqrt(variance) -// ====================================================================================== + //calculate mean radius, store in m1 + m1 = where(m0 > 0, + m1 / m0, 0.); + + res = where(m1 > 0 , m2 / m1 , 0); + typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); + res_tmp = iscloudy_rc_rico(snap); + res *= res_tmp; + prof_tmp = plotter.horizontal_sum(res_tmp); + + res_prof_hlpr = where(prof_tmp > 0, plotter.horizontal_sum(res) / prof_tmp, 0); + } + else if (plt == "sigma_r") + { + // std dev of droplet radius distribution averaged over cells + //typename Plotter_t::arr_t m0(plotter.nowall(typename Plotter_t::arr_t(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])), distance_from_walls)); + typename Plotter_t::arr_t m0(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"])); + typename Plotter_t::arr_t m1(plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"])); + typename Plotter_t::arr_t m2(plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"])); + + // calculate stddev of radius, store in m2 + m2 = where(m0 > 0, + m2 / m0 - m1 / m0 * m1 / m0 , 0.); + + // might be slightly negative due to numerical errors + m2 = where(m2 < 0, 0, m2); + res = sqrt(m2); // sqrt(variance) + + typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); + res_tmp = iscloudy_rc_rico(snap); + res *= res_tmp; + prof_tmp = plotter.horizontal_sum(res_tmp); + + res_prof_hlpr = where(prof_tmp > 0, plotter.horizontal_sum(res) / prof_tmp, 0); + } + else if (plt == "mean_r") + { + // relative dispersion (std dev / mean) of droplet radius distribution averaged over cells away from walls + //typename Plotter_t::arr_t m0(plotter.nowall(typename Plotter_t::arr_t(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])), distance_from_walls)); + typename Plotter_t::arr_t m0(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"])); + typename Plotter_t::arr_t m1(plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"])); + + //calculate mean radius, store in m1 + res = where(m0 > 0, + m1 / m0, 0.); + + typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); + res_tmp = iscloudy_rc_rico(snap); + res *= res_tmp; + prof_tmp = plotter.horizontal_sum(res_tmp); + + res_prof_hlpr = where(prof_tmp > 0, plotter.horizontal_sum(res) / prof_tmp, 0); + } + +// ====================================================================================== + if(normalize) { // interpolate profile to uniform vertical grid (height normalized by inversion height) diff --git a/drawbicyc/include/plot_series.hpp b/drawbicyc/include/plot_series.hpp index 66abcfd..b1f536e 100644 --- a/drawbicyc/include/plot_series.hpp +++ b/drawbicyc/include/plot_series.hpp @@ -1300,15 +1300,16 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) else if (plt == "cl_meanr") { // cloud droplets mean radius in cloudy grid cells + // modified for actrw_mom1. cloud -> actrw try { // cloud fraction (cloudy if N_c > 20/cm^3) - typename Plotter_t::arr_t snap(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); + typename Plotter_t::arr_t snap(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"])); snap *= rhod; // b4 it was specific moment snap /= 1e6; // per cm^3 snap = iscloudy(snap); // cloudiness mask - typename Plotter_t::arr_t snap_m0(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - typename Plotter_t::arr_t snap_m1(plotter.h5load_timestep("cloud_rw_mom1", at * n["outfreq"])*1e6); // in microns + typename Plotter_t::arr_t snap_m0(plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"])); + typename Plotter_t::arr_t snap_m1(plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"])*1e6); // in microns snap_m0 *= snap; snap_m1 *= snap; auto tot_m0 = blitz::sum(snap_m0);