Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

--sd_const_multi #104

Open
wants to merge 13 commits into
base: master
Choose a base branch
from
25 changes: 16 additions & 9 deletions plotting/Energy_spectrum/spectrum.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,16 @@
import h5py
import numpy as np
from sys import argv
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt


#velocities = ["u", "v", "w"]
#velocities = ["u", "v", "w", "cloud_rw_mom3", "rv", "th", "RH", "aerosol_rw_mom3"]
velocities = ["u", "cloud_rw_mom3", "w"]
#velocities = ["cloud_rw_mom3"]
velocities = ["w"]
#velocities = ["w"]

time_start = int(argv[1])
time_end = int(argv[2])
Expand Down Expand Up @@ -53,23 +57,26 @@
Exy_avg[vel] += Exy

K = np.fft.rfftfreq(nx - 1)
print nx
print K
# plt.loglog(K, Exy)
lmbd = 50. / K # assume dx=50m
# lmbd = 50. / K # assume dx=50m

if (t == time_start and lab==labels[0]):
plt.loglog(lmbd, 2e-1* K**(-5./3.) )

# if (t == time_start and lab==labels[0]):
# plt.loglog(lmbd, 2e-1* K**(-5./3.) )

f1=plt.figure(1)
for vel in velocities:
Exy_avg[vel] /= (time_end - time_start) / outfreq + 1
Exy_avg[vel] /= to_lvl+1 - from_lvl
#crudely scale
# Exy_avg[vel] /= Exy_avg[vel][len(Exy_avg[vel])-1]
plt.loglog(lmbd, Exy_avg[vel] , linewidth=2, label=lab+"_"+vel)

plt.xlim(10**4,10**2)
# plt.loglog(lmbd, Exy_avg[vel] , linewidth=2, label=lab+"_"+vel)
plt.loglog(K, Exy_avg[vel] , linewidth=2, label=lab+"_"+vel)
#plt.xlim(10**4,10**2)
plt.xlabel("l[m]")
plt.ylabel("PSD")
plt.legend()
plt.grid(True, which='both', linestyle='--')
plt.title("Mean PSD of w 322m<z<642m @3h")
plt.show()
f1.savefig('spectrum.png')
74 changes: 74 additions & 0 deletions plotting/Lasher_Trapp/Lasher_Trapp_series_comparison.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,74 @@
from matplotlib import rc
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator, MultipleLocator

import os
import sys
sys.path.append(os.path.dirname(os.path.abspath(__file__)) + "/../Matplotlib_common/")

from plot_ranges import xscaledict, yscaledict, xlimdict_series, ylimdict_series
from plot_series import *

# activate latex text rendering
rc('text', usetex=True)

Lasher_Trapp_vars = ["lwp", "rwp", "surf_precip", "acc_precip", "cl_nc"]

# init the plot
nplotx = 2
nploty= 3
fig, axarr = plt.subplots(nplotx,nploty)

plot_series( Lasher_Trapp_vars, 0, nplotx, nploty, axarr, xscaledict, yscaledict, xlimdict_series, ylimdict_series, xlabel='Time [h]')

# legend font size
plt.rcParams.update({'font.size': 8})

# hide axes on empty plots
if len( Lasher_Trapp_vars) % nploty == 0:
nemptyplots = 0
else:
nemptyplots = nploty - len( Lasher_Trapp_vars) % nploty
emptyplots = np.arange(nploty - nemptyplots, nploty)
for empty in emptyplots:
axarr[nplotx-1, empty].axis('off')

# hide hrzntl tic labels
x_empty_label = np.arange(0, nplotx-1)
y_empty_label = np.arange(nploty)
for x in x_empty_label:
for y in y_empty_label:
axarr[x,y].set_xticklabels([])

#axes = plt.gca()
#axes.tick_params(direction='in')
x_arr = np.arange(nplotx)
y_arr = np.arange(nploty)
for x in x_arr:
for y in y_arr:
#tics inside
axarr[x,y].tick_params(direction='in', which='both', top=1, right=1)
#minor tics
axarr[x,y].xaxis.set_minor_locator(AutoMinorLocator())
axarr[x,y].yaxis.set_minor_locator(AutoMinorLocator())
#labels and tics font size
for item in ([axarr[x,y].xaxis.label, axarr[x,y].yaxis.label] + axarr[x,y].get_xticklabels() + axarr[x,y].get_yticklabels()):
item.set_fontsize(8)
# subplot numbering
if y < nploty - nemptyplots or x < (nplotx - 1):
axarr[x,y].text(0.2, 0.9, labeldict[y + x*nploty], fontsize=8, transform=axarr[x,y].transAxes)

#single legend for the whole figure
handles, labels = axarr[0,0].get_legend_handles_labels()

lgd = fig.legend(handles, labels, handlelength=4, loc='lower center', bbox_to_anchor=(0.45,0))

#figure size
fig.set_size_inches(7.874, 5. + (len(labels) - 2) * 0.2)# 5.214)#20.75,13.74)
#distances between subplots and from bottom of the plot
fig.subplots_adjust(bottom=0.18 + (len(labels) - 2) * 0.03, hspace=0.1, wspace=0.4)

#fig.tight_layout(pad=0, w_pad=0, h_pad=0)

