diff --git a/Energy_spectrum/spectrum_refined.py b/Energy_spectrum/spectrum_refined.py index 13eeff7..526fa9a 100644 --- a/Energy_spectrum/spectrum_refined.py +++ b/Energy_spectrum/spectrum_refined.py @@ -34,6 +34,7 @@ lmbd = {} level_start_idx = {} level_end_idx = {} + exists = {} # read some constant parameters with h5py.File(directory + "/const.h5", 'r') as consth5: @@ -58,16 +59,21 @@ # initiliaze nx,ny,nz,dx and average energy for each variable for var in args.vars: filename = directory + "/timestep" + str(time_start_idx).zfill(10) + ".h5" - w3d = h5py.File(filename, "r")[var][:,:,:] - nx[var], ny[var], nz[var] = tuple(x for x in w3d.shape) - dx[var] = X / (nx[var] - 1) - dz[var] = Z / (nz[var] - 1) - Exy_avg[var] = np.zeros(int((nx[var]-1)/2 + 1)) - ref[var] = int(dx_adve / dx[var]) - assert(float(args.level_start / dz[var]).is_integer()) - assert(float(args.level_end / dz[var]).is_integer()) - level_start_idx[var] = int(args.level_start / dz[var]) - level_end_idx[var] = int(args.level_end / dz[var]) + 1 + + try: + w3d = h5py.File(filename, "r")[var][:,:,:] + nx[var], ny[var], nz[var] = tuple(x for x in w3d.shape) + dx[var] = X / (nx[var] - 1) + dz[var] = Z / (nz[var] - 1) + Exy_avg[var] = np.zeros(int((nx[var]-1)/2 + 1)) + ref[var] = int(dx_adve / dx[var]) + assert(float(args.level_start / dz[var]).is_integer()) + assert(float(args.level_end / dz[var]).is_integer()) + level_start_idx[var] = int(args.level_start / dz[var]) + level_end_idx[var] = int(args.level_end / dz[var]) + 1 + exists[var] = True + except: + exists[var] = False # time loop @@ -78,6 +84,8 @@ # variables loop for var in args.vars: print(var) + if not exists[var]: + continue print(nx[var],dx[var][0]) w3d = h5py.File(filename, "r")[var][0:nx[var]-1,0:ny[var]-1,:] # * 4. / 3. * 3.1416 * 1e3 @@ -85,7 +93,7 @@ for lvl in range(level_start_idx[var], level_end_idx[var]): w2d = w3d[:, :, lvl] #print w2d - + wkx = 1.0 / np.sqrt(nx[var] - 1) * np.fft.rfft(w2d, axis = 0) wky = 1.0 / np.sqrt(ny[var] - 1) * np.fft.rfft(w2d, axis = 1) @@ -107,15 +115,21 @@ # print(K, lmbd) if (t == time_start_idx and lab==args.labels[0]): - plt.loglog(lmbd[var], 2e-6* K**(-5./3.) ) + L = np.array([2e2, 2e3]) + plt.loglog(L, 2e-5 * L**(5./3.), label = "-5/3" , color="black", ls='dotted') + L = np.array([5e1, 3e2]) + plt.loglog(L, 2e-8 * L**(3.), label = "-3" , color="black", ls='dashed') for var in args.vars: + if not exists[var]: + continue + Exy_avg[var] /= (time_end_idx - time_start_idx) / outfreq + 1 Exy_avg[var] /= level_end_idx[var] - level_start_idx[var] #crudely scale #Exy_avg[var] /= Exy_avg[var][len(Exy_avg[var])-1] - Exy_avg[var] /= np.sum(Exy_avg[var]) +# Exy_avg[var] /= np.sum(Exy_avg[var]) #plot plt.loglog(lmbd[var], Exy_avg[var] , linewidth=2, label=lab+"_"+var) diff --git a/Matplotlib_common/latex_labels.py b/Matplotlib_common/latex_labels.py index 20c7265..1ee0471 100644 --- a/Matplotlib_common/latex_labels.py +++ b/Matplotlib_common/latex_labels.py @@ -31,16 +31,14 @@ "surf_precip" : 'Surf. precip. [mm/day]', "acc_precip" : 'Acc. precip. [mm]', "cl_nc" : '$N_c$ [cm$^{-3}$]', - "cl_nc_rico" : '$N_c$ [cm$^{-3}$]', # "cl_nr" : '$N_{r>25\mu m}$ [cm$^{-3}$]', "cl_nr" : '$N_{r}$ [cm$^{-3}$]', - "cl_nr_rico" : '$N_{r}$ [cm$^{-3}$]', "cl_gccn_conc" : '$N_{GCCN}$ [cm$^{-3}$]', "thl" : r'$\theta_l$ [K]', "00rtot" : '$q_{t}$ [g/kg]', "rliq" : '$q_{l}$ [g/kg]', "clfrac" : 'Cloud fraction', - "cloud_cover_rico" : 'Cloud cover', + "cloud_cover" : 'Cloud cover', "cloud_cover_dycoms" : 'Cloud cover', "prflux" : 'Precip. flux [W m$^{-2}$]', "wvar" : r'Var$\left(w\right)$ [m$^2$ s$^{-2}$]', @@ -53,7 +51,7 @@ "er" : 'Entrain. rate [cm s$^{-1}$]', "wvarmax" : 'Max. $w$ var. [m$^{2}$ s$^{-2}$]', "cloud_base_dycoms" : 'Cloud base [m]', - "min_cloud_base_rico" : 'Lowest cloud base [m]', + "min_cloud_base" : 'Lowest cloud base [m]', "inversion_height_rico" : 'Inversion height [m]', "gccn_rw_cl" : '$$ of GCCN droplets [um]', "non_gccn_rw_cl" : '$$ of CCN droplets [um]', @@ -62,13 +60,18 @@ "clb_bigrain_mean_conc" : 'conc. of (r$>$40um) @ clbase [1/cc]', "clb_bigrain_mean_inclt" : 'time since act. of (r$>$40um) @ clbase [s]', "clb_bigrain_mean_gccn_fraction" : 'frac. of (r$>$40um) on gccn @ clbase', - "cl_acnv25_rico" : 'acnv. rate [g m$^{-3}$ d$^{-1}$]', - "cl_accr25_rico" : 'accr. rate [g m$^{-3}$ d$^{-1}$]', - "cl_acnv25_dycoms" : 'acnv. rate [g m$^{-3}$ d$^{-1}$]', - "cl_accr25_dycoms" : 'accr. rate [g m$^{-3}$ d$^{-1}$]', - "cloud_avg_supersat" : '$$ in cloud [\%]', - "cl_meanr" : '$$ in cloud [um]', - "sd_conc" : '$$' + "cl_acnv25" : 'acnv. rate [g m$^{-3}$ d$^{-1}$]', + "cl_accr25" : 'accr. rate [g m$^{-3}$ d$^{-1}$]', + "cl_avg_supersat" : '$$ in cloud [\%]', + "cl_avg_th" : '$<\\theta>$ in cloud [K]', + "cl_avg_rv" : '$$ in cloud [1]', + "cl_avg_cloud_meanr" : '$$ of cloud droplets in cloud [um]', + "cl_avg_cloud_stddevr" : '$\sigma(r)$ of cloud droplets in cloud [um]', + "cl_avg_act_meanr" : '$$ of activated droplets in cloud [um]', + "cl_avg_act_stddevr" : '$\sigma(r)$ of activated droplets in cloud [um]', + "sd_conc" : '$$', + "surf_flux_latent" : 'latent surface flux [W/m$^2$]', + "surf_flux_sensible" : 'sensible surface flux [W/m$^2$]', } @@ -93,4 +96,10 @@ 17 : "(r)", 18 : "(s)", 19 : "(t)", + 20 : "(u)", + 21 : "(v)", + 22 : "(w)", + 23 : "(x)", + 24 : "(y)", + 25 : "(z)", } diff --git a/Matplotlib_common/plot_series.py b/Matplotlib_common/plot_series.py index 8194773..1f79bbe 100644 --- a/Matplotlib_common/plot_series.py +++ b/Matplotlib_common/plot_series.py @@ -1,17 +1,31 @@ import numpy as np -from sys import argv from math import floor +#from sys import argv +import argparse +import matplotlib.pyplot as plt from latex_labels import var_labels from read_UWLCM_arrays import * + def plot_series(var_list, plot_iter, nplotx, nploty, axarr, xscaledict, yscaledict, xlimdict, ylimdict, show_bin=False, suffix='', xlabel='', ylabeldict=var_labels, file_names=[], file_labels=[], linewidth=1): + + parser = argparse.ArgumentParser(description='Plot UWLCM series comparison') + parser.add_argument("-ts", "--time_start", type=float, required=False, help="start of the plotted period [s] (override default)") + parser.add_argument("-te", "--time_end", type=float, required=False, help="end of the plotted period [s] (override default)") + parser.add_argument("-d", "--dirs", action="extend", nargs="+", type=str, help="list of directories with the data", required=True) + parser.add_argument("-l", "--labels", action="extend", nargs="+", type=str, help="list of labels of the data (same order as --dirs)", required=True) + args, extra = parser.parse_known_args() + + if args.time_start is not None and args.time_end is not None: + xlimdict = {x: (args.time_start, args.time_end) for x in xlimdict} + + # if file names are not defined, read them and labels from command line if len(file_names)==0: - file_no = np.arange(1, len(argv)-1 , 2) - for no in file_no: - file_names.append(argv[no] + suffix) - file_labels.append(argv[no+1]) + for directory, lab in zip(args.dirs, args.labels): + file_names.append(directory + suffix) + file_labels.append(lab) for var in var_list: label_counter=0 @@ -20,6 +34,11 @@ def plot_series(var_list, plot_iter, nplotx, nploty, axarr, xscaledict, yscaledi series_file = open(file_name, "r") my_times = read_my_var(series_file, "position") my_res = read_my_var(series_file, var) + if len(my_res) == 0: # file does not contain this type of plot + print("skipping from: " + str(label_counter)) + label_counter+=1 + print("skipping to: " + str(label_counter)) + continue # rescale time to hours my_times = my_times / 3600. @@ -32,14 +51,15 @@ def plot_series(var_list, plot_iter, nplotx, nploty, axarr, xscaledict, yscaledi linestyles = ['--', '-.', ':'] dashList = [(3,1),(1,1),(4,1,1,1),(4,2)] - colorList = ['red', 'blue', 'green'] + # colorList = ['red', 'blue', 'green'] + colorList = plt.rcParams['axes.prop_cycle'].by_key()['color'] # default prop cycle colors # x label only on he lowest row xlabel_used = xlabel if plot_iter < nploty: xlabel_used = '' - plot_my_array(axarr, plot_iter, my_times, my_res, nploty, xlabel=xlabel_used, ylabel=ylabeldict[var], varlabel=file_labels[label_counter], dashes = dashList[label_counter % len(dashList)], xscale=xscaledict[var], yscale=yscaledict[var], xlim=xlimdict[var], ylim=ylimdict[var], linewidth=linewidth)#, color = colorList[int(floor(label_counter / len(dashList)))]) + plot_my_array(axarr, plot_iter, my_times, my_res, nploty, xlabel=xlabel_used, ylabel=ylabeldict[var], varlabel=file_labels[label_counter], dashes = dashList[label_counter % len(dashList)], xscale=xscaledict[var], yscale=yscaledict[var], xlim=xlimdict[var], ylim=ylimdict[var], linewidth=linewidth, color = colorList[label_counter % len(colorList)]) label_counter+=1 plot_iter = plot_iter + 1 return plot_iter diff --git a/Matplotlib_common/read_UWLCM_arrays.py b/Matplotlib_common/read_UWLCM_arrays.py index 82b4c3a..1441f55 100644 --- a/Matplotlib_common/read_UWLCM_arrays.py +++ b/Matplotlib_common/read_UWLCM_arrays.py @@ -18,6 +18,8 @@ def read_my_var(file_obj, var_name): arr_name = file_obj.readline() if(str(arr_name).strip() == str(var_name).strip()): return read_my_array(file_obj) + elif(len(arr_name)==0): + return np.empty(0) #def plot_my_array(axarr, plot_iter, time, val, nploty, xlabel=None, ylabel=None, varlabel=None , linestyle='--', dashes=(5,2), xlim=None, ylim=None, xscale="linear", yscale="linear"): # x = int(plot_iter / nploty) diff --git a/README.md b/README.md index 5ac1cdd..2d34b4d 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,28 @@ # UWLCM_plotting a set of scripts for plotting output of the UWLCM model + +CONTENTS + +plotting: + +- /drawbicyc - C++ program for processing UWLCM output; can calculate and plot series, profiles and snapshots for multiple parameters; also has C++ programs for averaging and comparing series or profiles from multiple runs, although this is probably better done in Python (e.g. /cases) + +- /cases - Python scripts for comparing series and profiles calculated using drawbicyc; specific to each LES case + +- /Energy_spectrum - Python scripts for plotting and comparing Fourier spectra of variables; works directly on UWLCM output + +- /histograms - Python scripts for plotting and comparing histograms with spatial and/or temporal distribution of variables; also capable of plotting scatter plots of correlation between two variables; works directly on UWLCM output + +- /NC_vs_AF - Python scripts for plotting scatter plots of number concentration of cloud droplets vs adiabatic fraction, or vs liquid water content; works directly on UWLCM output; NOTE: redundant to scatter plots in /histograms? + +- /papers - Python scripts for making figures used in our papers + +- /Size_spectra - Python scripts for comparing droplet size distributions; works directly on UWLCM output + + +helpers: + +- /Matplotlib_common - Python scripts for reading drawbicyc output; also plotting comparison of series and profiles from drawbicyc + +- /UWLCM_plotters - C++ classes for reading UWLCM output; used in drawbicyc + diff --git a/UWLCM_plotters/include/Plotter2d.hpp b/UWLCM_plotters/include/Plotter2d.hpp index 2bed638..f7ed60c 100644 --- a/UWLCM_plotters/include/Plotter2d.hpp +++ b/UWLCM_plotters/include/Plotter2d.hpp @@ -1,13 +1,13 @@ #pragma once #include "common.hpp" -#include "PlotterCommon.hpp" +#include "PlotterH5.hpp" template -class Plotter_t : public PlotterCommon {}; +class Plotter_t : public PlotterH5 {}; // 2d version template<> -class Plotter_t<2> : public PlotterCommon +class Plotter_t<2> : public PlotterH5 { public: static const int n_dims = 2; @@ -18,10 +18,10 @@ class Plotter_t<2> : public PlotterCommon arr_t dv; protected: - using parent_t = PlotterCommon; + using parent_t = PlotterH5; hsize_t n[2]; enum {x, z}; - arr_t tmp, tmp_srfc; + arr_t tmp, tmp_srfc, tmp_ref; public: @@ -202,6 +202,7 @@ class Plotter_t<2> : public PlotterCommon // other dataset are of the size x*z, resize tmp tmp.resize(n[0]-1, n[1]-1); tmp_srfc.resize(n[0]-1, 1); + tmp_ref.resize(tmp.shape()); } }; diff --git a/UWLCM_plotters/include/Plotter3d.hpp b/UWLCM_plotters/include/Plotter3d.hpp index 9d80fed..f62bcc2 100644 --- a/UWLCM_plotters/include/Plotter3d.hpp +++ b/UWLCM_plotters/include/Plotter3d.hpp @@ -4,7 +4,7 @@ // 3d version template<> -class Plotter_t<3> : public PlotterCommon +class Plotter_t<3> : public PlotterH5 { public: static const int n_dims = 3; @@ -15,10 +15,10 @@ class Plotter_t<3> : public PlotterCommon arr_t dv; protected: - using parent_t = PlotterCommon; + using parent_t = PlotterH5; hsize_t n[3]; enum {x, y, z}; - arr_t tmp, tmp_srfc; + arr_t tmp, tmp_srfc, tmp_ref; blitz::Range yrange; public: @@ -36,15 +36,17 @@ class Plotter_t<3> : public PlotterCommon cnt[3] = { n[x], n[y], srfc ? 1 : n[z] }, off[3] = { 0, 0, 0 }; this->h5s.selectHyperslab( H5S_SELECT_SET, cnt, off); + + arr_t &arr(tmp.extent(0) == n[x] ? tmp : tmp_ref); // crude check if it is refined or normal data; assume surface data is never refined hsize_t ext[3] = { - hsize_t(tmp.extent(0)), - hsize_t(tmp.extent(1)), - hsize_t(srfc ? tmp_srfc.extent(2) : tmp.extent(2)) + hsize_t(arr.extent(0)), + hsize_t(arr.extent(1)), + hsize_t(srfc ? tmp_srfc.extent(2) : arr.extent(2)) }; - this->h5d.read(srfc ? tmp_srfc.data() : tmp.data(), H5::PredType::NATIVE_FLOAT, H5::DataSpace(3, ext), h5s); + this->h5d.read(srfc ? tmp_srfc.data() : arr.data(), H5::PredType::NATIVE_FLOAT, H5::DataSpace(3, ext), h5s); - return blitz::safeToReturn((srfc ? tmp_srfc : tmp) + 0); + return blitz::safeToReturn((srfc ? tmp_srfc : arr) + 0); } auto h5load_timestep( @@ -235,6 +237,44 @@ class Plotter_t<3> : public PlotterCommon // other dataset are of the size x*z, resize tmp tmp.resize(n[0]-1, n[1]-1, n[2]-1); tmp_srfc.resize(n[0]-1, n[1]-1, 1); + + // init refined data + try + { + this->h5f.openDataSet("X refined").getSpace().getSimpleExtentDims(n, NULL); + this->map["refined x"] = n[0]-1; + this->map["refined y"] = n[1]-1; + this->map["refined z"] = n[2]-1; + tmp_ref.resize(n[0], n[1], n[2]); + h5load(file + "/const.h5", "X refined"); + this->map["refined dx"] = tmp_ref(1,0,0) - tmp_ref(0,0,0); + h5load(file + "/const.h5", "Y refined"); + this->map["refined dy"] = tmp_ref(0,1,0) - tmp_ref(0,0,0); + h5load(file + "/const.h5", "Z refined"); + this->map["refined dz"] = tmp_ref(0,0,1) - tmp_ref(0,0,0); + this->CellVol_ref = this->map["refined dx"] * this->map["refined dy"] * this->map["refined dz"]; + tmp_ref.resize(n[0]-1, n[1]-1, n[2]-1); + } + catch(...) // for pre-refinement simulations that didnt store refined stuff + { + this->map["refined x"] = this->map["x"]; + this->map["refined y"] = this->map["y"]; + this->map["refined z"] = this->map["z"]; + this->map["refined dx"] = this->map["dx"]; + this->map["refined dy"] = this->map["dy"]; + this->map["refined dz"] = this->map["dz"]; + this->CellVol_ref = this->CellVol; + tmp_ref.resize(tmp.shape()); + } + + for (auto const& x : this->map) +{ + std::cout << x.first // string (key) + << ':' + << x.second // string's value + << std::endl; +} + } }; diff --git a/UWLCM_plotters/include/PlotterCommon.hpp b/UWLCM_plotters/include/PlotterCommon.hpp index 1c75c83..77495ce 100644 --- a/UWLCM_plotters/include/PlotterCommon.hpp +++ b/UWLCM_plotters/include/PlotterCommon.hpp @@ -1,126 +1,269 @@ #pragma once +#include "Plotter2d.hpp" +#include "Plotter3d.hpp" +//#include -#include -#include -#include -#include - -class PlotterCommon +template +class PlotterCommon : public Plotter_t { - public: - using arr_prof_t = blitz::Array; + protected: + using parent_t = Plotter_t; - const string file; - std::map map; - std::map map_prof; - blitz::Array timesteps; - double CellVol, DomainSurf, DomainVol; + public: + using arr_t = typename parent_t::arr_t; + using parent_t::parent_t; protected: - H5::H5File h5f; - H5::DataSet h5d; - H5::Group h5g; - H5::DataSpace h5s; - - void h5load( - const string &file, - const string &dataset, - bool srfc = false - ) + const double L_evap = 2264.76e3; // latent heat of evaporation [J/kg] + + public: + + // ---- functions for diagnosing statistics ---- + // + arr_t& multiply_by_rhod(arr_t arr) { - if(h5f.getFileName() != file) - { - notice_macro("about to close current file: " << h5f.getFileName()) - h5f.close(); + if(arr.extent(NDims-1) == this->map_prof["rhod"].extent(0)) + arr *= this->map_prof["rhod"](this->LastIndex); + else if(arr.extent(NDims-1) == this->map_prof["refined rhod"].extent(0)) + arr *= this->map_prof["refined rhod"](this->LastIndex); + else + throw std::runtime_error("multiply_by_rhod: input array is neither normal grid size nor refined grid size"); + return arr; + } - notice_macro("about to open file: " << file) - h5f.openFile(file, H5F_ACC_RDONLY); + arr_t& multiply_by_CellVol(arr_t arr) + { + if(arr.extent(NDims-1) == this->map_prof["rhod"].extent(0)) + arr *= this->CellVol; + else if(arr.extent(NDims-1) == this->map_prof["refined rhod"].extent(0)) + arr *= this->CellVol_ref; + else + throw std::runtime_error("multiply_by_CellVol: input array is neither normal grid size nor refined grid size"); + return arr; + } + + // helper function that calculates staistics (mean and std_dev) of a field + std::pair hlpr(arr_t arr, int at) + { + std::pair res; + res.first = blitz::mean(arr); + arr = pow(arr - res.first, 2); + res.second = sqrt(blitz::mean(arr)); + return res; + } + + // mean and std_dev of supersaturation in cells with positive supersaturation [%] (characteristics of the spatial distribution at this timestep) + std::pair positive_supersat_stats_timestep(int at) + { + if(this->micro == "blk_1m") return {0,0}; + std::pair res; + + // read RH + arr_t RH(this->h5load_timestep("RH", at)); + RH -= 1.; + arr_t sat(RH.copy()); + sat = iscloudy_sat(RH); + RH *= sat; //apply the cloudiness mask + RH *= 100; // to get % + if(blitz::sum(sat) > 0) + res.first = blitz::sum(RH) / blitz::sum(sat); + else + res.first = 0; + + RH = pow(RH - res.first, 2); + RH *= sat; // apply filter + if(res.first>0) + res.second = sqrt(blitz::sum(RH) / blitz::sum(sat)); + else + res.second = 0.; + + return res; + } + + // mean and std_dev of supersaturation at droplet locations (i.e. supersat weighted by the number of droplets) + std::pair drop_supersat_stats_timestep(int at) + { + if(this->micro == "blk_1m") return {0,0}; + std::pair res; + + // read supersat + arr_t ssat(this->h5load_timestep("RH", at)); + ssat -= 1.; + + // read number of droplets per cell + arr_t nc(this->h5load_timestep("cloud_rw_mom0", at)); + multiply_by_rhod(nc); + multiply_by_CellVol(nc); + + + if(blitz::sum(nc) > 0) + { + res.first = blitz::sum(ssat * nc) / blitz::sum(nc); + ssat = pow(ssat - res.first, 2); + res.second = sqrt(blitz::sum(ssat * nc) / blitz::sum(nc)); } + else + { + res.first = 0; + res.second = 0; + } + + return res; + } + + // mean and std_dev of number of SDs (characteristics of the spatial distribution at this timestep) + std::pair sdconc_stats_timestep(int at) + { + if(this->micro == "blk_1m") return {0,0}; + arr_t sdconc(this->h5load_timestep("sd_conc", at)); + return hlpr(sdconc, at); + } + + // mean and std_dev of temperature [K] (characteristics of the spatial distribution at this timestep, without near-wall cells) + std::pair T_stats_timestep(int at) + { + // read theta away from walls + //arr_t tht(this->nowall(arr_t(this->h5load_timestep("th", at)), distance_from_walls)); + arr_t tht(arr_t(this->h5load_timestep("th", at))); + tht *= pow(this->map_prof["p_e"](this->LastIndex) / p_1000, R_d / c_pd); // tht -> T + return hlpr(tht, at); + } - notice_macro("about to read dataset: " << dataset) - h5d = h5f.openDataSet(dataset); - h5s = h5d.getSpace(); + // mean and std_dev of r_v [1] (characteristics of the spatial distribution at this timestep) + std::pair rv_stats_timestep(int at) + { + arr_t rv(arr_t(this->h5load_timestep("rv", at))); + return hlpr(rv, at); } - float h5load_attr(const string &file, const string &attr_name, const string &group_name) + // mean and std_dev of RH [1] (characteristics of the spatial distribution at this timestep) + std::pair RH_stats_timestep(int at) + { + if(this->micro == "blk_1m") return {0,0}; + std::pair res; + + arr_t RH(this->h5load_timestep("RH", at)); + return hlpr(RH, at); + } + + // surface precipitation since last output [mm/day] + double calc_surf_precip(double prec_vol_diff) + { + if(this->micro == "lgrngn") + return prec_vol_diff / this->DomainSurf / (double(this->map["outfreq"]) * this->map["dt"] / 3600. / 24.) * 1e3; // SDM + if(this->micro == "blk_1m" || this->micro == "blk_2m") + return prec_vol_diff / double(this->map["outfreq"]) // flux in [kg / m^3 / s] averaged over time since last output and over cells on the bottom + / (this->map["x"] * this->map["y"]) + * 3600. * 24. // per day + * this->map["dz"] // per m^2 + / 1e3 // to m^3 of water + * 1e3; // to mm + } + + // accumulated surface precipitation [mm] + double calc_acc_surf_precip(double prec_vol) + { + if(this->micro == "lgrngn") + return prec_vol / this->DomainSurf * 1e3; + if(this->micro == "blk_1m" || this->micro == "blk_2m") + return prec_vol * this->map["dt"] + / (this->map["x"] * this->map["y"]) + * this->map["dz"] // per m^2 + / 1e3 // to m^3 of water + * 1e3; // to mm + } + + // droplet removal rate (at boundaries) since last output [1/(cm^3 s)] + double calc_prtcl_removal(double prtcl_removal_diff) { - if(h5f.getFileName() != file) + if(this->micro == "lgrngn") + return prtcl_removal_diff / this->DomainVol / (double(this->map["outfreq"]) * this->map["dt"]) / 1e6; + if(this->micro == "blk_1m") + return 0; + } + + // heat flux thru boundary since last output [W/m^2] + double calc_heat_flux(double tot_th_diff, int z_idx) // input in [K] + { + if(this->micro == "lgrngn") { - notice_macro("about to close current file: " << h5f.getFileName()) - h5f.close(); - - notice_macro("about to open file: " << file) - h5f.openFile(file, H5F_ACC_RDONLY); + tot_th_diff *= pow(this->map_prof["p_e"](z_idx) / p_1000, R_d / c_pd); // tht -> T + double ret = tot_th_diff * c_pd * this->map_prof["rhod"](z_idx) // sum of th diff over boundary cells since last output (K) * c_pd * density + * this->map["dz"] / ((this->map["x"]-1) * (this->map["y"]-1)) // multiply by cell volume and divide by domain surface area (without walls) + * (double(this->map["outfreq"]) * this->map["dt"]); // divide by time since last output + return ret; } - - notice_macro(std::string("about to read group: " + group_name)) - h5g = h5f.openGroup(group_name); - - float ret; - notice_macro(std::string("about to open attribute: " + attr_name)) - auto attr = h5g.openAttribute(attr_name); - notice_macro(std::string("about to read attribute value")) - attr.read(attr.getDataType(), &ret); - return ret; + if(this->micro == "blk_1m") + return 0; } - template - void plot(gp_t &gp) + double calc_heat_flux_top(double mean_th_diff, bool errfix) { - //gp << "set cbtics format \"%.2tE%+03T\"\n"; - gp << "set cbtics font \", 8\"\n"; - // gp << "set rmargin 2cm\n"; + double tot_th_diff = mean_th_diff * this->map["x"] * this->map["y"]; + if(errfix) + { + // th_diff += (this->map["x"] * this->map["y"] - 1) * 280; // to counter to error in tot_th_diff calculation in UWLCM + tot_th_diff -= (2*this->map["x"] + 2*(this->map["y"] - 1)) * (280 - 285); // dont count side wall + } + return calc_heat_flux(tot_th_diff, this->map["z"]-1); } - public: + double calc_heat_flux_bot(double mean_th_diff, bool errfix) + { + double tot_th_diff = mean_th_diff * this->map["x"] * this->map["y"]; + if(errfix) + { + // th_diff += (this->map["x"] * this->map["y"] - 1) * 299; + tot_th_diff -= (2*this->map["x"] + 2*(this->map["y"] - 1)) * (299 - 285); + } + return calc_heat_flux(tot_th_diff, 0); + } - float h5load_attr_timestep(int at, const std::string attr_name, const std::string group_name = "/") + // kinematic rv flux thru boundary since last output [kg/kg * m / s] + double calc_moist_flux(double tot_rv_diff) // rv change summed over horizontal plane [kg/kg] { - string timestep_file = file + "/timestep" + zeropad(at, 10) + ".h5"; - return h5load_attr(timestep_file, attr_name, group_name); + if(this->micro == "lgrngn") + { + return tot_rv_diff * this->map["dz"] + * (double(this->map["outfreq"]) * this->map["dt"]); // divide by time since last output + } + if(this->micro == "blk_1m") + return 0; } - //ctor - PlotterCommon(const string &file): - file(file) + double calc_moist_flux_top(double mean_rv_diff, bool errfix) { - // init h5f - notice_macro("about to open file: " << file << "/const.h5") - h5f.openFile(file + "/const.h5", H5F_ACC_RDONLY); + // 3D assumed here! + double tot_rv_diff = mean_rv_diff * this->map["x"] * this->map["y"]; + if(errfix) + { +// rv_diff += (this->map["x"] * this->map["y"] - 1) * 0.0062192674278; // to counter to error in tot_rv_diff calculation in UWLCM + tot_rv_diff -= (2*this->map["x"] + 2*(this->map["y"] - 1)) * (0.0062192674278 - 0.00611718803008); // dont count side wall + } + return calc_moist_flux(tot_rv_diff); + } - // init dt and outfreq + double calc_moist_flux_bot(double mean_rv_diff, bool errfix) + { + // 3D assumed here! + double tot_rv_diff = mean_rv_diff * this->map["x"] * this->map["y"]; + if(errfix) { - map["dt"] = h5load_attr(file + "/const.h5", "dt", "advection"); - map["outfreq"] = h5load_attr(file + "/const.h5", "outfreq", "user_params"); - map["MPI_compiler"] = h5load_attr(file + "/const.h5", "MPI compiler (true/false)", "MPI details"); - - // read number of timesteps - hsize_t n; - h5load(file + "/const.h5", "T"); - h5s.getSimpleExtentDims(&n, NULL); - this->map["t"] = n; - // read timesteps - timesteps.resize(n); - h5d.read(timesteps.data(), H5::PredType::NATIVE_FLOAT); - - // read environmental pressure profile - h5load(file + "/const.h5", "p_e"); - h5s.getSimpleExtentDims(&n, NULL); - map_prof.emplace("p_e", arr_prof_t(n)); - h5d.read(map_prof["p_e"].data(), H5::PredType::NATIVE_FLOAT); - - // read SGS mixing length profile - h5load(file + "/const.h5", "mix_len"); - h5s.getSimpleExtentDims(&n, NULL); - map_prof.emplace("mix_len", arr_prof_t(n)); - h5d.read(map_prof["mix_len"].data(), H5::PredType::NATIVE_FLOAT); - - // read dry air density profile - h5load(file + "/const.h5", "rhod"); - h5s.getSimpleExtentDims(&n, NULL); - map_prof.emplace("rhod", arr_prof_t(n)); - h5d.read(map_prof["rhod"].data(), H5::PredType::NATIVE_FLOAT); +// rv_diff += (this->map["x"] * this->map["y"] - 1) * 0.0213489271007; // to counter to error in tot_rv_diff calculation in UWLCM + tot_rv_diff -= (2*this->map["x"] + 2*(this->map["y"] - 1)) * (0.0213489271007 - 0.00611718803008); // dont count side wall } + return calc_moist_flux(tot_rv_diff); + } + + //ctor + /* + PlotterCommon(const string &file, const string µ): + parent_t(file), + micro(micro), + res(this->tmp.shape()), + rhod(this->h5load(file + "/const.h5", "G")) + { } + */ }; diff --git a/UWLCM_plotters/include/PlotterH5.hpp b/UWLCM_plotters/include/PlotterH5.hpp new file mode 100644 index 0000000..ed62860 --- /dev/null +++ b/UWLCM_plotters/include/PlotterH5.hpp @@ -0,0 +1,140 @@ +#pragma once + +#include +#include +#include +#include + +class PlotterH5 +{ + public: + using arr_prof_t = blitz::Array; + + const string file; + std::map map; + std::map map_prof; + blitz::Array timesteps; + double CellVol, DomainSurf, DomainVol, CellVol_ref; + + protected: + H5::H5File h5f; + H5::DataSet h5d; + H5::Group h5g; + H5::DataSpace h5s; + + void h5load( + const string &file, + const string &dataset, + bool srfc = false + ) + { + if(h5f.getFileName() != file) + { + notice_macro("about to close current file: " << h5f.getFileName()) + h5f.close(); + + notice_macro("about to open file: " << file) + h5f.openFile(file, H5F_ACC_RDONLY); + } + + notice_macro("about to read dataset: " << dataset) + h5d = h5f.openDataSet(dataset); + h5s = h5d.getSpace(); + } + + float h5load_attr(const string &file, const string &attr_name, const string &group_name) + { + if(h5f.getFileName() != file) + { + notice_macro("about to close current file: " << h5f.getFileName()) + h5f.close(); + + notice_macro("about to open file: " << file) + h5f.openFile(file, H5F_ACC_RDONLY); + } + + notice_macro(std::string("about to read group: " + group_name)) + h5g = h5f.openGroup(group_name); + + float ret; + notice_macro(std::string("about to open attribute: " + attr_name)) + auto attr = h5g.openAttribute(attr_name); + notice_macro(std::string("about to read attribute value")) + attr.read(attr.getDataType(), &ret); + notice_macro(std::string("attribute value read: ") + std::to_string(ret)) + return ret; + } + + template + void plot(gp_t &gp) + { + //gp << "set cbtics format \"%.2tE%+03T\"\n"; + gp << "set cbtics font \", 8\"\n"; + // gp << "set rmargin 2cm\n"; + } + + public: + + float h5load_attr_timestep(int at, const std::string attr_name, const std::string group_name = "/") + { + string timestep_file = file + "/timestep" + zeropad(at, 10) + ".h5"; + return h5load_attr(timestep_file, attr_name, group_name); + } + + //ctor + PlotterH5(const string &file): + file(file) + { + // init h5f + notice_macro("about to open file: " << file << "/const.h5") + h5f.openFile(file + "/const.h5", H5F_ACC_RDONLY); + + // init dt and outfreq + { + map["dt"] = h5load_attr(file + "/const.h5", "dt", "advection"); + map["outfreq"] = h5load_attr(file + "/const.h5", "outfreq", "user_params"); + map["MPI_compiler"] = h5load_attr(file + "/const.h5", "MPI compiler (true/false)", "MPI details"); + + // read number of timesteps + hsize_t n; + h5load(file + "/const.h5", "T"); + h5s.getSimpleExtentDims(&n, NULL); + this->map["t"] = n; + // read timesteps + timesteps.resize(n); + h5d.read(timesteps.data(), H5::PredType::NATIVE_FLOAT); + + // read environmental pressure profile + h5load(file + "/const.h5", "p_e"); + h5s.getSimpleExtentDims(&n, NULL); + map_prof.emplace("p_e", arr_prof_t(n)); + h5d.read(map_prof["p_e"].data(), H5::PredType::NATIVE_FLOAT); + + // read SGS mixing length profile + h5load(file + "/const.h5", "mix_len"); + h5s.getSimpleExtentDims(&n, NULL); + map_prof.emplace("mix_len", arr_prof_t(n)); + h5d.read(map_prof["mix_len"].data(), H5::PredType::NATIVE_FLOAT); + + // read dry air density profile + h5load(file + "/const.h5", "rhod"); + h5s.getSimpleExtentDims(&n, NULL); + map_prof.emplace("rhod", arr_prof_t(n)); + h5d.read(map_prof["rhod"].data(), H5::PredType::NATIVE_FLOAT); + + // read dry air density profile on refined grid + try + { + h5load(file + "/const.h5", "refined rhod"); + h5s.getSimpleExtentDims(&n, NULL); + map_prof.emplace("refined rhod", arr_prof_t(n)); + h5d.read(map_prof["refined rhod"].data(), H5::PredType::NATIVE_FLOAT); + } + catch(...) // for pre-refinement simulations, use rhod as refined rhod + { + map_prof.emplace("refined rhod", map_prof["rhod"].copy()); + } + } + } +}; + diff --git a/UWLCM_plotters/include/PlotterMask.hpp b/UWLCM_plotters/include/PlotterMask.hpp new file mode 100644 index 0000000..4029c42 --- /dev/null +++ b/UWLCM_plotters/include/PlotterMask.hpp @@ -0,0 +1,255 @@ +#pragma once +#include "PlotterMicro.hpp" + +enum class mask_type_t{Rico11, Dycoms_rf02, unset}; + +template +class PlotterMask : public PlotterMicro +{ + public: + using parent_t = PlotterMicro; + using arr_t = typename parent_t::arr_t; + + private: + + mask_type_t mask_type; + arr_t mask; + + void calc_mask(int at) + { + if(mask_type == mask_type_t::Rico11) + { + mask = this->load_rc_timestep(at); + mask = iscloudy_rc_rico(mask); + } + else if(mask_type == mask_type_t::Dycoms_rf02) + { + mask = this->load_nc_timestep(at); + mask = iscloudy_nc_dycoms(mask); + } + else throw std::runtime_error("Invalid mask type"); + mask(this->hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux + } + + //lgrngn_droplet_prefix + + private: + // helper function that calculates staistics (mean and std_dev) of a field only in cloudy cells + std::pair cloud_hlpr(arr_t arr, int at) + { + std::pair res; + + // apply mask + calc_mask(at); + arr *= mask; + + if(blitz::sum(mask) > 0.) + res.first = blitz::sum(arr) / blitz::sum(mask); + else + res.first = 0.; + + arr = pow(arr - res.first, 2); + arr *= mask; // apply filter + if(res.first>0) + res.second = sqrt(blitz::sum(arr) / blitz::sum(mask)); + else + res.second = 0.; + + return res; + } + + + std::pair cloud_conc_stats_timestep_hlpr(int at, std::string lgrngn_prefix, std::string blk_2m_1, std::string blk_2m_2 = "") + { + arr_t conc; + // read concentration of activated droplets + if(this->micro == "blk_1m") return {0,0}; + // TODO: fix stupid copying of arrays + else if(this->micro == "lgrngn") { + arr_t tmp(this->h5load_timestep(lgrngn_prefix+"_mom0", at)); + conc.resize(tmp.shape()); + conc = tmp; + } + else if(this->micro == "blk_2m") { + arr_t tmp(arr_t(this->h5load_timestep(blk_2m_1, at))); + if(blk_2m_2 != "") + tmp += arr_t(this->h5load_timestep(blk_2m_2, at)); + conc.resize(tmp.shape()); + conc = tmp; + } + this->multiply_by_rhod(conc); // b4 it was specific moment + conc /= 1e6; // per cm^3 + return cloud_hlpr(conc, at); + } + + std::pair cloud_meanr_stats_timestep_hlpr(int at, std::string lgrngn_prefix) + { + if(this->micro == "blk_1m" || this->micro == "blk_2m") return {0,0}; + // read drop 0th raw moment / mass [1/kg] + arr_t m0th(this->h5load_timestep(lgrngn_prefix+"_mom0", at)); + // read drop 1st raw moment / mass [um/kg] + arr_t m1st(this->h5load_timestep(lgrngn_prefix+"_mom1", at) * 1e6); + // calculate mean radius, store in 1st + m1st = where(m0th > 0, m1st / m0th, 0.); + return cloud_hlpr(m1st, at); + } + + std::pair cloud_stddevr_stats_timestep_hlpr(int at, std::string lgrngn_prefix) + { + if(this->micro == "blk_1m" || this->micro == "blk_2m") return {0,0}; + // read drop 0th raw moment / mass [1/kg] + arr_t m0th(this->h5load_timestep(lgrngn_prefix+"_mom0", at)); + // read drop 1st raw moment / mass [um/kg] + arr_t m1st(this->h5load_timestep(lgrngn_prefix+"_mom1", at) * 1e6); + // read drop 2nd raw moment / mass [um^2/kg] + arr_t m2nd(this->h5load_timestep(lgrngn_prefix+"_mom2", at) * 1e12); + // calculate stddev of radius, store in 1st + m1st = where(m0th > 0, + m2nd / m0th - m1st / m0th * m1st / m0th, 0.); + // might be slightly negative due to numerical errors + m1st = where(m1st < 0, 0, m1st); + m1st = sqrt(m1st); + return cloud_hlpr(m1st, at); + } + + public: + + arr_t get_mask(int at) + { + calc_mask(at); + return mask; + } + + // functions for diagnosing statistics + // mean and std dev [g/kg] of the mixing ratio of activated dropelts in cloudy cells (characteristics of the spatial distribution at this timestep) + std::pair cloud_ract_stats_timestep(int at) + { + // read activated droplets mixing ratio + arr_t ract(this->load_ract_timestep(at)); + ract *= 1e3; // turn it into g/kg + return cloud_hlpr(ract, at); + } + + // same for supersat + std::pair cloud_supersat_stats_timestep(int at) + { + // read activated droplets mixing ratio + arr_t RH(this->load_RH_timestep(at)); + RH -= 1.; + RH *= 100; + return cloud_hlpr(RH, at); + } + // th + std::pair cloud_theta_stats_timestep(int at) + { + // read activated droplets mixing ratio + try + { + arr_t th(this->h5load_timestep("refined th", at)); + return cloud_hlpr(th, at); + } + catch (...) // should work for simulations without refinement + { + arr_t th(this->h5load_timestep("th", at)); + return cloud_hlpr(th, at); + } + } + // rv + std::pair cloud_rv_stats_timestep(int at) + { + // read activated droplets mixing ratio + try + { + arr_t rv(this->h5load_timestep("refined rv", at)); + return cloud_hlpr(rv, at); + } + catch (...) // should work for simulations without refinement + { + arr_t rv(this->h5load_timestep("rv", at)); + return cloud_hlpr(rv, at); + } + } + + // mean and std_dev of concentration of activated/cloud/rain droplets in cloudy cells [1/cm^3] (characteristics of the spatial distribution at this timestep) + std::pair cloud_actconc_stats_timestep(int at) + { + return cloud_conc_stats_timestep_hlpr(at, "actrw_rw", "nc", "nr"); + } + + std::pair cloud_cloudconc_stats_timestep(int at) + { + return cloud_conc_stats_timestep_hlpr(at, "cloud_rw", "nc"); + } + + std::pair cloud_rainconc_stats_timestep(int at) + { + return cloud_conc_stats_timestep_hlpr(at, "rain_rw", "nr"); + } + + + // mean and std_dev of number of SDs in cloudy cells (characteristics of the spatial distribution at this timestep) + std::pair cloud_sdconc_stats_timestep(int at) + { + if(this->micro == "blk_1m" || this->micro == "blk_2m") return {0,0}; + arr_t sdconc(this->h5load_timestep("sd_conc", at)); + return cloud_hlpr(sdconc, at); + } + + // mean and std_dev of number of activated SDs in cloudy cells (characteristics of the spatial distribution at this timestep) + std::pair cloud_sdconc_act_stats_timestep(int at) + { + if(this->micro == "blk_1m" || this->micro == "blk_2m") return {0,0}; + arr_t sdconc_act(this->h5load_timestep("sd_conc_act", at)); + return cloud_hlpr(sdconc_act, at); + } + + // mean and std_dev of mean radius of activated/cloud/rain droplets in cloudy cells [um] (characteristics of the spatial distribution at this timestep) + std::pair cloud_actmeanr_stats_timestep(int at) + { + return cloud_meanr_stats_timestep_hlpr(at, "actrw_rw"); + } + std::pair cloud_cloudmeanr_stats_timestep(int at) + { + return cloud_meanr_stats_timestep_hlpr(at, "cloud_rw"); + } + std::pair cloud_rainmeanr_stats_timestep(int at) + { + return cloud_meanr_stats_timestep_hlpr(at, "rain_rw"); + } + + // mean and std_dev of std_dev of radius of activated/cloud/rain droplets in cloudy cells [um] (characteristics of the spatial distribution at this timestep) + std::pair cloud_actstddevr_stats_timestep(int at) + { + return cloud_stddevr_stats_timestep_hlpr(at, "actrw_rw"); + } + std::pair cloud_cloudstddevr_stats_timestep(int at) + { + return cloud_stddevr_stats_timestep_hlpr(at, "cloud_rw"); + } + std::pair cloud_rainstddevr_stats_timestep(int at) + { + return cloud_stddevr_stats_timestep_hlpr(at, "rain_rw"); + } + + // height [m] of the center of mass of activated droplets + double act_com_z_timestep(int at) + { + arr_t ract(this->load_ract_timestep(at)); + arr_t weighted(ract.copy()); + weighted = weighted * this->LastIndex * this->map["dz"]; + if(blitz::sum(ract) > 1e-3) + return blitz::sum(weighted) / blitz::sum(ract); + else + return 0.; + } + + + //ctor + PlotterMask(const string &file, const string µ, const mask_type_t _mt): + parent_t(file, micro), + mask(this->tmp_ref.shape()), + mask_type(_mt) + { + } +}; + diff --git a/UWLCM_plotters/include/PlotterMicro.hpp b/UWLCM_plotters/include/PlotterMicro.hpp index 68dd71a..697ceb5 100644 --- a/UWLCM_plotters/include/PlotterMicro.hpp +++ b/UWLCM_plotters/include/PlotterMicro.hpp @@ -1,14 +1,12 @@ #pragma once -#include "Plotter2d.hpp" -#include "Plotter3d.hpp" -#include +#include "PlotterCommon.hpp" // TODO: make two: plotterlgrngn and plotter blk1m template -class PlotterMicro_t : public Plotter_t +class PlotterMicro : public PlotterCommon { protected: - using parent_t = Plotter_t; + using parent_t = PlotterCommon; public: using arr_t = typename parent_t::arr_t; @@ -20,57 +18,62 @@ class PlotterMicro_t : public Plotter_t const double L_evap = 2264.76e3; // latent heat of evaporation [J/kg] public: + // functions for diagnosing fields // // aerosol droplets mixing ratio - auto h5load_ra_timestep( + auto load_ra_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) //-> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") - res = this->h5load_timestep("aerosol_rw_mom3", at) * 4./3. * 3.1416 * 1e3; + return arr_t(this->h5load_timestep("aerosol_rw_mom3", at) * 4./3. * 3.1416 * 1e3); + else if(this->micro == "blk_1m") + { res = 0; - return blitz::safeToReturn(res + 0); + return res; + // return blitz::safeToReturn(res + 0); + } } // cloud droplets mixing ratio - auto h5load_rc_timestep( + auto load_rc_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) //-> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") - res = this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3; + return arr_t(this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3); else if(this->micro == "blk_1m") - res = this->h5load_timestep("rc", at); + return arr_t(this->h5load_timestep("rc", at)); else if(this->micro == "blk_2m") - res = this->h5load_timestep("rc", at); - return blitz::safeToReturn(res + 0); + return arr_t(this->h5load_timestep("rc", at)); } // rain droplets mixing ratio - auto h5load_rr_timestep( + auto load_rr_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) //-> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") - res = this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3; + return arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3); else if(this->micro == "blk_1m") - res = this->h5load_timestep("rr", at); + return arr_t(this->h5load_timestep("rr", at)); else if(this->micro == "blk_2m") - res = this->h5load_timestep("rr", at); - return blitz::safeToReturn(res + 0); + return arr_t(this->h5load_timestep("rr", at)); } // activated drops mixing ratio - auto h5load_ract_timestep( + auto load_ract_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) { if(this->micro == "lgrngn") { - res = this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3; - res += arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3); + return arr_t( + arr_t(this->h5load_timestep("cloud_rw_mom3", at) * 4./3. * 3.1416 * 1e3) + + arr_t(this->h5load_timestep("rain_rw_mom3", at) * 4./3. * 3.1416 * 1e3) + ); } else if(this->micro == "blk_1m") { @@ -82,143 +85,72 @@ class PlotterMicro_t : public Plotter_t res = this->h5load_timestep("rc", at); res += arr_t(this->h5load_timestep("rr", at)); } - return blitz::safeToReturn(res + 0); + // return blitz::safeToReturn(res + 0); + return res; } - // cloud droplets concentration [1/kg] - auto h5load_nc_timestep( + // cloud droplets concentration [1/cm^3] + auto load_nc_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) //-> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") - res = this->h5load_timestep("cloud_rw_mom0", at); + return arr_t(this->multiply_by_rhod(arr_t(this->h5load_timestep("cloud_rw_mom0", at) / 1e6))); else if(this->micro == "blk_1m") + { res = 0; + return res; + // return blitz::safeToReturn(res + 0); + } else if(this->micro == "blk_2m") - res = this->h5load_timestep("nc", at); - return blitz::safeToReturn(res + 0); + return arr_t(this->multiply_by_rhod(arr_t(this->h5load_timestep("nc", at) / 1e6))); } // precipitation flux [W/m2] - auto h5load_prflux_timestep( + auto load_prflux_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + )// -> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") { - res = this->h5load_timestep("precip_rate", at) + return arr_t(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; + * L_evap); } else if(this->micro == "blk_1m") try { res = this->h5load_timestep("precip_rate", at); // precip_rate is the difference between influx and outflux + // uppermost cell + res(this->hrzntl_slice(this->map["z"] - 1)) *= rhod(this->hrzntl_slice(this->map["z"] - 1)); + // cells below for(int z = this->map["z"] - 2; z>=0; --z) { - res(this->hrzntl_slice(z)) = res(this->hrzntl_slice(z+1)) - res(this->hrzntl_slice(z)); + res(this->hrzntl_slice(z)) = res(this->hrzntl_slice(z+1)) + res(this->hrzntl_slice(z)) * rhod(this->hrzntl_slice(z)); } - res *= rhod * this->map["dz"] * L_evap; + res *= -1 * this->map["dz"] * L_evap; } catch(...) { res = 0; } - return blitz::safeToReturn(res + 0); + // return blitz::safeToReturn(res + 0); + return res; } // RH - auto h5load_RH_timestep( + auto load_RH_timestep( int at - ) -> decltype(blitz::safeToReturn(arr_t() + 0)) + ) //-> decltype(blitz::safeToReturn(arr_t() + 0)) { if(this->micro == "lgrngn") - res = this->h5load_timestep("RH", at); + return arr_t(this->h5load_timestep("RH", at)); else if(this->micro == "blk_1m") res = 0; - return blitz::safeToReturn(res + 0); - } - - - // functions for diagnosing statistics - - // helper function that calculates staistics (mean and std_dev) of a field - std::pair hlpr(arr_t arr, int at) - { - std::pair res; - res.first = blitz::mean(arr); - arr = pow(arr - res.first, 2); - res.second = sqrt(blitz::mean(arr)); - return res; + // return blitz::safeToReturn(res + 0); + return res; } - - // helper function that calculates staistics (mean and std_dev) of a field only in cloudy cells - std::pair cloud_hlpr(arr_t arr, int at) - { - std::pair res; - // read activated droplets mixing ratio - arr_t mask(h5load_ract_timestep(at)); - mask = iscloudy_rc(mask); - arr *= mask; // apply filter - - if(blitz::sum(mask) > 0.) - res.first = blitz::sum(arr) / blitz::sum(mask); - else - res.first = 0.; - - arr = pow(arr - res.first, 2); - arr *= mask; // apply filter - if(res.first>0) - res.second = sqrt(blitz::sum(arr) / blitz::sum(mask)); - else - res.second = 0.; - - return res; - } - - // height [m] of the center of mass of activated droplets - double act_com_z_timestep(int at) - { - arr_t ract(h5load_ract_timestep(at)); - arr_t weighted(ract.copy()); - weighted = weighted * this->LastIndex * this->map["dz"]; - if(blitz::sum(ract) > 1e-3) - return blitz::sum(weighted) / blitz::sum(ract); - else - return 0.; - } - - // mean and std dev [g/kg] of the mixing ratio of activated dropelts in cloudy cells (characteristics of the spatial distribution at this timestep) - std::pair cloud_ract_stats_timestep(int at) - { - // read activated droplets mixing ratio - arr_t ract(h5load_ract_timestep(at)); - ract *= 1e3; // turn it into g/kg - return cloud_hlpr(ract, at); - } - - // mean and std_dev of concentration of activated droplets in cloudy cells [1/cm^3] (characteristics of the spatial distribution at this timestep) - std::pair cloud_actconc_stats_timestep(int at) - { - arr_t actconc; - // read concentration of activated droplets - if(this->micro == "blk_1m") return {0,0}; - // TODO: fix stupid copying of arrays - else if(this->micro == "lgrngn") { - arr_t tmp(this->h5load_timestep("actrw_rw_mom0", at)); - actconc.resize(tmp.shape()); - actconc = tmp; - } - else if(this->micro == "blk_2m") { - arr_t tmp(arr_t(this->h5load_timestep("nc", at)) + arr_t(this->h5load_timestep("nr", at))); - actconc.resize(tmp.shape()); - actconc = tmp; - } - actconc *= rhod; // b4 it was specific moment - actconc /= 1e6; // per cm^3 - return cloud_hlpr(actconc, at); - } // mean and std_dev of supersaturation in cells with positive supersaturation [%] (characteristics of the spatial distribution at this timestep) std::pair positive_supersat_stats_timestep(int at) @@ -284,55 +216,7 @@ class PlotterMicro_t : public Plotter_t { if(this->micro == "blk_1m") return {0,0}; arr_t sdconc(this->h5load_timestep("sd_conc", at)); - return hlpr(sdconc, at); - } - - // mean and std_dev of number of SDs in cloudy cells (characteristics of the spatial distribution at this timestep) - std::pair cloud_sdconc_stats_timestep(int at) - { - if(this->micro == "blk_1m") return {0,0}; - arr_t sdconc(this->h5load_timestep("sd_conc", at)); - return cloud_hlpr(sdconc, at); - } - - // mean and std_dev of number of activated SDs in cloudy cells (characteristics of the spatial distribution at this timestep) - std::pair cloud_sdconc_act_stats_timestep(int at) - { - if(this->micro == "blk_1m") return {0,0}; - arr_t sdconc_act(this->h5load_timestep("sd_conc_act", at)); - return cloud_hlpr(sdconc_act, at); - } - - // mean and std_dev of mean radius of activated droplets in cloudy cells [um] (characteristics of the spatial distribution at this timestep) - std::pair cloud_meanr_stats_timestep(int at) - { - if(this->micro == "blk_1m") return {0,0}; - // read act drop 0th raw moment / mass [1/kg] - arr_t act0th(this->h5load_timestep("actrw_rw_mom0", at)); - // read act drop 1st raw moment / mass [um/kg] - arr_t act1st(this->h5load_timestep("actrw_rw_mom1", at) * 1e6); - // calculate mean radius, store in act1st - act1st = where(act0th > 0, act1st / act0th, 0.); - return cloud_hlpr(act1st, at); - } - - // mean and std_dev of std_dev of radius of activated droplets in cloudy cells [um] (characteristics of the spatial distribution at this timestep) - std::pair cloud_stddevr_stats_timestep(int at) - { - if(this->micro == "blk_1m") return {0,0}; - // read act drop 0th raw moment / mass [1/kg] - arr_t act0th(this->h5load_timestep("actrw_rw_mom0", at)); - // read act drop 1st raw moment / mass [um/kg] - arr_t act1st(this->h5load_timestep("actrw_rw_mom1", at) * 1e6); - // read act drop 2nd raw moment / mass [um^2/kg] - arr_t act2nd(this->h5load_timestep("actrw_rw_mom2", at) * 1e12); - // calculate stddev of radius, store in act1st - act1st = where(act0th > 0, - act2nd / act0th - act1st / act0th * act1st / act0th, 0.); - // might be slightly negative due to numerical errors - act1st = where(act1st < 0, 0, act1st); - act1st = sqrt(act1st); - return cloud_hlpr(act1st, at); + return this->hlpr(sdconc, at); } // mean and std_dev of temperature [K] (characteristics of the spatial distribution at this timestep, without near-wall cells) @@ -342,14 +226,14 @@ class PlotterMicro_t : public Plotter_t //arr_t tht(this->nowall(arr_t(this->h5load_timestep("th", at)), distance_from_walls)); arr_t tht(arr_t(this->h5load_timestep("th", at))); tht *= pow(this->map_prof["p_e"](this->LastIndex) / p_1000, R_d / c_pd); // tht -> T - return hlpr(tht, at); + return this->hlpr(tht, at); } // mean and std_dev of r_v [1] (characteristics of the spatial distribution at this timestep) std::pair rv_stats_timestep(int at) { arr_t rv(arr_t(this->h5load_timestep("rv", at))); - return hlpr(rv, at); + return this->hlpr(rv, at); } // mean and std_dev of RH [1] (characteristics of the spatial distribution at this timestep) @@ -359,7 +243,7 @@ class PlotterMicro_t : public Plotter_t std::pair res; arr_t RH(this->h5load_timestep("RH", at)); - return hlpr(RH, at); + return this->hlpr(RH, at); } // surface precipitation since last output [mm/day] @@ -479,11 +363,11 @@ class PlotterMicro_t : public Plotter_t } //ctor - PlotterMicro_t(const string &file, const string µ): + PlotterMicro(const string &file, const string µ): parent_t(file), micro(micro), - res(this->tmp.shape()), - rhod(this->h5load(file + "/const.h5", "G")) + res(this->tmp.shape()) + //,rhod(this->h5load(file + "/const.h5", "G")) { } }; diff --git a/UWLCM_plotters/include/common_filters.hpp b/UWLCM_plotters/include/common_filters.hpp index 1ad0b2e..9b8dd04 100644 --- a/UWLCM_plotters/include/common_filters.hpp +++ b/UWLCM_plotters/include/common_filters.hpp @@ -26,11 +26,11 @@ double is_th_prtrb(double x) } BZ_DECLARE_FUNCTION(is_th_prtrb) -double iscloudy(double x) +double iscloudy_nc_dycoms(double x) { return x > 20. ? 1. : 0.; } -BZ_DECLARE_FUNCTION(iscloudy) +BZ_DECLARE_FUNCTION(iscloudy_nc_dycoms) double iscloudy_sat(double x) { diff --git a/cases/RICO11/Rico_series_comparison.py b/cases/RICO11/Rico_series_comparison.py index 3caf101..a604b01 100644 --- a/cases/RICO11/Rico_series_comparison.py +++ b/cases/RICO11/Rico_series_comparison.py @@ -2,6 +2,8 @@ import matplotlib.pyplot as plt import matplotlib.ticker as ticker from matplotlib.ticker import AutoMinorLocator, MultipleLocator +import argparse + import os import sys @@ -16,14 +18,20 @@ # activate latex text rendering rc('text', usetex=True) -rico_vars = ["lwp", "rwp", "cloud_cover_rico", "min_cloud_base_rico", "inversion_height_rico", "cl_nc_rico", "cl_nr_rico", "surf_precip", "acc_precip", "RH_max", "cloud_avg_supersat", "wvarmax", "cl_meanr", "sd_conc"]#, "cl_acnv25_rico", "cl_accr25_rico"]#, "surf_flux_latent", "surf_flux_sensible" ] +rico_vars = ["lwp", "rwp", "cloud_cover", "min_cloud_base", "inversion_height_rico", "cl_nc", "cl_nr", "surf_precip", "acc_precip", "RH_max", "wvarmax", "sd_conc", "cl_avg_act_meanr", "cl_avg_act_stddevr", "cl_avg_cloud_meanr", "cl_avg_cloud_stddevr", "cl_avg_supersat", "cl_avg_th", "cl_avg_rv", "surf_flux_latent", "surf_flux_sensible"]#, "cl_acnv25", "cl_accr25"]# ] #rico_vars = ["clb_bigrain_mean_rd","clb_bigrain_mean_kappa","clb_bigrain_mean_conc","clb_bigrain_mean_inclt", "cl_nr"] # variables that need rescaling of the yrange to the limited x range of 1-6h -rescale_vars = ["lwp", "cloud_cover_rico", "min_cloud_base_rico", "inversion_height_rico", "cl_nc"]# rico_vars +rescale_vars = ["lwp", "cloud_cover", "min_cloud_base", "inversion_height_rico", "cl_nc"]# rico_vars + +# arguments +parser = argparse.ArgumentParser(description='Plot UWLCM series comparison for RICO simulations.') +parser.add_argument("-of", "--outfig", help="output file name", required=True) +args, extra = parser.parse_known_args() + # init the plot -nplotx = 4 +nplotx = 6 nploty= 4 fig, axarr = plt.subplots(nplotx,nploty) x_arr = np.arange(nplotx) @@ -64,7 +72,7 @@ var = rico_vars[x*nploty + y] if var in rescale_vars: autoscale_y(axarr[x,y], margin=0.1) - if(var == "cloud_cover_rico" or var == "cl_nr"): + if(var == "cloud_cover" or var == "cl_nr"): axarr[x,y].yaxis.set_major_formatter(ticker.FormatStrFormatter('%0.2f')) # hide hrzntl tic labels @@ -78,11 +86,11 @@ 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) +fig.set_size_inches(7.874, 1.5 * nplotx + (len(labels) - 2) * 0.1)# 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.subplots_adjust(bottom=0.18 + (len(labels) - 2) * 0.01, 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,)) +fig.savefig(args.outfig, bbox_inches='tight', dpi=300)#, bbox_extra_artists=(lgd,)) diff --git a/cases/RICO11/plot_ranges.py b/cases/RICO11/plot_ranges.py index 2db01e2..ae530e5 100644 --- a/cases/RICO11/plot_ranges.py +++ b/cases/RICO11/plot_ranges.py @@ -5,8 +5,6 @@ "prflux" : "linear", "cl_nc" : "linear", "cl_nr" : "linear", - "cl_nc_rico" : "linear", - "cl_nr_rico" : "linear", "clfrac" : "linear", "wvar" : "linear", "w3rd" : "linear", @@ -17,6 +15,8 @@ "er" : "linear", "wvarmax" : "linear", "surf_precip" : "linear", + "surf_flux_latent" : "linear", + "surf_flux_sensible" : "linear", "acc_precip" : "linear", "cloud_base" : "linear", "gccn_rw_cl" : "linear", @@ -29,13 +29,19 @@ "clb_bigrain_mean_conc" : "linear", "clb_bigrain_mean_inclt" : "linear", "clb_bigrain_mean_gccn_fraction" : "linear", - "cloud_cover_rico" : "linear", - "min_cloud_base_rico" : "linear", + "cloud_cover" : "linear", + "min_cloud_base" : "linear", "inversion_height_rico" : "linear", - "cl_acnv25_rico" : "linear", - "cl_accr25_rico" : "linear", + "cl_acnv25" : "linear", + "cl_accr25" : "linear", "RH_max" : "linear", - "cloud_avg_supersat" : "linear", + "cl_avg_supersat" : "linear", + "cl_avg_th" : "linear", + "cl_avg_rv" : "linear", + "cl_avg_cloud_meanr" : "linear", + "cl_avg_cloud_stddevr" : "linear", + "cl_avg_act_meanr" : "linear", + "cl_avg_act_stddevr" : "linear", "cl_meanr" : "linear", "sd_conc" : "linear", } @@ -47,8 +53,6 @@ "prflux" : "linear", "cl_nc" : "linear", "cl_nr" : "linear", - "cl_nc_rico" : "linear", - "cl_nr_rico" : "linear", "clfrac" : "linear", "wvar" : "linear", "w3rd" : "linear", @@ -59,6 +63,8 @@ "er" : "linear", "wvarmax" : "linear", "surf_precip" : "linear", + "surf_flux_latent" : "linear", + "surf_flux_sensible" : "linear", "acc_precip" : "linear", "cloud_base" : "linear", "gccn_rw_cl" : "linear", @@ -71,13 +77,19 @@ "clb_bigrain_mean_conc" : "linear", "clb_bigrain_mean_inclt" : "linear", "clb_bigrain_mean_gccn_fraction" : "linear", - "cloud_cover_rico" : "linear", - "min_cloud_base_rico" : "linear", + "cloud_cover" : "linear", + "min_cloud_base" : "linear", "inversion_height_rico" : "linear", - "cl_acnv25_rico" : "linear", - "cl_accr25_rico" : "linear", + "cl_acnv25" : "linear", + "cl_accr25" : "linear", "RH_max" : "linear", - "cloud_avg_supersat" : "linear", + "cl_avg_supersat" : "linear", + "cl_avg_th" : "linear", + "cl_avg_rv" : "linear", + "cl_avg_cloud_meanr" : "linear", + "cl_avg_cloud_stddevr" : "linear", + "cl_avg_act_meanr" : "linear", + "cl_avg_act_stddevr" : "linear", "cl_meanr" : "linear", "sd_conc" : "linear", } @@ -89,8 +101,6 @@ "prflux" : None,#(0,20), "cl_nc" : None,#(0,90), "cl_nr" : None,#(0,90), - "cl_nc_rico" : None,#(0,90), - "cl_nr_rico" : None,#(0,90), "clfrac" : None, "wvar" : None, "w3rd" : None, @@ -101,7 +111,13 @@ "base_prflux_vs_clhght" : None,#(1,10000) "base_prflux_vs_clhght number of occurances" : None,#(1,10000) "RH_max" : None, - "cloud_avg_supersat" : None, + "cl_avg_supersat" : None, + "cl_avg_th" : None, + "cl_avg_rv" : None, + "cl_avg_cloud_meanr" : None, + "cl_avg_cloud_stddevr" : None, + "cl_avg_act_meanr" : None, + "cl_avg_act_stddevr" : None, "cl_meanr" : None, "sd_conc" : None, } @@ -113,8 +129,6 @@ "prflux" : (0,3000), "cl_nc" : (0,3000), "cl_nr" : (0,3000), - "cl_nc_rico" : (0,3000), - "cl_nr_rico" : (0,3000), "clfrac" : (0,3000), "wvar" : (0,3000), "w3rd" : (0,3000), @@ -125,52 +139,64 @@ "base_prflux_vs_clhght" : (0,2500), "base_prflux_vs_clhght number of occurances" : (0,2500), "RH_max" : (0,3000), - "cloud_avg_supersat" : (0,3000), + "cl_avg_supersat" : (0,3000), + "cl_avg_th" : (0,3000), + "cl_avg_rv" : (0,3000), + "cl_avg_cloud_meanr" : (0,3000), + "cl_avg_cloud_stddevr" : (0,3000), + "cl_avg_act_meanr" : (0,3000), + "cl_avg_act_stddevr" : (0,3000), "cl_meanr" : (0,3000), "sd_conc" : (0,3000), } xlimdict_series = { - "clfrac" : (0,5), - "cl_nc" : (0,5), - "cl_nr" : (0,5), - "cl_nc_rico" : (0,5), - "cl_nr_rico" : (0,5), - "lwp" : (0,5), - "rwp" : (0,5), - "er" : (0,5), - "wvarmax" : (0,5), - "surf_precip" : (0,5), - "acc_precip" : (0,5), - "cl_gccn_conc" : (0,5), - "cloud_base" : (0,5), - "clb_bigrain_mean_rd" : (0,5), - "clb_bigrain_mean_kappa" : (0,5), - "clb_bigrain_mean_conc" : (0,5), - "clb_bigrain_mean_inclt" : (0,5), - "clb_bigrain_mean_gccn_fraction" : (0,5), - "cloud_cover_rico" : (0,5), - "min_cloud_base_rico" : (0,5), - "inversion_height_rico" : (0,5), - "cl_acnv25_rico" : (0,5), - "cl_accr25_rico" : (0,5), - "RH_max" : (0,5), - "cloud_avg_supersat" : (0,5), - "cl_meanr" : (0,5), - "sd_conc" : (0,5), + "clfrac" : (0,24), + "cl_nc" : (0,24), + "cl_nr" : (0,24), + "lwp" : (0,24), + "rwp" : (0,24), + "er" : (0,24), + "wvarmax" : (0,24), + "surf_precip" : (0,24), + "surf_flux_latent" : (0,24), + "surf_flux_sensible" : (0,24), + "acc_precip" : (0,24), + "cl_gccn_conc" : (0,24), + "cloud_base" : (0,24), + "clb_bigrain_mean_rd" : (0,24), + "clb_bigrain_mean_kappa" : (0,24), + "clb_bigrain_mean_conc" : (0,24), + "clb_bigrain_mean_inclt" : (0,24), + "clb_bigrain_mean_gccn_fraction" : (0,24), + "cloud_cover" : (0,24), + "min_cloud_base" : (0,24), + "inversion_height_rico" : (0,24), + "cl_acnv25" : (0,24), + "cl_accr25" : (0,24), + "RH_max" : (0,24), + "cl_avg_supersat" : (0,24), + "cl_avg_th" : (0,24), + "cl_avg_rv" : (0,24), + "cl_avg_cloud_meanr" : (0,24), + "cl_avg_cloud_stddevr" : (0,24), + "cl_avg_act_meanr" : (0,24), + "cl_avg_act_stddevr" : (0,24), + "cl_meanr" : (0,24), + "sd_conc" : (0,24), } ylimdict_series = { "clfrac" : None, "cl_nc" : None,#(-1,155), "cl_nr" : None,#(-0.02,.55), - "cl_nc_rico" : None,#(-1,155), - "cl_nr_rico" : None,#(-0.02,.55), "lwp" : None,#(-1,45), "rwp" : None,#(-1,15), "er" : None, "wvarmax" : None, "surf_precip" : None,# (-0.05,1.2), + "surf_flux_latent" : (120,190), + "surf_flux_sensible" : None, "acc_precip" : None,#(-0.002,0.08), "cl_gccn_conc" : None,#(1e-6, 1), "cloud_base" : None, @@ -179,13 +205,19 @@ "clb_bigrain_mean_conc" : None, "clb_bigrain_mean_inclt" : None, "clb_bigrain_mean_gccn_fraction" : None, - "cloud_cover_rico" : None,#(-0.1,0.6), - "min_cloud_base_rico" : None,#(-1,750), + "cloud_cover" : None,#(-0.1,0.6), + "min_cloud_base" : None,#(-1,750), "inversion_height_rico" : None,#(1000,3000), - "cl_acnv25_rico" : None, - "cl_accr25_rico" : None, + "cl_acnv25" : None, + "cl_accr25" : None, "RH_max" : None, - "cloud_avg_supersat" : None, + "cl_avg_supersat" : None, + "cl_avg_th" : (296,302), + "cl_avg_rv" : (1.2e-2,1.6e-2), + "cl_avg_cloud_meanr" : None, + "cl_avg_cloud_stddevr" : None, + "cl_avg_act_meanr" : None, + "cl_avg_act_stddevr" : None, "cl_meanr" : None, "sd_conc" : None, } diff --git a/drawbicyc/drawbicyc.cpp b/drawbicyc/drawbicyc.cpp index 2bb53ba..a814446 100644 --- a/drawbicyc/drawbicyc.cpp +++ b/drawbicyc/drawbicyc.cpp @@ -42,7 +42,7 @@ int main(int argc, char** argv) // parse dir name std::string dir = vm["dir"].as(), - h5 = dir + "out_" + micro; + h5 = dir;// + "out_" + micro; // reading required plot types bool flag_series = vm["series"].as(), @@ -56,9 +56,24 @@ int main(int argc, char** argv) H5::DataSet h5d = h5f.openDataSet("G"); H5::DataSpace h5s = h5d.getSpace(); int NDims = h5s.getSimpleExtentNdims(); + + // selecting type of cloud mask + mask_type_t mask_type = mask_type_t::unset; + if(type == "dycoms") + { + std::cout << "Using Dycoms_rf02 cloud mask." << std::endl; + mask_type = mask_type_t::Dycoms_rf02; + } + else + { + std::cout << "Using Rico11 cloud mask." << std::endl; + mask_type = mask_type_t::Rico11; + } // detecting if subgrid model was on - bool sgs = true; +// std::cerr << "checking of sgs model was used... "; + bool sgs = h5f.nameExists("sgs");// true; + /* try { auto h5g = h5f.openGroup("sgs"); @@ -67,26 +82,28 @@ int main(int argc, char** argv) { sgs = false; } + */ + // std::cerr << sgs << std::endl; Plots plots(type, sgs); if(NDims == 2) { - 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)); + if(flag_series) plot_series(PlotterMask<2>(h5, micro, mask_type_t::Rico11), plots, type); + if(flag_profiles) plot_profiles(PlotterMask<2>(h5, micro, mask_type_t::Rico11), plots, type, normalize_prof); +// if(flag_fields) plot_fields(PlotterMask<2>(h5, micro), plots, type); +// if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMask<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_series) plot_series(PlotterMask<3>(h5, micro, mask_type), plots, type); + if(flag_profiles) plot_profiles(PlotterMask<3>(h5, micro, mask_type), plots, type, normalize_prof); +// if(flag_fields) plot_fields(PlotterMask<3>(h5, micro), plots, type); +// if(flag_qv_qc_2_6_10_min) plot_qv_qc_2_6_10_min(PlotterMask<2>(h5, micro)); // if(flag_lgrngn_spec) { -// plot_lgrngn_spec_positions(PlotterMicro_t<3>(h5, "lgrngn")); -// plot_lgrngn_spec(PlotterMicro_t<3>(h5, "lgrngn")); +// plot_lgrngn_spec_positions(PlotterMask<3>(h5, "lgrngn")); +// plot_lgrngn_spec(PlotterMask<3>(h5, "lgrngn")); // } } else diff --git a/drawbicyc/include/cases/Dycoms_RF02/plots.hpp b/drawbicyc/include/cases/Dycoms_RF02/plots.hpp index fae9ce4..33f30bb 100644 --- a/drawbicyc/include/cases/Dycoms_RF02/plots.hpp +++ b/drawbicyc/include/cases/Dycoms_RF02/plots.hpp @@ -4,12 +4,15 @@ const std::vector series_dycoms({ // "cl_nr", //"cl_acnv25_dycoms", //"cl_accr25_dycoms", -//"wvarmax", "cloud_cover_dycoms", +"wvarmax", "cloud_cover_dycoms", "lwp", - "er", + "er_dycoms", "surf_precip", // "acc_precip", "cl_nc" + , "surf_flux_latent" + , "surf_flux_sensible" + , "RH" // "cloud_base_dycoms", // "cloud_base_precip_dycoms" diff --git a/drawbicyc/include/cases/RICO11/plots.hpp b/drawbicyc/include/cases/RICO11/plots.hpp index e841b4f..df2aba2 100644 --- a/drawbicyc/include/cases/RICO11/plots.hpp +++ b/drawbicyc/include/cases/RICO11/plots.hpp @@ -1,30 +1,40 @@ #pragma once const std::vector series_rico({ -// "clb_bigrain_mean_rd", -// "clb_bigrain_mean_kappa", -// "clb_bigrain_mean_conc", -// "clb_bigrain_mean_inclt", -//, "clb_bigrain_mean_gccn_fraction" - -//"cl_acnv25_rico", -//"cl_accr25_rico", "RH_max", -"cloud_cover_rico", -"min_cloud_base_rico", +"cloud_cover", +"min_cloud_base", "inversion_height_rico", "tot_water", "lwp", "rwp", "surf_precip", -"acc_precip" -,"nc" -,"nr" -,"cl_nc_rico" -,"cl_nr_rico" -,"cloud_avg_supersat" -,"wvarmax" -,"cl_meanr" //TODO: zmienic maske na rico +"acc_precip", +"cl_nc", +"cl_nr", +"cl_avg_supersat", +"cl_avg_th", +"cl_avg_rv", +"cl_avg_cloud_meanr", +"cl_avg_cloud_stddevr", +"cl_avg_act_meanr", +"cl_avg_act_stddevr", +"wvarmax", + "surf_flux_latent", + "surf_flux_sensible", +"sd_conc" + + +// "clb_bigrain_mean_rd", +// "clb_bigrain_mean_kappa", +// "clb_bigrain_mean_conc", +// "clb_bigrain_mean_inclt", +//, "clb_bigrain_mean_gccn_fraction" + +//"cl_acnv25", +//"cl_accr25", +//,"nc" +//,"nr" //TODO (po usprawnieniu cloud mask i ujednoliceniu tego: /* ,"cl_avg_cloud_rad" @@ -34,12 +44,9 @@ const std::vector series_rico({ //,"rd_geq_0.8um_conc" //,"cl_rd_lt_0.8um_conc" //,"cl_rd_geq_0.8um_conc" -,"sd_conc" /* "cloud_base", - "surf_flux_latent", - "surf_flux_sensible" ,"cl_sd_conc" //"mass_dry", ,"cl_gccn_conc", "gccn_conc" diff --git a/drawbicyc/include/gnuplot_series_set_labels.hpp b/drawbicyc/include/gnuplot_series_set_labels.hpp index 49757f1..ace1ae0 100644 --- a/drawbicyc/include/gnuplot_series_set_labels.hpp +++ b/drawbicyc/include/gnuplot_series_set_labels.hpp @@ -14,7 +14,7 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt) gp << "set ylabel ''\n"; gp << "set title 'autoconersion rate (r>25um) (cloudy) [g/(m3*s}]'\n"; } - if (plt == "cloud_cover_dycoms" || plt == "cloud_cover_rico") + if (plt == "cloud_cover_dycoms" || plt == "cloud_cover") { // res_pos *= 60.; gp << "set xlabel ''\n"; @@ -63,13 +63,25 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt) gp << "set xlabel 'time [min]'\n"; gp << "set title 'relative std dev N_c'\n"; } - else if (plt == "cloud_avg_supersat") + else if (plt == "cl_avg_supersat") { gp << "set ylabel ' [%]'\n"; gp << "set xlabel 'time [min]'\n"; - gp << "set title 'avg supersaturation'\n"; + gp << "set title 'avg supersaturation in cloudy cells'\n"; } - else if (plt == "cloud_std_dev_supersat") + else if (plt == "cl_avg_th") + { + gp << "set ylabel 'theta [K]'\n"; + gp << "set xlabel 'time [min]'\n"; + gp << "set title 'avg theta in cloudy cells'\n"; + } + else if (plt == "cl_avg_rv") + { + gp << "set ylabel 'r_v [1]'\n"; + gp << "set xlabel 'time [min]'\n"; + gp << "set title 'avg rv in cloudy cells'\n"; + } + else if (plt == "cl_std_dev_supersat") { gp << "set ylabel 'sigma(S) / '\n"; gp << "set xlabel 'time [min]'\n"; @@ -87,6 +99,30 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt) gp << "set xlabel 'time [min]'\n"; gp << "set title 'average std dev of radius of activated droplets'\n"; } + else if (plt == "cl_avg_act_meanr") + { + gp << "set ylabel ' [um]'\n"; + gp << "set xlabel 'time [min]'\n"; + gp << "set title 'mean radius of activated droplets in cloud'\n"; + } + else if (plt == "cl_avg_act_stddevr") + { + gp << "set ylabel ' [um]'\n"; + gp << "set xlabel 'time [min]'\n"; + gp << "set title 'std. dev. of radius of activated droplets in cloud'\n"; + } + else if (plt == "cl_avg_cloud_meanr") + { + gp << "set ylabel ' [um]'\n"; + gp << "set xlabel 'time [min]'\n"; + gp << "set title 'mean radius of cloud droplets in cloud'\n"; + } + else if (plt == "cl_avg_cloud_stddevr") + { + gp << "set ylabel ' [um]'\n"; + gp << "set xlabel 'time [min]'\n"; + gp << "set title 'std. dev. of radius of cloud droplets in cloud'\n"; + } else if (plt == "sd_conc") { gp << "set ylabel ''\n"; @@ -253,7 +289,7 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt) gp << "set xlabel ''\n"; gp << "set ylabel ''\n"; } - else if (plt == "min_cloud_base_rico") + else if (plt == "min_cloud_base") { gp << "set title 'lowest cloud base [m]'\n"; gp << "set xlabel ''\n"; @@ -327,7 +363,7 @@ void gnuplot_series_set_labels(Gnuplot &gp, std::string plt) { gp << "set title 'latent surf flux [W/m^2]'\n"; } - else if (plt == "er") + else if (plt == "er_dycoms") { gp << "set title 'entrainment rate [cm / s]'\n"; gp << "set xlabel ''\n"; diff --git a/drawbicyc/include/plot_fields.hpp b/drawbicyc/include/plot_fields.hpp index 81f66c6..0c26edf 100644 --- a/drawbicyc/include/plot_fields.hpp +++ b/drawbicyc/include/plot_fields.hpp @@ -1,7 +1,7 @@ #include #include #include -#include +#include #include "plots.hpp" #include "gnuplot.hpp" @@ -36,7 +36,7 @@ void plot_fields(Plotter_t plotter, Plots plots, std::string type) { try{ // cloud water content - auto tmp = plotter.h5load_ract_timestep(at * n["outfreq"]) * 1e3; + auto tmp = plotter.load_ract_timestep(at * n["outfreq"]) * 1e3; std::string title = "cloud water mixing ratio [g/kg]"; gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(at) * n["outfreq"] * n["dt"] / 60.) << "min'\n"; diff --git a/drawbicyc/include/plot_lgrngn_spec.hpp b/drawbicyc/include/plot_lgrngn_spec.hpp index eea73e3..8eb9095 100644 --- a/drawbicyc/include/plot_lgrngn_spec.hpp +++ b/drawbicyc/include/plot_lgrngn_spec.hpp @@ -1,5 +1,5 @@ #include -#include +#include #include "plots.hpp" #include "focus.hpp" #include "gnuplot.hpp" @@ -26,7 +26,7 @@ void plot_lgrngn_spec_positions(Plotter_t plotter) try{ // cloud water content - auto tmp = plotter.h5load_ract_timestep(spectra_step * n["outfreq"]) * 1e3; + auto tmp = plotter.load_ract_timestep(spectra_step * n["outfreq"]) * 1e3; std::string title = "cloud water mixing ratio [g/kg]"; gp << "set title '" + title + " t = " << std::fixed << std::setprecision(2) << (double(spectra_step) * n["outfreq"] * n["dt"] / 60.) << "min'\n"; @@ -121,7 +121,7 @@ void plot_lgrngn_spec(Plotter_t plotter) // calc ratio of water content to adiabatic water content { - auto tmp = plotter.h5load_ract_timestep(spectra_step * n["outfreq"]) * 1e3; + auto tmp = plotter.load_ract_timestep(spectra_step * n["outfreq"]) * 1e3; double ratio = mean(tmp(focusBox)) / r_c_adiab; gp << "set label 4 'AF = " << std::fixed << std::setprecision(2) << ratio << "' at graph .2, .63 font \",15\"\n"; } diff --git a/drawbicyc/include/plot_prof.hpp b/drawbicyc/include/plot_prof.hpp index 2e7d086..3515190 100644 --- a/drawbicyc/include/plot_prof.hpp +++ b/drawbicyc/include/plot_prof.hpp @@ -1,4 +1,4 @@ -#include +#include #include //#include #include "plots.hpp" @@ -8,6 +8,7 @@ template void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool normalize) { + using arr_t = typename Plotter_t::arr_t; // read opts po::options_description opts("profile plotting options"); @@ -87,9 +88,9 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool if (plt == "rliq") { // liquid water content - // res += plotter.h5load_ra_timestep(at * n["outfreq"]) * 1e3; // aerosol - res += plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud - res += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; // rain + // res += plotter.load_ra_timestep(at * n["outfreq"]) * 1e3; // aerosol + res += plotter.load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud + res += plotter.load_rr_timestep(at * n["outfreq"]) * 1e3; // rain res_prof_hlpr = plotter.horizontal_mean(res); // average in x } if (plt == "gccn_conc") @@ -132,8 +133,7 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res_tmp = where(res_tmp > 0 , res_tmp / snap, res_tmp); } { - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - res_tmp2 = iscloudy_rc_rico(snap); + res_tmp2 = arr_t(plotter.get_mask(at * n["outfreq"])); res_tmp *= res_tmp2; } // mean only over downdraught cells @@ -149,8 +149,7 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res_tmp = where(res_tmp > 0 , res_tmp / snap, res_tmp); // mean radius } { - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - res_tmp2 = iscloudy_rc_rico(snap); + res_tmp2 = arr_t(plotter.get_mask(at * n["outfreq"])); res_tmp *= res_tmp2; } // mean only over downdraught cells @@ -206,8 +205,7 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res_tmp2 = isdowndraught(snap); } { // cloudy - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - res_tmp = iscloudy_rc_rico(snap); + res_tmp = arr_t(plotter.get_mask(at * n["outfreq"])); res_tmp2 *= res_tmp; } // mean rw @@ -235,8 +233,7 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res_tmp2 = isdowndraught(snap); } { // cloudy - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - res_tmp = iscloudy_rc_rico(snap); + res_tmp = arr_t(plotter.get_mask(at * n["outfreq"])); res_tmp2 *= res_tmp; } // mean rw @@ -264,8 +261,7 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res_tmp2 = isupdraught(snap); } { // cloudy - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - res_tmp = iscloudy_rc_rico(snap); + res_tmp = arr_t(plotter.get_mask(at * n["outfreq"])); res_tmp2 *= res_tmp; } // mean rw @@ -401,8 +397,7 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res_tmp = where(res_tmp > 0 , res_tmp / snap, res_tmp); } { - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - res_tmp2 = iscloudy_rc_rico(snap); + res_tmp2 = arr_t(plotter.get_mask(at * n["outfreq"])); res_tmp *= res_tmp2; } // mean only over downdraught cells @@ -423,8 +418,7 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool res_tmp = where(res_tmp > 0 , res_tmp / snap, res_tmp); } { - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - res_tmp2 = iscloudy_rc_rico(snap); + res_tmp2 = arr_t(plotter.get_mask(at * n["outfreq"])); res_tmp *= res_tmp2; } // mean only over downdraught cells @@ -443,10 +437,9 @@ 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 { - auto tmp = plotter.h5load_nc_timestep(at * n["outfreq"]); + auto tmp = plotter.load_nc_timestep(at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); res_tmp = snap; - res_tmp *= rhod / 1e6; // per cm^3 } // updraft only res_tmp *= res_tmp2; @@ -462,15 +455,12 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool } { - auto tmp = plotter.h5load_nc_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - res_tmp = snap; - snap = iscloudy(snap); // cloudiness mask + arr_t snap(plotter.get_mask(at * n["outfreq"])); res_tmp2 *= snap; // cloudy updrafts only } + res_tmp = typename Plotter_t::arr_t(plotter.load_nc_timestep(at * n["outfreq"])); + // mean only over cloudy updraught cells prof_tmp = plotter.horizontal_sum(res_tmp2); // number of downdraft cells on a given level @@ -490,10 +480,9 @@ 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 { - auto tmp = plotter.h5load_nc_timestep(at * n["outfreq"]); + auto tmp = plotter.load_nc_timestep(at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); res_tmp = snap; - res_tmp *= rhod / 1e6; // per cm^3 } // updraft only res_tmp *= res_tmp2; @@ -605,19 +594,19 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool } else if (plt == "sat_RH") { - res = plotter.h5load_RH_timestep(at * n["outfreq"]); + res = plotter.load_RH_timestep(at * n["outfreq"]); res = (res -1) * 100; res_prof_hlpr = plotter.horizontal_mean(res); // average in x } else if (plt == "RH") { - res = plotter.h5load_RH_timestep(at * n["outfreq"]); + res = plotter.load_RH_timestep(at * n["outfreq"]); res_prof_hlpr = plotter.horizontal_mean(res) * 100; // average in x; [%] } else if (plt == "sat_RH_up") { { - auto tmp = plotter.h5load_RH_timestep(at * n["outfreq"]); + auto tmp = plotter.load_RH_timestep(at * n["outfreq"]); typename Plotter_t::arr_t snap(tmp); res_tmp = (snap -1) * 100; } @@ -637,9 +626,9 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool } else if (plt == "00rtot") { - // res = plotter.h5load_ra_timestep(at * n["outfreq"]) * 1e3; // aerosol - res = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud - res += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; // rain + // res = plotter.load_ra_timestep(at * n["outfreq"]) * 1e3; // aerosol + res = plotter.load_rc_timestep(at * n["outfreq"]) * 1e3; // cloud + res += plotter.load_rr_timestep(at * n["outfreq"]) * 1e3; // rain res += plotter.h5load_timestep("rv", at * n["outfreq"]) * 1e3; // vapour res_prof_hlpr = plotter.horizontal_mean(res); // average in x @@ -659,7 +648,7 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool else if (plt == "N_c") { // cloud drops concentration [1/cm^3] - res = plotter.h5load_nc_timestep(at * n["outfreq"]) * rhod / 1e6; // from sepcific to normal moment + per cm^3 + res = plotter.load_nc_timestep(at * n["outfreq"]); // from sepcific to normal moment + per cm^3 res_prof_hlpr = plotter.horizontal_mean(res); // average in x } else if (plt == "rd_geq_0.8um_conc") @@ -679,15 +668,8 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool // cloud droplet (0.5um < r < 25 um) concentration in cloudy grid cells try { - // cloud fraction (cloudy if N_c > 20/cm^3) - auto tmp = plotter.h5load_nc_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - typename Plotter_t::arr_t snap2; - snap2.resize(snap.shape()); - snap2=snap; - snap = iscloudy(snap); // cloudiness mask + typename Plotter_t::arr_t snap2(plotter.load_nc_timestep(at * n["outfreq"])); + arr_t snap(plotter.get_mask(at * n["outfreq"])); snap2 *= snap; // mean only over cloudy cells @@ -701,9 +683,9 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool // liquid potential temp [K] { auto &ql(res_tmp2); - ql = plotter.h5load_rc_timestep(at * n["outfreq"]); // cloud -// ql += plotter.h5load_ra_timestep(at * n["outfreq"]); // aerosol - ql += plotter.h5load_rr_timestep(at * n["outfreq"]); // rain + ql = plotter.load_rc_timestep(at * n["outfreq"]); // cloud +// ql += plotter.load_ra_timestep(at * n["outfreq"]); // aerosol + ql += plotter.load_rr_timestep(at * n["outfreq"]); // rain // ql is now q_l (liq water content) // auto tmp = plotter.h5load_timestep("th", at * n["outfreq"]); // typename Plotter_t::arr_t th_d(tmp); @@ -737,11 +719,7 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool { // cloud fraction (cloudy if N_c > 20/cm^3) { - auto tmp = plotter.h5load_nc_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - snap = iscloudy(snap); + arr_t snap(plotter.get_mask(at * n["outfreq"])); res += snap; } res_prof_hlpr = plotter.horizontal_mean(res); // average in x @@ -750,15 +728,13 @@ void plot_profiles(Plotter_t plotter, Plots plots, std::string type, const bool // TODO: 'normalize' messes with this plot { res_prof_hlpr = 0; - // cloudy cells (cloudy if q_c > 0.01 g/kg as in RICO paper. NOTE: also add q_r ?) - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - snap = iscloudy_rc_rico(snap); // cloudiness mask + arr_t snap(plotter.get_mask(at * n["outfreq"])); plotter.k_i = blitz::sum(snap, plotter.LastIndex); // sum in the vertical, assumes that all cloudy cells in a column belong to the same cloud plotter.tmp_int_hrzntl_slice = blitz::first(snap > 0, plotter.LastIndex); // cloud base hgt over dz // precipitation flux(doesnt include vertical velocity w!) - res = plotter.h5load_prflux_timestep(at * n["outfreq"]); + res = plotter.load_prflux_timestep(at * n["outfreq"]); plotter.tmp_float_hrzntl_slice = plotter.get_value_at_hgt(res, plotter.tmp_int_hrzntl_slice); // precip flux at cloud base // NOTE: we assume that k_i and tmp_float_hr... is contiguous in memory @@ -776,7 +752,7 @@ 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.load_prflux_timestep(at * n["outfreq"]); res_prof_hlpr = plotter.horizontal_mean(res); // average in x } // add vertical velocity to precipitation flux (3rd mom of cloud drops * w) diff --git a/drawbicyc/include/plot_qv_qc_2_6_10_min.hpp b/drawbicyc/include/plot_qv_qc_2_6_10_min.hpp index 8542acf..8bfb01a 100644 --- a/drawbicyc/include/plot_qv_qc_2_6_10_min.hpp +++ b/drawbicyc/include/plot_qv_qc_2_6_10_min.hpp @@ -1,7 +1,7 @@ #include #include -#include +#include #include #include "plots.hpp" #include "gnuplot.hpp" @@ -103,7 +103,7 @@ void plot_qv_qc_2_6_10_min(Plotter_t plotter) gp << "set title offset 0, -0.8 '$q_c$ [g/kg], t = 2 min'\n"; try{ // cloud water content - typename Plotter_t::arr_t tmp(plotter.h5load_ract_timestep(60) * 1e3); + typename Plotter_t::arr_t tmp(plotter.load_ract_timestep(60) * 1e3); std::cout << tmp; plotter.plot(gp, tmp); } @@ -118,7 +118,7 @@ void plot_qv_qc_2_6_10_min(Plotter_t plotter) gp << "set title offset 0, -0.8 '$q_c$ [g/kg], t = 6 min'\n"; try{ // cloud water content - typename Plotter_t::arr_t tmp(plotter.h5load_ract_timestep(360) * 1e3); + typename Plotter_t::arr_t tmp(plotter.load_ract_timestep(360) * 1e3); std::cout << tmp; plotter.plot(gp, tmp); } @@ -133,7 +133,7 @@ void plot_qv_qc_2_6_10_min(Plotter_t plotter) gp << "set title offset 0, -0.8 '$q_c$ [g/kg], t = 10 min'\n"; try{ // cloud water content - typename Plotter_t::arr_t tmp (plotter.h5load_ract_timestep(600) * 1e3); + typename Plotter_t::arr_t tmp (plotter.load_ract_timestep(600) * 1e3); std::cout << tmp; plotter.plot(gp, tmp); } diff --git a/drawbicyc/include/plot_series.hpp b/drawbicyc/include/plot_series.hpp index 2167e77..f01e27c 100644 --- a/drawbicyc/include/plot_series.hpp +++ b/drawbicyc/include/plot_series.hpp @@ -1,4 +1,4 @@ -#include +#include #include #include "plots.hpp" #include "gnuplot_series_set_labels.hpp" @@ -7,6 +7,8 @@ template void plot_series(Plotter_t plotter, Plots plots, std::string type) { + using arr_t = typename Plotter_t::arr_t; + auto& n = plotter.map; auto& n_prof = plotter.map_prof; for(auto elem : n) @@ -24,15 +26,15 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // read in density auto tmp = plotter.h5load(plotter.file + "/const.h5", "G"); - typename Plotter_t::arr_t rhod(tmp); - typename Plotter_t::arr_t rtot(rhod.shape()); + arr_t rhod(tmp); + arr_t rtot(rhod.shape()); - typename Plotter_t::arr_t res_tmp(rhod.shape()); + 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; + std::vector prev_u_vec, prev_w_vec; // container for the running sum - typename Plotter_t::arr_t run_sum_u(rhod.shape()), run_sum_w(rhod.shape()); + arr_t run_sum_u(rhod.shape()), run_sum_w(rhod.shape()); run_sum_u = 0; run_sum_w = 0; // number of timesteps over which the running avg of u and w is calculated (1 min interval) @@ -122,29 +124,27 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) for (auto &plt : plots.series) { + std::cerr << plt << std::endl; if (plt == "cloud_cover_dycoms") { try { - auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; //g/kg - typename Plotter_t::arr_t snap(tmp); - snap += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; //g/kg - snap *= rhod; // water per cubic metre (should be wet density...) - plotter.k_i = blitz::sum(snap, plotter.LastIndex) * n["dz"]; // LWP [g/m2] in the column + auto tmp = plotter.load_rc_timestep(at * n["outfreq"]) * 1e3; //g/kg + arr_t snap(tmp); + snap += plotter.load_rr_timestep(at * n["outfreq"]) * 1e3; //g/kg + plotter.multiply_by_rhod(snap); + plotter.k_i = blitz::sum(snap, plotter.LastIndex) * n["refined dz"]; // LWP [g/m2] in the column plotter.k_i = where(plotter.k_i > 20 , 1 , 0); // cloudiness as in Ackermann et al. res_series[plt](at) = blitz::mean(plotter.k_i); } catch(...) {if(at==first_timestep) data_found[plt]=0;} } - else if (plt == "cloud_cover_rico") + else if (plt == "cloud_cover") { try { // cloud fraction (fraction of columns with at least one cloudy cell, i.e. cell with q_c > 0.01 g/kg) - auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap = iscloudy_rc_rico(snap); // find cells with rc>1e-5 - snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux + arr_t snap(plotter.get_mask(at * n["outfreq"])); plotter.k_i = blitz::first((snap == 1), plotter.LastIndex); plotter.k_i = blitz::sum(snap, plotter.LastIndex); plotter.k_i = where(plotter.k_i > 0, 1, 0); @@ -160,13 +160,13 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // read RH auto tmp = plotter.h5load_timestep("RH", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); + arr_t snap(tmp); res_series[plt](at) = blitz::max(snap); } catch(...) {if(at==first_timestep) data_found[plt]=0;} } - // r_act averaged over cloudy cells - else if (plt == "ract_avg") + // mixing ratio of acivated droplets averaged over cloudy cells + else if (plt == "cl_ract") { try { @@ -177,14 +177,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) catch(...) {if(at==first_timestep) data_found[plt]=0;} } // cloud top height - else if (plt =="cloud_top_height") + else if (plt =="cl_top_height") { try { - // cloud fraction (cloudy if ql > 1e-5)) - auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap = iscloudy_rc_rico(snap); + arr_t snap(plotter.get_mask(at * n["outfreq"])); plotter.k_i = blitz::last((snap == 1), plotter.LastIndex); auto cloudy_column = plotter.k_i.copy(); cloudy_column = blitz::sum(snap, plotter.LastIndex); @@ -237,27 +234,27 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) /* { auto tmp = plotter.h5load_timestep("aerosol_rw_mom3", at * n["outfreq"]) * 4./3. * 3.1416 * 1e3; - typename Plotter_t::arr_t snap(tmp); + arr_t snap(tmp); snap *= rhod; res_series[plt](at) = blitz::mean(snap); } */ { auto tmp = plotter.h5load_timestep("cloud_rw_mom3", at * n["outfreq"]) * 4./3. * 3.1416 * 1e3; - typename Plotter_t::arr_t snap(tmp); - snap *= rhod; + arr_t snap(tmp); + plotter.multiply_by_rhod(snap); res_series[plt](at) += blitz::mean(snap); } { auto tmp = plotter.h5load_timestep("rain_rw_mom3", at * n["outfreq"]) * 4./3. * 3.1416 * 1e3; - typename Plotter_t::arr_t snap(tmp); - snap *= rhod; + arr_t snap(tmp); + plotter.multiply_by_rhod(snap); res_series[plt](at) += blitz::mean(snap); } { auto tmp = plotter.h5load_timestep("rv", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap *= rhod; + arr_t snap(tmp); + plotter.multiply_by_rhod(snap); res_series[plt](at) += blitz::mean(snap); } } @@ -277,10 +274,10 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // vertical velocity at the center of mass of activated droplets try { - auto tmp = plotter.h5load_ract_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - typename Plotter_t::arr_t snap2(tmp); - typename Plotter_t::arr_t snap3(tmp); + auto tmp = plotter.load_ract_timestep(at * n["outfreq"]); + arr_t snap(tmp); + arr_t snap2(tmp); + arr_t snap3(tmp); snap2 = snap2 * plotter.LastIndex; snap3 = snap3 * blitz::tensor::i; @@ -289,7 +286,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) int z_idx = blitz::sum(snap2) / blitz::sum(snap); int x_idx = blitz::sum(snap3) / blitz::sum(snap); auto tmp2 = plotter.h5load_timestep("w", at * n["outfreq"]); - typename Plotter_t::arr_t snap_mom(tmp2); + arr_t snap_mom(tmp2); res_series[plt](at) = snap_mom(x_idx, z_idx); } else @@ -302,10 +299,10 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // supersaturation at the center of mass of activated droplets try { - auto tmp = plotter.h5load_ract_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - typename Plotter_t::arr_t snap2(tmp); - typename Plotter_t::arr_t snap3(tmp); + auto tmp = plotter.load_ract_timestep(at * n["outfreq"]); + arr_t snap(tmp); + arr_t snap2(tmp); + arr_t snap3(tmp); snap2 = snap2 * plotter.LastIndex; snap3 = snap3 * blitz::tensor::i; @@ -314,7 +311,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) int z_idx = blitz::sum(snap2) / blitz::sum(snap); int x_idx = blitz::sum(snap3) / blitz::sum(snap); auto tmp2 = plotter.h5load_timestep("RH", at * n["outfreq"]); - typename Plotter_t::arr_t snap_mom(tmp2); + arr_t snap_mom(tmp2); res_series[plt](at) = snap_mom(x_idx, z_idx) - 1; } else @@ -327,10 +324,10 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // 0th moment of rw distribution at the center of mass of activated droplets (particles concentration), 2D only try { - auto tmp = plotter.h5load_ract_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - typename Plotter_t::arr_t snap2(tmp); - typename Plotter_t::arr_t snap3(tmp); + auto tmp = plotter.load_ract_timestep(at * n["outfreq"]); + arr_t snap(tmp); + arr_t snap2(tmp); + arr_t snap3(tmp); snap2 = snap2 * plotter.LastIndex; snap3 = snap3 * blitz::tensor::i; @@ -340,9 +337,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) com_x_idx(at) = blitz::sum(snap3) / blitz::sum(snap); std::cout << at << ": (" << com_x_idx(at) << "," << com_z_idx(at) << ")" << std::endl; auto tmp2 = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]); - typename Plotter_t::arr_t snap_mom(tmp2); + arr_t snap_mom(tmp2); com_N_c(at) = snap_mom(com_x_idx(at), com_z_idx(at)); // 0th raw moment / mass [1/kg] - snap_mom *= rhod; // now per m^3 + plotter.multiply_by_rhod(snap); res_series[plt](at) = snap_mom(com_x_idx(at), com_z_idx(at)); } else @@ -359,7 +356,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { auto tmp = plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); // 1st raw moment / mass [m / kg] + arr_t snap(tmp); // 1st raw moment / mass [m / kg] if(com_N_c(at) > 0) res_series[plt](at) = snap(com_x_idx(at), com_z_idx(at)) / com_N_c(at); else @@ -374,13 +371,13 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at * n["outfreq"]); - typename Plotter_t::arr_t zeroth_raw_mom(tmp); // 0th raw moment / mass [1 / kg] + arr_t zeroth_raw_mom(tmp); // 0th raw moment / mass [1 / kg] tmp = plotter.h5load_timestep("actrw_rw_mom1", at * n["outfreq"]); - typename Plotter_t::arr_t first_raw_mom(tmp); // 1st raw moment / mass [m / kg] + arr_t first_raw_mom(tmp); // 1st raw moment / mass [m / kg] tmp = plotter.h5load_timestep("actrw_rw_mom2", at * n["outfreq"]); - typename Plotter_t::arr_t second_raw_mom(tmp); // 2nd raw moment / mass [m^2 / kg] + arr_t second_raw_mom(tmp); // 2nd raw moment / mass [m^2 / kg] tmp = plotter.h5load_timestep("sd_conc", at * n["outfreq"]); - typename Plotter_t::arr_t sd_conc(tmp); // number of SDs + arr_t sd_conc(tmp); // number of SDs if(com_N_c(at) > 0) { double SD_no = sd_conc(com_x_idx(at), com_z_idx(at)); @@ -413,7 +410,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { tmp = plotter.h5load_timestep("sd_conc", at * n["outfreq"]); - typename Plotter_t::arr_t sd_conc(tmp); // number of SDs + arr_t sd_conc(tmp); // number of SDs res_series[plt](at) = sd_conc(com_x_idx(at), com_z_idx(at)); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -424,7 +421,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { auto tmp = plotter.h5load_timestep("th", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); + arr_t snap(tmp); res_tmp = is_th_prtrb(snap); // find cells with th>300.1 snap *= res_tmp; // apply filter @@ -442,7 +439,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // height of cells with largest gradient of theta try { - typename Plotter_t::arr_t th(plotter.h5load_timestep("th", at * n["outfreq"])); + arr_t th(plotter.h5load_timestep("th", at * n["outfreq"])); auto grad = plotter.cent_diff_vert(th); auto max_index = blitz::maxIndex(grad, plotter.LastIndex); res_series[plt](at) = (blitz::mean(max_index) + 1) * n["dz"]; @@ -455,9 +452,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); + arr_t snap(tmp); snap /= 1e6; // per cm^3 - snap *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap); res_series[plt](at) = blitz::mean(snap); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -468,9 +465,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { auto tmp = plotter.h5load_timestep("all_rw_mom0", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); + arr_t snap(tmp); snap /= 1e6; // per cm^3 - snap *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap); res_series[plt](at) = blitz::mean(snap); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -481,9 +478,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { auto tmp = plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); + arr_t snap(tmp); snap /= 1e6; // per cm^3 - snap *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap); res_series[plt](at) = blitz::mean(snap); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -493,43 +490,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // cloud droplet (0.5um < r < 25 um) concentration in cloudy grid cells try { - // cloud fraction (cloudy if N_c > 20/cm^3) - auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - typename Plotter_t::arr_t snap2; - snap2.resize(snap.shape()); - snap2=snap; - snap = iscloudy(snap); // cloudiness mask - snap2 *= snap; - if(blitz::sum(snap) > 0) - res_series[plt](at) = blitz::sum(snap2) / blitz::sum(snap); - else - res_series[plt](at) = 0; - } - catch(...) {if(at==first_timestep) data_found[plt]=0;} - } - else if (plt == "cl_nc_rico") - { - // cloud droplet (0.5um < r < 25 um) concentration in cloudy grid cells - try - { - // rico cloud mask - auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap = iscloudy_rc_rico(snap); // find cells with rc>1e-5 - snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux - - typename Plotter_t::arr_t snap2(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - snap2 *= rhod; // b4 it was specific moment - snap2 /= 1e6; // per cm^3 - - snap2 *= snap; - if(blitz::sum(snap) > 0) - res_series[plt](at) = blitz::sum(snap2) / blitz::sum(snap); - else - res_series[plt](at) = 0; + auto stats = plotter.cloud_cloudconc_stats_timestep(at * n["outfreq"]); + res_series[plt](at) = stats.first; + res_series_std_dev[plt](at) = stats.second; } catch(...) {if(at==first_timestep) data_found[plt]=0;} } @@ -538,45 +501,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // rain drop (25um < r) concentration in cloudy grid cells try { - // cloud fraction (cloudy if N_c > 20/cm^3) - auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - snap = iscloudy(snap); // cloudiness mask - auto tmp2 = plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"]); - typename Plotter_t::arr_t snap2(tmp2); - snap2 *= rhod; // b4 it was specific moment - snap2 /= 1e6; // per cm^3 - snap2 *= snap; - if(blitz::sum(snap) > 0) - res_series[plt](at) = blitz::sum(snap2) / blitz::sum(snap); - else - res_series[plt](at) = 0; - } - catch(...) {if(at==first_timestep) data_found[plt]=0;} - } - else if (plt == "cl_nr_rico") - { - // rain drop (25um < r) concentration in cloudy grid cells - try - { - // rico cloud mask - auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap = iscloudy_rc_rico(snap); // find cells with rc>1e-5 - snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux - - auto tmp2 = plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"]); - typename Plotter_t::arr_t snap2(tmp2); - snap2 *= rhod; // b4 it was specific moment - snap2 /= 1e6; // per cm^3 - - snap2 *= snap; - if(blitz::sum(snap) > 0) - res_series[plt](at) = blitz::sum(snap2) / blitz::sum(snap); - else - res_series[plt](at) = 0; + auto stats = plotter.cloud_rainconc_stats_timestep(at * n["outfreq"]); + res_series[plt](at) = stats.first; + res_series_std_dev[plt](at) = stats.second; } catch(...) {if(at==first_timestep) data_found[plt]=0;} } @@ -586,48 +513,36 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { auto tmp = plotter.h5load_timestep("actrw_rw_mom0", at *n["outfreq"]) * rhod; - typename Plotter_t::arr_t snap(tmp); + arr_t snap(tmp); res_series[plt](at) = blitz::sum(snap)*plotter.CellVol; } catch(...) {if(at==first_timestep) data_found[plt]=0;} } - else if (plt == "cloud_base_dycoms") + else if (plt == "cloud_base") { // average cloud base in the domain try { - // cloud fraction (cloudy if N_c > 20/cm^3) - auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - snap = iscloudy(snap); // cloudiness mask - snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux + arr_t snap(plotter.get_mask(at * n["outfreq"])); plotter.k_i = blitz::first((snap == 1), plotter.LastIndex); auto cloudy_column = plotter.k_i.copy(); cloudy_column = blitz::sum(snap, plotter.LastIndex); cloudy_column = where(cloudy_column > 0, 1, 0); plotter.k_i = where(cloudy_column == 0, 0, plotter.k_i); if(blitz::sum(cloudy_column) > 0) - res_series[plt](at) = double(blitz::sum(plotter.k_i)) / blitz::sum(cloudy_column) * n["dz"]; + res_series[plt](at) = double(blitz::sum(plotter.k_i)) / blitz::sum(cloudy_column) * n["refined dz"]; else res_series[plt](at) = 0; } catch(...) {if(at==first_timestep) data_found[plt]=0;} } - else if (plt == "cloud_base_precip_dycoms") + else if (plt == "cloud_base_precip") { // domain-averaged (all columns) precipitation at the average cloud base height [mm/day] try { - // -- average cloud base, almost exactly as in "cloud_base_dycoms"... -- - // cloud fraction (cloudy if N_c > 20/cm^3) - auto tmp = plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - snap = iscloudy(snap); // cloudiness mask - snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux + // -- average cloud base, almost exactly as in "cloud_base"... -- + arr_t snap(plotter.get_mask(at * n["outfreq"])); plotter.k_i = blitz::first((snap == 1), plotter.LastIndex); auto cloudy_column = plotter.k_i.copy(); cloudy_column = blitz::sum(snap, plotter.LastIndex); @@ -644,35 +559,32 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) else { // -- precipitation at this height averaged over all cells, cloudy or not -- - auto prflux = plotter.h5load_prflux_timestep(at * n["outfreq"]); // prflux in [W/m^2] + auto prflux = plotter.load_prflux_timestep(at * n["outfreq"]); // prflux in [W/m^2] res_series[plt](at) = blitz::mean(prflux(plotter.hrzntl_slice(cloud_base_idx))) / 2264.705 * 3.6 * 24; // convert to [mm/day] } } catch(...) {if(at==first_timestep) data_found[plt]=0;} } - else if (plt == "min_cloud_base_rico") + else if (plt == "min_cloud_base") { // lowest cloud base in the domain try { - // cloud fraction (cloudy if r_c > 1e-5) - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - snap = iscloudy_rc_rico(snap); - snap(plotter.hrzntl_slice(0)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux + arr_t snap(plotter.get_mask(at * n["outfreq"])); plotter.k_i = blitz::first((snap == 1), plotter.LastIndex); auto cloudy_column = plotter.k_i.copy(); cloudy_column = blitz::sum(snap, plotter.LastIndex); cloudy_column = where(cloudy_column > 0, 1, 0); plotter.k_i = where(cloudy_column == 0, 1e6, plotter.k_i); // 1e6 denotes no clouds in the column if(blitz::sum(cloudy_column) > 0) - res_series[plt](at) = blitz::min(plotter.k_i) * n["dz"]; + res_series[plt](at) = blitz::min(plotter.k_i) * n["refined dz"]; else res_series[plt](at) = 0; } catch(...) {if(at==first_timestep) data_found[plt]=0;} } // average concentration of activated droplets in cloudy cells - else if (plt == "cloud_avg_act_conc") + else if (plt == "cl_avg_act_conc") { try { @@ -682,34 +594,77 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) } catch(...) {if(at==first_timestep) data_found[plt]=0;} } - // average supersaturation in cells with S>0 - else if (plt == "cloud_avg_supersat") + // average supersaturation in cloudy cells + else if (plt == "cl_avg_supersat") + { + try + { + auto stats = plotter.cloud_supersat_stats_timestep(at * n["outfreq"]); + res_series[plt](at) = stats.first; + res_series_std_dev[plt](at) = stats.second; + } + catch(...) {if(at==first_timestep) data_found[plt]=0;} + } + // average theta in cloudy cells + else if (plt == "cl_avg_th") + { + try + { + auto stats = plotter.cloud_theta_stats_timestep(at * n["outfreq"]); + res_series[plt](at) = stats.first; + res_series_std_dev[plt](at) = stats.second; + } + catch(...) {if(at==first_timestep) data_found[plt]=0;} + } + // average rv in cloudy cells + else if (plt == "cl_avg_rv") { try { - auto stats = plotter.positive_supersat_stats_timestep(at * n["outfreq"]); + auto stats = plotter.cloud_rv_stats_timestep(at * n["outfreq"]); res_series[plt](at) = stats.first; res_series_std_dev[plt](at) = stats.second; } catch(...) {if(at==first_timestep) data_found[plt]=0;} } // spatial average of mean radius of activated droplets in cloudy cells - else if (plt == "cl_avg_cloud_rad") + else if (plt == "cl_avg_act_meanr") + { + try + { + auto stats = plotter.cloud_actmeanr_stats_timestep(at * n["outfreq"]); + res_series[plt](at) = stats.first; + } + catch(...) {if(at==first_timestep) data_found[plt]=0;} + } + + // spatial average of standard deviation of activated droplet radius distribution in cloudy cells + else if (plt == "cl_avg_act_stddevr") { try { - auto stats = plotter.cloud_meanr_stats_timestep(at * n["outfreq"]); + auto stats = plotter.cloud_actstddevr_stats_timestep(at * n["outfreq"]); + res_series[plt](at) = stats.first; + } + catch(...) {if(at==first_timestep) data_found[plt]=0;} + } + // spatial average of mean radius of cloud droplets in cloudy cells + else if (plt == "cl_avg_cloud_meanr") + { + try + { + auto stats = plotter.cloud_cloudmeanr_stats_timestep(at * n["outfreq"]); res_series[plt](at) = stats.first; } catch(...) {if(at==first_timestep) data_found[plt]=0;} } - // spatial average of standard deviation of the acttivated droplets radius distribution in cloudy cells - else if (plt == "cloud_avg_std_dev_act_rad") + // spatial average of standard deviation of cloud droplet radius distribution in cloudy cells + else if (plt == "cl_avg_cloud_stddevr") { try { - auto stats = plotter.cloud_stddevr_stats_timestep(at * n["outfreq"]); + auto stats = plotter.cloud_cloudstddevr_stats_timestep(at * n["outfreq"]); res_series[plt](at) = stats.first; } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -722,8 +677,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { auto tmp = plotter.h5load_timestep("rd_rng000_mom3", at * n["outfreq"]) * 4./3. * 3.14 * rho_dry * 1e3; - typename Plotter_t::arr_t snap(tmp); - snap *= rhod * plotter.CellVol; // turn mixing ratio in g/kg to total mass in g + arr_t snap(tmp); + plotter.multiply_by_rhod(snap); + plotter.multiply_by_CellVol(snap); res_series[plt](at) = blitz::sum(snap); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -755,71 +711,17 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) } catch(...) {if(at==first_timestep) data_found[plt]=0;} } - else if (plt == "cl_acnv25_dycoms") - { - // autconversion rate with rain threshold r=25um (cloudy cells as in Dycoms) [g/(m3*s)] - // rather coarse estimate, sum of acnv accumulated over ALL cells since the last output - // is divided by the instantaneous volume of all cloudy cells - // TODO: output instantaneous acnv rate in libcloud, not the accumulated one? - 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"])); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - snap = iscloudy(snap); // cloudiness mask - - typename Plotter_t::arr_t acc_acnv(plotter.h5load_timestep("acc_acnv25", at * n["outfreq"])); - auto tot_acc_acnv = blitz::sum(acc_acnv); - - if(blitz::sum(snap) > 0) - res_series[plt](at) = 4./3. * 3.1416 * 1e6 * (tot_acc_acnv - tot_acc_acnv_prev) / ((blitz::sum(snap) * plotter.CellVol) * (n["outfreq"] * n["dt"])); - else - res_series[plt](at) = 0; - - tot_acc_acnv_prev = tot_acc_acnv; - } - catch(...) {if(at==first_timestep) data_found[plt]=0;} - } - else if (plt == "cl_accr25_dycoms") + else if (plt == "cl_acnv25") { - // accretion rate with rain threshold r=25um (cloudy cells as in Dycoms) [g/(m3*s)] - // rather coarse estimate, sum of accr accumulated over ALL cells since the last output - // is divided by the instantaneous volume of all cloudy cells - // TODO: output instantaneous accr rate in libcloud, not the accumulated one? - 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"])); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - snap = iscloudy(snap); // cloudiness mask - - typename Plotter_t::arr_t acc_accr(plotter.h5load_timestep("acc_accr25", at * n["outfreq"])); - double tot_acc_accr = blitz::sum(acc_accr); - - if(blitz::sum(snap) > 0) - res_series[plt](at) = 4./3. * 3.14166 * 1e6 * (tot_acc_accr - tot_acc_accr_prev) / ((blitz::sum(snap) * plotter.CellVol) * (n["outfreq"] * n["dt"])); - else - res_series[plt](at) = 0; - - tot_acc_accr_prev = tot_acc_accr; - } - catch(...) {if(at==first_timestep) data_found[plt]=0;} - } - else if (plt == "cl_acnv25_rico") - { - // autconversion rate with rain threshold r=25um (cloudy cells as in Rico) [g/(m3*s)] + // autconversion rate with rain threshold r=25um [g/(m3*s)] // rather coarse estimate, sum of acnv accumulated over ALL cells since the last output // is divided by the instantaneous volume of all cloudy cells // TODO: output instantaneous acnv rate in libcloud, not the accumulated one? try { - // cloud fraction (cloudy if r_c > 1e-5) - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - snap = iscloudy_rc_rico(snap); + arr_t snap(plotter.get_mask(at * n["outfreq"])); // cloud mask - typename Plotter_t::arr_t acc_acnv(plotter.h5load_timestep("acc_acnv25", at * n["outfreq"])); + arr_t acc_acnv(plotter.h5load_timestep("acc_acnv25", at * n["outfreq"])); auto tot_acc_acnv = blitz::sum(acc_acnv); if(blitz::sum(snap) > 0) @@ -831,19 +733,17 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) } catch(...) {if(at==first_timestep) data_found[plt]=0;} } - else if (plt == "cl_accr25_rico") + else if (plt == "cl_accr25") { - // accretion rate with rain threshold r=25um (cloudy cells as in Rico) [g/(m3*s)] + // accretion rate with rain threshold r=25um [g/(m3*s)] // rather coarse estimate, sum of accr accumulated over ALL cells since the last output // is divided by the instantaneous volume of all cloudy cells // TODO: output instantaneous accr rate in libcloud, not the accumulated one? try { - // cloud fraction (cloudy if r_c > 1e-5) - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - snap = iscloudy_rc_rico(snap); + arr_t snap(plotter.get_mask(at * n["outfreq"])); // cloud mask - typename Plotter_t::arr_t acc_accr(plotter.h5load_timestep("acc_accr25", at * n["outfreq"])); + arr_t acc_accr(plotter.h5load_timestep("acc_accr25", at * n["outfreq"])); double tot_acc_accr = blitz::sum(acc_accr); if(blitz::sum(snap) > 0) @@ -861,11 +761,11 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { { - auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; //g/kg - typename Plotter_t::arr_t snap(tmp); - snap += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; //g/kg - snap *= rhod; // water per cubic metre (should be wet density...) - res_series[plt](at) = blitz::mean(snap); + arr_t rl(plotter.load_rc_timestep(at * n["outfreq"])); + arr_t rr(plotter.load_rr_timestep(at * n["outfreq"])); + rl = (rl + rr) * 1e3; // g/kg + plotter.multiply_by_rhod(rl); + res_series[plt](at) = blitz::mean(rl); } } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -876,8 +776,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { { - typename Plotter_t::arr_t snap(plotter.h5load_rr_timestep(at * n["outfreq"])); - snap *= rhod * 1e3; // water per cubic metre (should be wet density...) & g/kg + arr_t snap(plotter.load_rr_timestep(at * n["outfreq"])); + snap *= 1e3; + plotter.multiply_by_rhod(snap); res_series[plt](at) = blitz::mean(snap); } } @@ -889,7 +790,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { { - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); + arr_t snap(plotter.load_rc_timestep(at * n["outfreq"])); snap *= rhod * 1e3; // water per cubic metre (should be wet density...) & g/kg res_series[plt](at) = blitz::mean(snap); } @@ -902,9 +803,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) 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; + auto tmp = plotter.load_rc_timestep(at * n["outfreq"]) * rhod; + arr_t snap(tmp); + snap += plotter.load_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; @@ -919,8 +820,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { { - auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); + auto tmp = plotter.load_rc_timestep(at * n["outfreq"]); + 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; @@ -935,8 +836,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { { - auto tmp = plotter.h5load_rr_timestep(at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); + auto tmp = plotter.load_rr_timestep(at * n["outfreq"]); + 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; @@ -950,7 +851,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { { - typename Plotter_t::arr_t snap(plotter.h5load_timestep("latent surface flux", at * n["outfreq"], true)); + arr_t snap(plotter.h5load_timestep("latent surface flux", at * n["outfreq"], true)); res_series[plt](at) = blitz::mean(snap); } } @@ -961,27 +862,27 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { { - typename Plotter_t::arr_t snap(plotter.h5load_timestep("sensible surface flux", at * n["outfreq"], true)); + arr_t snap(plotter.h5load_timestep("sensible surface flux", at * n["outfreq"], true)); res_series[plt](at) = blitz::mean(snap); } } catch(...) {if(at==first_timestep) data_found[plt]=0;} } - else if (plt == "er") + else if (plt == "er_dycoms") { //entrainment rate as in the 2009 Ackerman paper // to store total mixingg ratio try { { - auto tmp = plotter.h5load_rc_timestep(at * n["outfreq"]) * 1e3; //g/kg - typename Plotter_t::arr_t snap(tmp); - snap += plotter.h5load_rr_timestep(at * n["outfreq"]) * 1e3; //g/kg + auto tmp = plotter.load_rc_timestep(at * n["outfreq"]) * 1e3; //g/kg + arr_t snap(tmp); + snap += plotter.load_rr_timestep(at * n["outfreq"]) * 1e3; //g/kg rtot = snap; } { auto tmp = plotter.h5load_timestep("rv", at * n["outfreq"]) * 1e3; - typename Plotter_t::arr_t snap(tmp); // vapor mixing ratio [g/kg] + arr_t snap(tmp); // vapor mixing ratio [g/kg] rtot += snap; } plotter.k_i = 0; @@ -996,7 +897,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { auto tmp = plotter.h5load_timestep("w", at * n["outfreq"]); - typename Plotter_t::arr_t snap(tmp); + arr_t snap(tmp); Array mean(n["z"]); snap = snap * snap; // 2nd power, w_mean = 0 // mean variance of w in horizontal @@ -1010,19 +911,19 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t u(plotter.h5load_timestep("u", at * n["outfreq"])); + arr_t u(plotter.h5load_timestep("u", at * n["outfreq"])); plotter.subtract_horizontal_mean(u); u = u * u; res_series[plt](at) = blitz::mean(plotter.horizontal_mean(u)); - typename Plotter_t::arr_t w(plotter.h5load_timestep("w", at * n["outfreq"])); + arr_t w(plotter.h5load_timestep("w", at * n["outfreq"])); plotter.subtract_horizontal_mean(w); w = w * w; res_series[plt](at) += blitz::mean(plotter.horizontal_mean(w)); if (Plotter_t::n_dims > 2) { - typename Plotter_t::arr_t v(plotter.h5load_timestep("v", at * n["outfreq"])); + arr_t v(plotter.h5load_timestep("v", at * n["outfreq"])); plotter.subtract_horizontal_mean(v); v = v * v; res_series[plt](at) += blitz::mean(plotter.horizontal_mean(v)); @@ -1030,8 +931,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) res_series[plt](at) *= 0.5; // * n["dz"]; - typename Plotter_t::arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); - typename Plotter_t::arr_t snap; + arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); + arr_t snap; snap.reference(tke); res_series[plt](at) += blitz::mean(plotter.horizontal_mean(snap)); @@ -1042,19 +943,19 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t u(plotter.h5load_timestep("u", at * n["outfreq"])); + arr_t u(plotter.h5load_timestep("u", at * n["outfreq"])); plotter.subtract_horizontal_mean(u); u = u * u; res_series[plt](at) = blitz::mean(plotter.horizontal_mean(u)); - typename Plotter_t::arr_t w(plotter.h5load_timestep("w", at * n["outfreq"])); + arr_t w(plotter.h5load_timestep("w", at * n["outfreq"])); plotter.subtract_horizontal_mean(w); w = w * w; res_series[plt](at) += blitz::mean(plotter.horizontal_mean(w)); res_series[plt](at) *= 0.5;// * n["dz"]; - typename Plotter_t::arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); + arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); if (Plotter_t::n_dims == 3) { // assume that sgs tke is isotropic, hence 2/3 are in the uw plane @@ -1072,8 +973,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t u(plotter.h5load_timestep("u", at * n["outfreq"])); - typename Plotter_t::arr_t w(plotter.h5load_timestep("w", at * n["outfreq"])); + arr_t u(plotter.h5load_timestep("u", at * n["outfreq"])); + arr_t w(plotter.h5load_timestep("w", at * n["outfreq"])); int step_no = at - first_timestep; if(step_no == 0) // the first step { @@ -1123,7 +1024,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) res_series[plt](at) *= 0.5;// * n["dz"]; - typename Plotter_t::arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); + arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); if (Plotter_t::n_dims == 3) { // assume that sgs tke is isotropic, hence 2/3 are in the uw plane @@ -1140,7 +1041,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); + arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); res_series[plt](at) = blitz::mean(plotter.horizontal_mean(tke)); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -1149,15 +1050,15 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t tot_m0(plotter.h5load_timestep("aerosol_rw_mom0", at * n["outfreq"])); - tot_m0 += typename Plotter_t::arr_t(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - tot_m0 += typename Plotter_t::arr_t(plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"])); + arr_t tot_m0(plotter.h5load_timestep("aerosol_rw_mom0", at * n["outfreq"])); + tot_m0 += arr_t(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); + tot_m0 += arr_t(plotter.h5load_timestep("rain_rw_mom0", at * n["outfreq"])); - typename Plotter_t::arr_t tke(plotter.h5load_timestep("all_up_mom2", at * n["outfreq"])); - tke += typename Plotter_t::arr_t(plotter.h5load_timestep("all_wp_mom2", at * n["outfreq"])); + arr_t tke(plotter.h5load_timestep("all_up_mom2", at * n["outfreq"])); + tke += arr_t(plotter.h5load_timestep("all_wp_mom2", at * n["outfreq"])); if (Plotter_t::n_dims > 2) - tke += typename Plotter_t::arr_t(plotter.h5load_timestep("all_vp_mom2", at * n["outfreq"])); + tke += arr_t(plotter.h5load_timestep("all_vp_mom2", at * n["outfreq"])); tke = blitz::where(tot_m0 > 0., 0.5 * tke / tot_m0, 0); // tke in each cell @@ -1170,13 +1071,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // gccn (r_d > 2 um) concentration in cloudy grid cells 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"])); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - snap = iscloudy(snap); // cloudiness mask - typename Plotter_t::arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"])); - snap2 *= rhod; // b4 it was specific moment + arr_t snap(plotter.get_mask(at * n["outfreq"])); + arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"])); + plotter.multiply_by_rhod(snap2); snap2 /= 1e6; // per cm^3 snap2 *= snap; if(blitz::sum(snap) > 0) @@ -1191,13 +1088,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // gccn (r_d < 2 um) concentration in cloudy grid cells 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"])); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - snap = iscloudy(snap); // cloudiness mask - typename Plotter_t::arr_t snap2(plotter.h5load_timestep("non_gccn_rw_mom0", at * n["outfreq"])); - snap2 *= rhod; // b4 it was specific moment + arr_t snap(plotter.get_mask(at * n["outfreq"])); + arr_t snap2(plotter.h5load_timestep("non_gccn_rw_mom0", at * n["outfreq"])); + plotter.multiply_by_rhod(snap2); snap2 /= 1e6; // per cm^3 snap2 *= snap; if(blitz::sum(snap) > 0) @@ -1211,13 +1104,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { 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"])); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - snap = iscloudy(snap); // cloudiness mask - typename Plotter_t::arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"])); - typename Plotter_t::arr_t snap3(plotter.h5load_timestep("non_gccn_rw_mom0", at * n["outfreq"])); + arr_t snap(plotter.get_mask(at * n["outfreq"])); + arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"])); + arr_t snap3(plotter.h5load_timestep("non_gccn_rw_mom0", at * n["outfreq"])); snap2 = where(snap3 > 0, snap2 / snap3, 0); // even if snap3=0, they are noncloudy anyway snap2 *= snap; @@ -1234,13 +1123,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // gccn (r_d > 2 um) mean radius in cloudy grid cells 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"])); - 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("gccn_rw_mom0", at * n["outfreq"])); - typename Plotter_t::arr_t snap_m1(plotter.h5load_timestep("gccn_rw_mom1", at * n["outfreq"]) * 1e6); // in microns + arr_t snap(plotter.get_mask(at * n["outfreq"])); + arr_t snap_m0(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"])); + arr_t snap_m1(plotter.h5load_timestep("gccn_rw_mom1", at * n["outfreq"]) * 1e6); // in microns snap_m0 *= snap; snap_m1 *= snap; auto tot_gccn_m0 = blitz::sum(snap_m0); @@ -1256,9 +1141,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // gccn (r_d > 2 um) concentration try { - typename Plotter_t::arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"])); + arr_t snap2(plotter.h5load_timestep("gccn_rw_mom0", at * n["outfreq"])); snap2 /= 1e6; // per cm^3 - snap2 *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap2); res_series[plt](at) = blitz::mean(snap2); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -1268,9 +1153,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // gccn (r_d > 2 um) concentration try { - typename Plotter_t::arr_t snap2(plotter.h5load_timestep("non_gccn_rw_mom0", at * n["outfreq"])); + arr_t snap2(plotter.h5load_timestep("non_gccn_rw_mom0", at * n["outfreq"])); snap2 /= 1e6; // per cm^3 - snap2 *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap2); res_series[plt](at) = blitz::mean(snap2); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -1282,13 +1167,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // gccn (r_d >= 0.8 um) concentration in cloudy grid cells 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"])); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - snap = iscloudy(snap); // cloudiness mask - typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_geq_0.8um_rw_mom0", at * n["outfreq"])); - snap2 *= rhod; // b4 it was specific moment + arr_t snap(plotter.get_mask(at * n["outfreq"])); + arr_t snap2(plotter.h5load_timestep("rd_geq_0.8um_rw_mom0", at * n["outfreq"])); + plotter.multiply_by_rhod(snap2); snap2 /= 1e6; // per cm^3 snap2 *= snap; if(blitz::sum(snap) > 0) @@ -1303,9 +1184,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // gccn (r_d >= 0.8 um) concentration try { - typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_geq_0.8um_rw_mom0", at * n["outfreq"])); + arr_t snap2(plotter.h5load_timestep("rd_geq_0.8um_rw_mom0", at * n["outfreq"])); snap2 /= 1e6; // per cm^3 - snap2 *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap2); res_series[plt](at) = blitz::mean(snap2); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -1316,13 +1197,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { 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"])); - snap *= rhod; // b4 it was specific moment - snap /= 1e6; // per cm^3 - snap = iscloudy(snap); // cloudiness mask - typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_lt_0.8um_rw_mom0", at * n["outfreq"])); - snap2 *= rhod; // b4 it was specific moment + arr_t snap(plotter.get_mask(at * n["outfreq"])); + arr_t snap2(plotter.h5load_timestep("rd_lt_0.8um_rw_mom0", at * n["outfreq"])); + plotter.multiply_by_rhod(snap2); snap2 /= 1e6; // per cm^3 snap2 *= snap; if(blitz::sum(snap) > 0) @@ -1336,55 +1213,29 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t snap2(plotter.h5load_timestep("rd_lt_0.8um_rw_mom0", at * n["outfreq"])); + arr_t snap2(plotter.h5load_timestep("rd_lt_0.8um_rw_mom0", at * n["outfreq"])); snap2 /= 1e6; // per cm^3 - snap2 *= rhod; // b4 it was per milligram + plotter.multiply_by_rhod(snap2); res_series[plt](at) = blitz::mean(snap2); } catch(...) {if(at==first_timestep) data_found[plt]=0;} } - - else if (plt == "cl_meanr") - { - // cloud droplets mean radius in cloudy grid cells - 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"])); - 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 - snap_m0 *= snap; - snap_m1 *= snap; - auto tot_m0 = blitz::sum(snap_m0); - if(tot_m0 > 0) - res_series[plt](at) = blitz::sum(snap_m1) / tot_m0; - else - res_series[plt](at) = 0; - } - catch(...) {if(at==first_timestep) data_found[plt]=0;} - } // cloud base mean incloud time of bigrain (r>40um) else if (plt == "clb_bigrain_mean_inclt") { try { - // find cloud base (cloudy if q_c > 0.1 g/kg) - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - snap = iscloudy_rc(snap); // cloudiness mask - //for(int i=0;i<10;++i) - // snap(plotter.hrzntl_slice(i)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux + // find cloud base + arr_t snap(plotter.get_mask(at * n["outfreq"])); plotter.k_i = blitz::first((snap == 1), plotter.LastIndex); // 0-th specific mom of bigrain cloud drops - typename Plotter_t::arr_t bigrain_conc(plotter.h5load_timestep("bigrain_rw_mom0", at * n["outfreq"])); + arr_t bigrain_conc(plotter.h5load_timestep("bigrain_rw_mom0", at * n["outfreq"])); // concentration of bigrain cloud drops at cloud base plotter.tmp_float_hrzntl_slice = plotter.get_value_at_hgt(bigrain_conc, plotter.k_i) * plotter.get_value_at_hgt(rhod, plotter.k_i); // we need to multiply by rhod here, because different cloud bases can mean different rhod // 1st specific mom of incloud time of bigrain drops - typename Plotter_t::arr_t bigrain_inclt_mom1(plotter.h5load_timestep("bigrain_incloud_time_mom1", at * n["outfreq"])); + arr_t bigrain_inclt_mom1(plotter.h5load_timestep("bigrain_incl_time_mom1", at * n["outfreq"])); // 1st mom of incloud time at cloud base plotter.tmp_float_hrzntl_slice2 = plotter.get_value_at_hgt(bigrain_inclt_mom1, plotter.k_i) * plotter.get_value_at_hgt(rhod, plotter.k_i); // same as above @@ -1401,19 +1252,18 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { // find cloud base (cloudy if q_c > 0.1 g/kg) - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - snap = iscloudy_rc(snap); // cloudiness mask + arr_t snap(plotter.get_mask(at * n["outfreq"])); //for(int i=0;i<10;++i) // snap(plotter.hrzntl_slice(i)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux plotter.k_i = blitz::first((snap == 1), plotter.LastIndex); // 0-th specific mom of bigrain cloud drops - typename Plotter_t::arr_t bigrain_conc(plotter.h5load_timestep("bigrain_rw_mom0", at * n["outfreq"])); + arr_t bigrain_conc(plotter.h5load_timestep("bigrain_rw_mom0", at * n["outfreq"])); // concentration of bigrain cloud drops at cloud base plotter.tmp_float_hrzntl_slice = plotter.get_value_at_hgt(bigrain_conc, plotter.k_i) * plotter.get_value_at_hgt(rhod, plotter.k_i); // we need to multiply by rhod here, because different cloud bases can mean different rhod // 1st specific mom of rd of bigrain drops - typename Plotter_t::arr_t bigrain_rd_mom1(plotter.h5load_timestep("bigrain_rd_mom1", at * n["outfreq"])); + arr_t bigrain_rd_mom1(plotter.h5load_timestep("bigrain_rd_mom1", at * n["outfreq"])); // 1st mom of rd at cloud base plotter.tmp_float_hrzntl_slice2 = plotter.get_value_at_hgt(bigrain_rd_mom1, plotter.k_i) * plotter.get_value_at_hgt(rhod, plotter.k_i); // same as above @@ -1430,19 +1280,18 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { // find cloud base (cloudy if q_c > 0.1 g/kg) - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - snap = iscloudy_rc(snap); // cloudiness mask + arr_t snap(plotter.get_mask(at * n["outfreq"])); //for(int i=0;i<10;++i) // snap(plotter.hrzntl_slice(i)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux plotter.k_i = blitz::first((snap == 1), plotter.LastIndex); // 0-th specific mom of bigrain cloud drops - typename Plotter_t::arr_t bigrain_conc(plotter.h5load_timestep("bigrain_rw_mom0", at * n["outfreq"])); + arr_t bigrain_conc(plotter.h5load_timestep("bigrain_rw_mom0", at * n["outfreq"])); // concentration of bigrain cloud drops at cloud base plotter.tmp_float_hrzntl_slice = plotter.get_value_at_hgt(bigrain_conc, plotter.k_i) * plotter.get_value_at_hgt(rhod, plotter.k_i); // we need to multiply by rhod here, because different cloud bases can mean different rhod // 1st specific mom of rd of bigrain drops - typename Plotter_t::arr_t bigrain_kappa_mom1(plotter.h5load_timestep("bigrain_kappa_mom1", at * n["outfreq"])); + arr_t bigrain_kappa_mom1(plotter.h5load_timestep("bigrain_kappa_mom1", at * n["outfreq"])); // 1st mom of rd at cloud base plotter.tmp_float_hrzntl_slice2 = plotter.get_value_at_hgt(bigrain_kappa_mom1, plotter.k_i) * plotter.get_value_at_hgt(rhod, plotter.k_i); // same as above @@ -1459,14 +1308,13 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { // find cloud base (cloudy if q_c > 0.1 g/kg) - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - snap = iscloudy_rc(snap); // cloudiness mask + arr_t snap(plotter.get_mask(at * n["outfreq"])); //for(int i=0;i<10;++i) // snap(plotter.hrzntl_slice(i)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux plotter.k_i = blitz::first((snap == 1), plotter.LastIndex); // 0-th specific mom of bigrain cloud drops - typename Plotter_t::arr_t bigrain_conc(plotter.h5load_timestep("bigrain_rw_mom0", at * n["outfreq"])); + arr_t bigrain_conc(plotter.h5load_timestep("bigrain_rw_mom0", at * n["outfreq"])); // concentration of bigrain cloud drops at cloud base [1/m^3] plotter.tmp_float_hrzntl_slice = plotter.get_value_at_hgt(bigrain_conc, plotter.k_i) * plotter.get_value_at_hgt(rhod, plotter.k_i); // we need to multiply by rhod here, because different cloud bases can mean different rhod @@ -1485,19 +1333,18 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { // find cloud base (cloudy if q_c > 0.1 g/kg) - typename Plotter_t::arr_t snap(plotter.h5load_rc_timestep(at * n["outfreq"])); - snap = iscloudy_rc(snap); // cloudiness mask + arr_t snap(plotter.get_mask(at * n["outfreq"])); //for(int i=0;i<10;++i) // snap(plotter.hrzntl_slice(i)) = 0; // cheat to avoid occasional "cloudy" cell at ground level due to activation from surf flux plotter.k_i = blitz::first((snap == 1), plotter.LastIndex); // 0-th specific mom of bigrain cloud drops - typename Plotter_t::arr_t bigrain_conc(plotter.h5load_timestep("bigrain_rw_mom0", at * n["outfreq"])); + arr_t bigrain_conc(plotter.h5load_timestep("bigrain_rw_mom0", at * n["outfreq"])); // concentration of bigrain cloud drops at cloud base [1/m^3] plotter.tmp_float_hrzntl_slice = plotter.get_value_at_hgt(bigrain_conc, plotter.k_i) * plotter.get_value_at_hgt(rhod, plotter.k_i); // we need to multiply by rhod here, because different cloud bases can mean different rhod // 0-th specific mom of bigrain cloud drops formed on gccn - typename Plotter_t::arr_t bigrain_gccn_conc(plotter.h5load_timestep("bigrain_gccn_rw_mom0", at * n["outfreq"])); + arr_t bigrain_gccn_conc(plotter.h5load_timestep("bigrain_gccn_rw_mom0", at * n["outfreq"])); // concentration of bigrain cloud drops at cloud base [1/m^3] plotter.tmp_float_hrzntl_slice2 = plotter.get_value_at_hgt(bigrain_gccn_conc, plotter.k_i) * plotter.get_value_at_hgt(rhod, plotter.k_i); // we need to multiply by rhod here, because different cloud bases can mean different rhod @@ -1545,8 +1392,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t lwc(plotter.h5load_rc_timestep(at * n["outfreq"])); // cloud water, no rain water in pi chamber icmw - //res_series[plt](at) = blitz::mean(typename Plotter_t::arr_t(plotter.nowall(lwc, distance_from_walls))); + arr_t lwc(plotter.load_rc_timestep(at * n["outfreq"])); // cloud water, no rain water in pi chamber icmw + //res_series[plt](at) = blitz::mean(arr_t(plotter.nowall(lwc, distance_from_walls))); res_series[plt](at) = blitz::mean(lwc); } catch(...) {if(at==first_timestep) data_found[plt]=0;} @@ -1556,7 +1403,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t lwc(plotter.h5load_rc_timestep(at * n["outfreq"])); // cloud water, no rain water in pi chamber icmw + arr_t lwc(plotter.load_rc_timestep(at * n["outfreq"])); // cloud water, no rain water in pi chamber icmw lwc *= rhod; res_series[plt](at) = blitz::mean(lwc); } @@ -1577,7 +1424,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t nc(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); + arr_t nc(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); nc *= rhod; // 1/kg -> 1/m^3 res_series[plt](at) = blitz::mean(nc); } @@ -1588,7 +1435,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t na(plotter.h5load_timestep("aerosol_rw_mom0", at * n["outfreq"])); + arr_t na(plotter.h5load_timestep("aerosol_rw_mom0", at * n["outfreq"])); na *= rhod; // 1/kg -> 1/m^3 res_series[plt](at) = blitz::mean(na); } @@ -1599,19 +1446,19 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t u(plotter.h5load_timestep("u", at * n["outfreq"])); + arr_t u(plotter.h5load_timestep("u", at * n["outfreq"])); plotter.subtract_horizontal_mean(u); u = u * u; res_series[plt](at) = blitz::mean(u); - typename Plotter_t::arr_t w(plotter.h5load_timestep("w", at * n["outfreq"])); + arr_t w(plotter.h5load_timestep("w", at * n["outfreq"])); plotter.subtract_horizontal_mean(w); w = w * w; res_series[plt](at) += blitz::mean(w); if (Plotter_t::n_dims > 2) { - typename Plotter_t::arr_t v(plotter.h5load_timestep("v", at * n["outfreq"])); + arr_t v(plotter.h5load_timestep("v", at * n["outfreq"])); plotter.subtract_horizontal_mean(v); v = v * v; res_series[plt](at) += blitz::mean(v); @@ -1625,19 +1472,19 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { try { - typename Plotter_t::arr_t u(plotter.h5load_timestep("u", at * n["outfreq"])); + arr_t u(plotter.h5load_timestep("u", at * n["outfreq"])); plotter.subtract_horizontal_mean(u); u = u * u; res_series[plt](at) = blitz::mean(u); - typename Plotter_t::arr_t w(plotter.h5load_timestep("w", at * n["outfreq"])); + arr_t w(plotter.h5load_timestep("w", at * n["outfreq"])); plotter.subtract_horizontal_mean(w); w = w * w; res_series[plt](at) += blitz::mean(w); if (Plotter_t::n_dims > 2) { - typename Plotter_t::arr_t v(plotter.h5load_timestep("v", at * n["outfreq"])); + arr_t v(plotter.h5load_timestep("v", at * n["outfreq"])); plotter.subtract_horizontal_mean(v); v = v * v; res_series[plt](at) += blitz::mean(v); @@ -1645,8 +1492,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) res_series[plt](at) *= 0.5; // * n["dz"]; - typename Plotter_t::arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); - typename Plotter_t::arr_t snap; + arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); + arr_t snap; snap.reference(tke); res_series[plt](at) += blitz::mean(plotter.horizontal_mean(snap)); @@ -1658,9 +1505,9 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // droplets mean radius away from walls try { - typename Plotter_t::arr_t m0(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); - typename Plotter_t::arr_t m1(plotter.h5load_timestep("cloud_rw_mom1", at * n["outfreq"])); - //auto tot_m0 = blitz::sum(typename Plotter_t::arr_t(plotter.nowall(m0, distance_from_walls))); + arr_t m0(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); + arr_t m1(plotter.h5load_timestep("cloud_rw_mom1", at * n["outfreq"])); + //auto tot_m0 = blitz::sum(arr_t(plotter.nowall(m0, distance_from_walls))); auto tot_m0 = blitz::sum(m0); if(tot_m0 > 0) res_series[plt](at) = blitz::sum(m1) / tot_m0; @@ -1674,8 +1521,8 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // droplets effective radius away from walls try { - typename Plotter_t::arr_t m2(plotter.h5load_timestep("cloud_rw_mom2", at * n["outfreq"])); - typename Plotter_t::arr_t m3(plotter.h5load_timestep("cloud_rw_mom3", at * n["outfreq"])); + arr_t m2(plotter.h5load_timestep("cloud_rw_mom2", at * n["outfreq"])); + arr_t m3(plotter.h5load_timestep("cloud_rw_mom3", at * n["outfreq"])); auto tot_m2 = blitz::sum(m2); if(tot_m2 > 0) res_series[plt](at) = blitz::sum(m3) / tot_m2; @@ -1729,10 +1576,10 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) // relative dispersion (std dev / mean) of droplet radius distribution averaged over cells away from walls try { - //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("cloud_rw_mom0", at * n["outfreq"])); - typename Plotter_t::arr_t m1(plotter.h5load_timestep("cloud_rw_mom1", at * n["outfreq"])); - typename Plotter_t::arr_t m2(plotter.h5load_timestep("cloud_rw_mom2", at * n["outfreq"])); + //arr_t m0(plotter.nowall(arr_t(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])), distance_from_walls)); + arr_t m0(plotter.h5load_timestep("cloud_rw_mom0", at * n["outfreq"])); + arr_t m1(plotter.h5load_timestep("cloud_rw_mom1", at * n["outfreq"])); + arr_t m2(plotter.h5load_timestep("cloud_rw_mom2", at * n["outfreq"])); // calculate stddev of radius, store in m2 m2 = where(m0 > 0, m2 / m0 - m1 / m0 * m1 / m0, 0.); @@ -1758,7 +1605,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) try { const float C_E = 0.845; - typename Plotter_t::arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); + arr_t tke(plotter.h5load_timestep("tke", at * n["outfreq"])); tke = pow(tke, 3./2.); tke /= n_prof["mix_len"](plotter.LastIndex); // divide by SGS mixing length // res_series[plt](at) = blitz::mean(plotter.nowall(tke, distance_from_walls)) * C_E; @@ -1836,7 +1683,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) res_series[plt] /= 1000.; res_pos *= 60.; } - else if (plt == "ract_avg") + else if (plt == "cl_ract") { plot_std_dev = true; res_pos *= 60.; @@ -1845,30 +1692,62 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { res_pos *= 60.; } - else if (plt == "cloud_avg_act_conc") + else if (plt == "cl_nr") + { + plot_std_dev = true; + res_pos *= 60.; + } + else if (plt == "cl_nc") + { + plot_std_dev = true; + res_pos *= 60.; + } + else if (plt == "cl_avg_act_conc") + { + plot_std_dev = true; + res_pos *= 60.; + } + else if (plt == "cl_std_dev_act_conc") + { + res_pos *= 60.; + } + else if (plt == "cl_avg_supersat") + { + plot_std_dev = true; + res_pos *= 60.; + } + else if (plt == "cl_avg_th") + { + plot_std_dev = true; + res_pos *= 60.; + } + else if (plt == "cl_avg_rv") { plot_std_dev = true; res_pos *= 60.; } - else if (plt == "cloud_std_dev_act_conc") + else if (plt == "cl_std_dev_supersat") { res_pos *= 60.; } - else if (plt == "cloud_avg_supersat") + else if (plt == "cl_avg_act_meanr") { plot_std_dev = true; res_pos *= 60.; } - else if (plt == "cloud_std_dev_supersat") + else if (plt == "cl_avg_act_stddevr") { + plot_std_dev = true; res_pos *= 60.; } - else if (plt == "cloud_avg_act_rad") + else if (plt == "cl_avg_cloud_meanr") { + plot_std_dev = true; res_pos *= 60.; } - else if (plt == "cloud_avg_std_dev_act_rad") + else if (plt == "cl_avg_cloud_stddevr") { + plot_std_dev = true; res_pos *= 60.; } else if (plt == "sd_conc") @@ -1923,7 +1802,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) { res_pos *= 60.; } - else if (plt == "cloud_top_height") + else if (plt == "cl_top_height") { res_pos *= 60.; } @@ -1933,17 +1812,17 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) } else if (plt == "lwp") { - res_series[plt] *= (n["z"] - 1) * n["dz"]; // top and bottom cells are smaller + res_series[plt] *= (n["refined z"] - 1) * n["refined dz"]; // top and bottom cells are smaller } else if (plt == "rwp") { - res_series[plt] *= (n["z"] - 1) * n["dz"]; // top and bottom cells are smaller + res_series[plt] *= (n["refined z"] - 1) * n["refined dz"]; // top and bottom cells are smaller } else if (plt == "cwp") { res_series[plt] *= (n["z"] - 1) * n["dz"]; // top and bottom cells are smaller } - else if (plt == "er") + else if (plt == "er_dycoms") { // central difference, in cm Range nofirstlast = Range(1, last_timestep-1); @@ -2028,7 +1907,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) gp << "plot '-' with l"; if(plot_std_dev) gp << ", '-' w l, '-' w l"; - else if(plt == "cl_acnv25_dycoms" || plt == "cl_acnv25_rico" || plt == "cl_accr25_dycoms" || plt == "cl_accr25_rico") + else if(plt == "cl_acnv25" || plt == "cl_accr25" ) gp << ", '-' w l"; gp << " \n"; @@ -2044,7 +1923,7 @@ void plot_series(Plotter_t plotter, Plots plots, std::string type) res_series[plt] = res_series[plt] - 2*res_series_std_dev[plt]; gp.send1d(boost::make_tuple(res_pos, res_series[plt])); } - if(plt == "cl_acnv25_dycoms" || plt == "cl_acnv25_rico" || plt == "cl_accr25_dycoms" || plt == "cl_accr25_rico") + if(plt == "cl_acnv25" || plt == "cl_accr25") { // acnv/accr rate averaged since the start of the simulation int nt = last_timestep - first_timestep + 1; diff --git a/histograms/histogram.py b/histograms/histogram.py index 430512a..2adf2f3 100644 --- a/histograms/histogram.py +++ b/histograms/histogram.py @@ -2,31 +2,49 @@ import argparse import h5py import numpy as np -from sys import argv +from sys import path import matplotlib.pyplot as plt from collections import OrderedDict -c_pd = 1005.7 -R_d = 287. +path.append('/home/piotr/singu_built_libraries/usr/lib/python3/dist-packages/') +from libcloudphxx import common -def exner(p): - return (p / 1e5)**(R_d/c_pd) -v_exner = np.vectorize(exner) +#c_pd = 1005.7 +#R_d = 287. +# +#def exner(p): +# return (p / 1e5)**(R_d/c_pd) +#v_exner = np.vectorize(exner) +# def calc_T(th, p): - return th * v_exner(p) -v_calc_T = np.vectorize(calc_T) + return np.float64(th) * common.exner(np.float64(p)) +#v_calc_T = np.vectorize(calc_T) -# Tetens: r_vs=3.8/(p*exp(-17.2693882*(T-273.15)/(T-35.86))-6.109) p in mb, T in Kelvins -def calc_rv_s(th, rv, p): - T = v_calc_T(th, p) - return 3.8 / (p*np.exp(-17.2693882*(T-273.15)/(T-35.86))-6.109) -v_calc_rv_s = np.vectorize(calc_rv_s) +# +## Tetens: r_vs=3.8/(p*exp(-17.2693882*(T-273.15)/(T-35.86))-6.109) p in mb, T in Kelvins +#def calc_rv_s(th, rv, p): +# T = v_calc_T(th, p) +# return 3.8 / (p*np.exp(-17.2693882*(T-273.15)/(T-35.86))-6.109) + +def calc_r_vs_T(T, p): + return common.r_vs(np.float64(T), np.float64(p)) +v_calc_r_vs_T = np.vectorize(calc_r_vs_T) + +def calc_r_vs(th, p): + T = calc_T(th, p) + return common.r_vs(T, p) +v_calc_r_vs = np.vectorize(calc_r_vs) def calc_RH(th, rv, p): - return rv / v_calc_rv_s(th, rv, p) / 100. + return rv / v_calc_r_vs(th, p) v_calc_RH = np.vectorize(calc_RH) +def calc_th_l(th, rl, p): + T = calc_T(th, p) + return th - 1. / common.exner(np.float64(p)) * common.l_v(np.float64(T)) / common.c_pd * rl +v_calc_th_l = np.vectorize(calc_th_l) + #mpl.rcParams['figure.figsize'] = 10, 10 plt.rcParams.update({'font.size': 10}) plt.figure(1, figsize=(10,10)) # histogram plot @@ -45,9 +63,11 @@ def calc_RH(th, rv, p): parser.add_argument("--outfreq", type=int, required=False, help="output frequency of the simulation [number of time steps], if not specified it will be read from const.h5 (if possible)") parser.add_argument('--normalize', action='store_true', help="normalize the histogram") parser.add_argument('--no_histogram', action='store_true', help="dont save the histogram plot") -parser.add_argument('--no_correlations', action='store_true', help="dont calculate correlations and dot save the scatter plot") +parser.add_argument('--no_scatter', action='store_true', help="dont calculate correlations and dot save the scatter plot") parser.add_argument('--mask_rico', action='store_true', help="compute histogram only within cloud cells (using the rico cloud mask)") -parser.set_defaults(normalie=False) +parser.add_argument('--mask_rico_rl', action='store_true', help="compute histogram only within cloud cells (using the rico cloud mask based on the r_l array)") +parser.add_argument('--scatter_saturation', action='store_true', help="plot saturation line on th vs rv scatter plots") +parser.set_defaults(normalize=False, scatter_saturation=False) args = parser.parse_args() @@ -58,6 +78,8 @@ def calc_RH(th, rv, p): dz = {} ref = {} +scatter_saturation_plotted = False + # directories loop for directory, lab in zip(args.dirs, args.labels): @@ -104,14 +126,21 @@ def calc_RH(th, rv, p): # init variable-specific array parameters based on the first timestep filename = directory + "/timestep" + str(time_start_idx).zfill(10) + ".h5" - # special case of RH calculated from th and rv - # NOTE: RH_derived is calculated with r_v / r_vs, where r_vs comes from the Tetens formula - # in UWLCM r_vs comes from clausius-clapeyron. - # RH_derived is shifted to the right with respect to RH from UWLCM - maybe due to this difference? + # special cases of variables that need processing first if(var == "RH_derived"): w3d = h5py.File(filename, "r")["th"][:,:,:] elif(can_plot_refined_RH_derived and var == "refined RH_derived"): w3d = h5py.File(filename, "r")["refined th"][:,:,:] + elif(var == "r_tot"): # total resolved water mixing ratio + try: + w3d = h5py.File(filename, "r")["rv"][:,:,:] + except: + continue + elif(var == "th_l"): # liquid water potential temperature + try: + w3d = h5py.File(filename, "r")["th"][:,:,:] + except: + continue else: try: w3d = h5py.File(filename, "r")[var][:,:,:] @@ -144,16 +173,28 @@ def calc_RH(th, rv, p): lab_var = lab + '_' + str(var) + assert(not (args.mask_rico and args.mask_rico_rl)) + if(args.mask_rico): try: mask = h5py.File(filename, "r")["cloud_rw_mom3"][:,:,:] if(mask.shape != w3d.shape): - print("Cloud mask shape is different than "+var+" shape. Skipping the plot.") + print("Cloud mask shape: " + str(mask.shape) + " is different than "+var+" shape: " + str(w3d.shape) + ". Skipping the plot.") continue except: print("Can't find cloud_rw_mom3 data, so can't use RICO cloud mask. Skipping the plot.") continue + if(args.mask_rico_rl): + try: + mask = h5py.File(filename, "r")["r_l"][:,:,:] + if(mask.shape != w3d.shape): + print("Cloud mask shape: " + str(mask.shape) + " is different than "+var+" shape: " + str(w3d.shape) + ". Skipping the plot.") + continue + except: + print("Can't find r_l data, so can't use RICO r_l cloud mask. Skipping the plot.") + continue + total_arr[lab_var] = np.zeros(0) plot_labels[lab_var] = lab_var @@ -182,6 +223,18 @@ def calc_RH(th, rv, p): p[i,j] = refined_p_e[level_start_idx:level_end_idx] w3d = v_calc_RH(th, rv, p) + elif(var == "r_tot"): # total resolved water mixing ratio + w3d = h5py.File(filename, "r")["rv"][:,:,level_start_idx:level_end_idx] + h5py.File(filename, "r")["r_l"][:,:,level_start_idx:level_end_idx] + + elif(var == "th_l"): # liquid water potential temp + th = h5py.File(filename, "r")["th"][:,:,level_start_idx:level_end_idx] + rl = h5py.File(filename, "r")["r_l"][:,:,level_start_idx:level_end_idx] + p = np.empty(th.shape) + for i in np.arange(p.shape[0]): + for j in np.arange(p.shape[1]): + p[i,j] = p_e[level_start_idx:level_end_idx] + w3d = v_calc_th_l(th, rl, p) + else: w3d = h5py.File(filename, "r")[var][0:nx-1, 0:ny-1, level_start_idx:level_end_idx] # * 4. / 3. * 3.1416 * 1e3 @@ -191,6 +244,12 @@ def calc_RH(th, rv, p): mask = np.where(mask > 1.e-5, 1., 0.) w3d = w3d[(mask == 1)] + # read and apply cloud mask + if(args.mask_rico_rl): + mask = h5py.File(filename, "r")["r_l"][0:nx, 0:ny, level_start_idx:level_end_idx] + mask = np.where(mask > 1.e-5, 1., 0.) + w3d = w3d[(mask == 1)] + total_arr[lab_var] = np.append(total_arr[lab_var], w3d) @@ -217,9 +276,21 @@ def calc_RH(th, rv, p): plt.scatter(total_arr[lab_var1].flatten(), total_arr[lab_var2].flatten(), marker='.', s=1, label=lab) plt.xlabel(var1) plt.ylabel(var2) - # plt.clf() - + # add saturation line on a rv vs th plot + # TODO: add this line also to refined rv vs refined th plots + if(args.scatter_saturation and var1 == "th" and var2 == "rv" and args.level_start == args.level_end and not scatter_saturation_plotted): + print("plotting a saturation line on the scatter plot") + press = p_e[level_start_idx] + min_th = total_arr[lab_var1].min() + max_th = total_arr[lab_var1].max() + min_T = calc_T(min_th, press) + max_T = calc_T(max_th, press) + v_th = np.linspace(min_th, max_th, 100) + v_T = np.linspace(min_T, max_T, 100) + v_r_vs = v_calc_r_vs_T(v_T, press) + plt.plot(v_th, v_r_vs, c='black', ls='--', label='r_vs') + scatter_saturation_plotted = True @@ -271,5 +342,8 @@ def calc_RH(th, rv, p): plt.figure(2) plt.legend()#loc = 'lower center') -if(not args.no_correlations): - plt.savefig(args.outfig + 'scatter.svg') +plt.grid(True, which='both', linestyle='--') +plt.title("z=["+str(args.level_start)+"m, "+str(args.level_end)+"m] @["+str(args.time_start)+"s, "+str(args.time_end)+"s]") + +if(not args.no_scatter): + plt.savefig(args.outfig + 'scatter.png') # .svg