Skip to content

Commit 6f41434

Browse files
quantheoryquantheory
committed
Random changes and new scripts.
1 parent b9f29f5 commit 6f41434

9 files changed

Lines changed: 3364 additions & 167 deletions

‎cluster_twmt_plot.py‎

Lines changed: 550 additions & 0 deletions
Large diffs are not rendered by default.

‎lwp_prect.py‎

Lines changed: 254 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,254 @@
1+
#!/usr/bin/env python
2+
3+
import numpy as np
4+
import scipy.linalg as la
5+
import scipy.stats as stats
6+
import matplotlib
7+
matplotlib.use('Agg')
8+
import matplotlib.pyplot as plt
9+
import netCDF4 as nc4
10+
import numdifftools as ndt
11+
import sklearn.cluster as clstr
12+
13+
from mg2_constants import *
14+
15+
HIST_FILE_NAME = "/p/lscratchh/santos36/ACME/short_timestep_ctrl_diags/run/short_timestep_ctrl_diags.cam.h1.0001-01-01-72000.nc"
16+
SHORT_HIST_FILE_NAME = "/p/lscratchh/santos36/ACME/short_timestep_diags/run/short_timestep_diags.cam.h1.0001-01-01-72000.nc"
17+
18+
nbins = 20
19+
qsmall = 1.e-18
20+
21+
cmap = plt.get_cmap('coolwarm')
22+
23+
file = nc4.Dataset(HIST_FILE_NAME, 'r')
24+
25+
ncol = len(file.dimensions['ncol'])
26+
lev = len(file.dimensions['lev'])
27+
ilev = len(file.dimensions['ilev'])
28+
29+
lwp = file.variables["MPLWPI"][0,:]
30+
twp = file.variables["MPLWPI"][0,:] + file.variables["MPIWPI"][0,:]
31+
precl = file.variables["PRECL"][0,:]
32+
prect = file.variables["PRECT"][0,:]
33+
qrsedten = file.variables["QRSEDTEN"][0,:,:]
34+
evaprain = file.variables["EVAPPREC"][0,:,:] - file.variables["EVAPSNOW"][0,:,:]
35+
prc = file.variables["PRCO"][0,:,:]
36+
pra = file.variables["PRAO"][0,:,:]
37+
pmids = file.variables['lev'][:] * 100.
38+
pints = file.variables['ilev'][:] * 100.
39+
40+
precl_long = precl
41+
prect_long = prect
42+
43+
PLT_FILE = "lwp_precl"
44+
45+
plt.plot(lwp, precl, '.')
46+
plt.xlabel("Water path (kg/m^2)")
47+
plt.ylabel("PRECL (m/s)")
48+
plt.savefig(PLT_FILE+".eps")
49+
plt.savefig(PLT_FILE+".png")
50+
plt.close()
51+
52+
plt.plot(lwp, precl, '.')
53+
plt.xlabel("Water path (kg/m^2)")
54+
plt.ylabel("PRECL (m/s)")
55+
plt.axis([0., 6., 0., 8.e-7])
56+
plt.savefig(PLT_FILE+"_scaled.eps")
57+
plt.savefig(PLT_FILE+"_scaled.png")
58+
plt.close()
59+
60+
PLT_FILE = "lwp_precl_hist"
61+
62+
plt.hist2d(lwp, precl, bins=nbins, range=((0.05, 1), (7.e-9, 1.e-7)), cmap=cmap)
63+
plt.xlabel("Water path (kg/m^2)")
64+
plt.ylabel("PRECL (m/s)")
65+
plt.axis('tight')
66+
plt.savefig(PLT_FILE+".eps")
67+
plt.savefig(PLT_FILE+".png")
68+
plt.close()
69+
70+
qrsedtot = np.zeros((ncol,))
71+
evapraintot = np.zeros((ncol,))
72+
for i in range(lev):
73+
qrsedtot[:] = qrsedtot[:] + np.abs(qrsedten[i,:])*(pints[i] - pints[i+1])/(lev*gravit)
74+
evapraintot[:] = evapraintot[:] + evaprain[i,:]*(pints[i] - pints[i+1])/gravit
75+
76+
evaprainavg = np.zeros((lev,))
77+
prcavg = np.zeros((lev,))
78+
praavg = np.zeros((lev,))
79+
80+
for n in range(ncol):
81+
evaprainavg[:] = evaprainavg[:] + evaprain[:,n]*(pints[:-1] - pints[1:])/gravit
82+
prcavg[:] = prcavg[:] + prc[:,n]*(pints[:-1] - pints[1:])/gravit
83+
praavg[:] = praavg[:] + pra[:,n]*(pints[:-1] - pints[1:])/gravit
84+
85+
PLT_FILE = "evaprain_qrsedtot"
86+
87+
plt.plot(evapraintot, qrsedtot, '.')
88+
plt.xlabel("Rain evaporation rate (kg/m^2)")
89+
plt.ylabel("Sedimentation column average rate (kg/m^2)")
90+
plt.savefig(PLT_FILE+".eps")
91+
plt.savefig(PLT_FILE+".png")
92+
plt.close()
93+
94+
PLT_FILE = "precl_hist"
95+
96+
bins = np.linspace(qsmall, 2.e-8, nbins+1)
97+
plt.hist(precl, bins=bins)
98+
plt.xlabel("PRECL (m/s)")
99+
plt.axis([0., 2.e-8, 0, 30000])
100+
plt.savefig(PLT_FILE+".eps")
101+
plt.savefig(PLT_FILE+".png")
102+
plt.close()
103+
104+
PLT_FILE = "evaprain_hist"
105+
106+
bins = np.linspace(-5.e-5, -qsmall, nbins+1)
107+
plt.hist(evapraintot, bins=bins)
108+
plt.xlabel("Rain evaporation rate (kg/m^2)")
109+
plt.savefig(PLT_FILE+".eps")
110+
plt.savefig(PLT_FILE+".png")
111+
plt.close()
112+
113+
file = nc4.Dataset(SHORT_HIST_FILE_NAME, 'r')
114+
115+
ncol = len(file.dimensions['ncol'])
116+
lev = len(file.dimensions['lev'])
117+
ilev = len(file.dimensions['ilev'])
118+
119+
lwp = file.variables["MPLWPI"][0,:]
120+
twp = file.variables["MPLWPI"][0,:] + file.variables["MPIWPI"][0,:]
121+
precl = file.variables["PRECL"][0,:]
122+
prect = file.variables["PRECT"][0,:]
123+
qrsedten = file.variables["QRSEDTEN"][0,:,:]
124+
evaprain = file.variables["EVAPPREC"][0,:,:] - file.variables["EVAPSNOW"][0,:,:]
125+
prc = file.variables["PRCO"][0,:,:]
126+
pra = file.variables["PRAO"][0,:,:]
127+
pmids = file.variables['lev'][:] * 100.
128+
pints = file.variables['ilev'][:] * 100.
129+
130+
PLT_FILE = "lwp_precl_short"
131+
132+
plt.plot(lwp, precl, '.')
133+
plt.xlabel("Water path (kg/m^2)")
134+
plt.ylabel("PRECT (m/s)")
135+
plt.savefig(PLT_FILE+".eps")
136+
plt.savefig(PLT_FILE+".png")
137+
plt.close()
138+
139+
plt.plot(lwp, precl, '.')
140+
plt.xlabel("Water path (kg/m^2)")
141+
plt.ylabel("PRECT (m/s)")
142+
plt.axis([0., 6., 0., 8.e-7])
143+
plt.savefig(PLT_FILE+"_scaled.eps")
144+
plt.savefig(PLT_FILE+"_scaled.png")
145+
plt.close()
146+
147+
PLT_FILE = "lwp_precl_hist_short"
148+
149+
plt.hist2d(lwp, precl, bins=nbins, range=((0.05, 1.), (5.e-9, 1.e-7)), cmap=cmap)
150+
plt.xlabel("Water path (kg/m^2)")
151+
plt.ylabel("PRECL (m/s)")
152+
plt.axis('tight')
153+
plt.savefig(PLT_FILE+".eps")
154+
plt.savefig(PLT_FILE+".png")
155+
plt.close()
156+
157+
qrsedtot = np.zeros((ncol,))
158+
evapraintot = np.zeros((ncol,))
159+
for i in range(lev):
160+
qrsedtot[:] = qrsedtot[:] + np.abs(qrsedten[i,:])*(pints[i] - pints[i+1])/(lev*gravit)
161+
evapraintot[:] = evapraintot[:] + evaprain[i,:]*(pints[i] - pints[i+1])/gravit
162+
163+
evaprainavg_short = np.zeros((lev,))
164+
prcavg_short = np.zeros((lev,))
165+
praavg_short = np.zeros((lev,))
166+
167+
for n in range(ncol):
168+
evaprainavg_short[:] = evaprainavg_short[:] + evaprain[:,n]*(pints[:-1] - pints[1:])/gravit
169+
prcavg_short[:] = prcavg_short[:] + prc[:,n]*(pints[:-1] - pints[1:])/gravit
170+
praavg_short[:] = praavg_short[:] + pra[:,n]*(pints[:-1] - pints[1:])/gravit
171+
172+
PLT_FILE = "evaprain_qrsedtot_short"
173+
174+
plt.plot(evapraintot, qrsedtot, '.')
175+
plt.xlabel("Rain evaporation rate (kg/m^2)")
176+
plt.ylabel("Sedimentation column average rate (kg/m^2)")
177+
plt.savefig(PLT_FILE+".eps")
178+
plt.savefig(PLT_FILE+".png")
179+
plt.close()
180+
181+
PLT_FILE = "evaprain_vertical"
182+
183+
plt.plot(evaprainavg, pmids, label='300s timestep')
184+
plt.plot(evaprainavg_short, pmids, label='1s timestep')
185+
plt.xlabel("Rain evaporation rate (kg/m^2)")
186+
plt.ylabel("Pressure (Pa)")
187+
plt.axis('tight')
188+
plt.gca().invert_yaxis()
189+
plt.legend(loc='best')
190+
plt.savefig(PLT_FILE+".eps")
191+
plt.savefig(PLT_FILE+".png")
192+
plt.close()
193+
194+
PLT_FILE = "prc_vertical"
195+
196+
plt.plot(prcavg, pmids, label='300s timestep')
197+
plt.plot(prcavg_short, pmids, label='1s timestep')
198+
plt.xlabel("Autoconversion rate (kg/m^2)")
199+
plt.ylabel("Pressure (Pa)")
200+
plt.axis('tight')
201+
plt.gca().invert_yaxis()
202+
plt.legend(loc='best')
203+
plt.savefig(PLT_FILE+".eps")
204+
plt.savefig(PLT_FILE+".png")
205+
plt.close()
206+
207+
PLT_FILE = "pra_vertical"
208+
209+
plt.plot(praavg, pmids, label='300s timestep')
210+
plt.plot(praavg_short, pmids, label='1s timestep')
211+
plt.xlabel("Accretion rate (kg/m^2)")
212+
plt.ylabel("Pressure (Pa)")
213+
plt.axis('tight')
214+
plt.gca().invert_yaxis()
215+
plt.legend(loc='best')
216+
plt.savefig(PLT_FILE+".eps")
217+
plt.savefig(PLT_FILE+".png")
218+
plt.close()
219+
220+
PLT_FILE = "precl_hist_short"
221+
222+
bins = np.linspace(qsmall, 2.e-8, nbins+1)
223+
plt.hist(precl, bins=bins)
224+
plt.xlabel("PRECL (m/s)")
225+
plt.axis([0., 2.e-8, 0, 30000])
226+
plt.savefig(PLT_FILE+".eps")
227+
plt.savefig(PLT_FILE+".png")
228+
plt.close()
229+
230+
PLT_FILE = "precl_hist_overlap"
231+
232+
bins = np.linspace(qsmall, 2.e-8, nbins+1)
233+
plt.hist(precl, bins=bins, label='1s timestep')
234+
plt.hist(precl_long, bins=bins, label='300s timestep')
235+
plt.xlabel("PRECL (m/s)")
236+
plt.axis([0., 2.e-8, 0, 30000])
237+
plt.legend(loc='best')
238+
plt.savefig(PLT_FILE+".eps")
239+
plt.savefig(PLT_FILE+".png")
240+
plt.close()
241+
242+
PLT_FILE = "evaprain_hist_short"
243+
244+
bins = np.linspace(-1.5e-5, -qsmall, nbins+1)
245+
plt.hist(evapraintot, bins=bins)
246+
plt.xlabel("Rain evaporation rate (kg/m^2)")
247+
plt.savefig(PLT_FILE+".eps")
248+
plt.savefig(PLT_FILE+".png")
249+
plt.close()
250+
251+
print('Mean PRECL @ 300s = ', np.mean(precl_long))
252+
print('Mean PRECL @ 1s = ', np.mean(precl))
253+
print('Mean PRECT @ 300s = ', np.mean(prect_long))
254+
print('Mean PRECT @ 1s = ', np.mean(prect))

‎mg2_constants.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -30,5 +30,5 @@
3030
prc_exp1 = -1.2
3131
cld_sed = 1.
3232
mg_prc_coeff_fix = True
33-
alpha_grad = 1.
33+
alpha_grad = 2.
3434
beta_grad = 1.

0 commit comments

Comments
 (0)