#plt.show()
fig.savefig(argv[len(sys.argv)-1], bbox_inches='tight', dpi=300)#, bbox_extra_artists=(lgd,))
105 changes: 105 additions & 0 deletions plotting/Lasher_Trapp/plot_ranges.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,105 @@
xscaledict = {
"thl" : "linear",
"00rtot" : "linear",
"rliq" : "linear",
"prflux" : "linear",
"cl_nc" : "linear",
"clfrac" : "linear",
"wvar" : "linear",
"w3rd" : "linear",
"sat_RH" : "linear",
"rad_flx" : "linear",
"lwp" : "linear",
"rwp" : "linear",
"er" : "linear",
"wvarmax" : "linear",
"surf_precip" : "linear",
"acc_precip" : "linear",
"cloud_base" : "linear",
"gccn_rw_cl" : "linear",
"non_gccn_rw_cl" : "linear",
"base_prflux_vs_clhght" : "log",
"cl_gccn_conc" : "linear"
}

yscaledict = {
"thl" : "linear",
"00rtot" : "linear",
"rliq" : "linear",
"prflux" : "linear",
"cl_nc" : "linear",
"clfrac" : "linear",
"wvar" : "linear",
"w3rd" : "linear",
"sat_RH" : "linear",
"rad_flx" : "linear",
"lwp" : "linear",
"rwp" : "linear",
"er" : "linear",
"wvarmax" : "linear",
"surf_precip" : "linear",
"acc_precip" : "linear",
"cloud_base" : "linear",
"gccn_rw_cl" : "linear",
"non_gccn_rw_cl" : "linear",
"base_prflux_vs_clhght" : "linear",
"cl_gccn_conc" : "log"
}

xlimdict_profs = {
"thl" : None,
"00rtot" : None,
"rliq" : None,
"prflux" : (0,20),
"cl_nc" : (0,90),
"clfrac" : None,
"wvar" : None,
"w3rd" : None,
"sat_RH" : None,
"rad_flx" : None,
"gccn_rw_cl" : (0,90),
"non_gccn_rw_cl" : (0,12),
"base_prflux_vs_clhght" : (1,10000)
}

ylimdict_profs = {
"thl" : (0,3000),
"00rtot" : (0,3000),
"rliq" : (0,3000),
"prflux" : (0,3000),
"cl_nc" : (0,3000),
"clfrac" : (0,3000),
"wvar" : (0,3000),
"w3rd" : (0,3000),
"sat_RH" : (0,3000),
"rad_flx" : (0,3000),
"gccn_rw_cl" : (0,3000),
"non_gccn_rw_cl" : (0,3000),
"base_prflux_vs_clhght" : (0,2500)
}

xlimdict_series = {
"clfrac" : None,
"cl_nc" : None,
"lwp" : None,
"rwp" : None,
"er" : None,
"wvarmax" : None,
"surf_precip" : None,
"acc_precip" : None,
"cl_gccn_conc" : None,
"cloud_base" : None
}

ylimdict_series = {
"clfrac" : None,
"cl_nc" : None,
"lwp" : None,
"rwp" : None,
"er" : None,
"wvarmax" : None,
"surf_precip" : None,
"acc_precip" : None,#(0,0.07),
"cl_gccn_conc" : (1e-6, 1),
"cloud_base" : None
}
Binary file added plotting/Lasher_Trapp/plot_ranges.pyc
Binary file not shown.
10 changes: 5 additions & 5 deletions plotting/drawbicyc/drawbicyc.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ int main(int argc, char** argv)
("qv_qc_2_6_10_min", po::value<bool>()->default_value(false) , "plot comparison of qv and qc fields at 2, 6 and 10 min?")
("dir", po::value<std::string>()->required() , "directory containing out_lgrngn")
("micro", po::value<std::string>()->required(), "one of: blk_1m, blk_2m, lgrngn")
("type", po::value<std::string>()->required(), "one of: dycoms, moist_thermal, rico")//, base_prflux_vs_clhght")
("type", po::value<std::string>()->required(), "one of: dycoms, moist_thermal, rico, Lasher_Trapp")//, base_prflux_vs_clhght")
;

po::variables_map vm;
Expand All @@ -33,8 +33,8 @@ int main(int argc, char** argv)

// handling the "type" option
std::string type = vm["type"].as<std::string>();
if(type != "dycoms" && type != "moist_thermal" && type != "rico")// && type != "base_prflux_vs_clhght")
throw std::runtime_error("Unrecognized 'type' option, only dycoms, rico, moist_thermal available now");//, base_prflux_vs_clhght available now");
if(type != "dycoms" && type != "moist_thermal" && type != "rico" && type != "Lasher_Trapp")// && type != "base_prflux_vs_clhght")
throw std::runtime_error("Unrecognized 'type' option, only dycoms, rico, Lasher_Trapp, moist_thermal available now");//, base_prflux_vs_clhght available now");

// should profiles be normalized by inversion height
const bool normalize_prof = type == "dycoms";
Expand All @@ -56,10 +56,10 @@ int main(int argc, char** argv)
H5::DataSet h5d = h5f.openDataSet("G");
H5::DataSpace h5s = h5d.getSpace();
int NDims = h5s.getSimpleExtentNdims();

// detecting if subgrid model was on
bool sgs = true;
try
try
{
auto h5g = h5f.openGroup("sgs");
}
Expand Down
Loading