Skip to content

Commit

Permalink
nc vs af plots: python3 and af profiles plot
Browse files Browse the repository at this point in the history
  • Loading branch information
@pdziekan
pdziekan committed Oct 8, 2020
1 parent 2409a68 commit ec1dd8f
Show file tree
Hide file tree
Showing 2 changed files with 43 additions and 43 deletions.
20 changes: 10 additions & 10 deletions NC_vs_AF/nc_vs_af.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@

for timestep in timesteps:
plt.clf()
print timestep
print(timestep)
filename = input_dir + "/timestep" + str(timestep).zfill(10) + ".h5"
#rl = (h5py.File(filename, "r")["cloud_rw_mom3"][:,:,:] + h5py.File(filename, "r")["rain_rw_mom3"][:,:,:]) * 4. / 3. * 3.1416 * 1e3; # kg/kg
rl = (h5py.File(filename, "r")["cloud_rw_mom3"][:,:,:]) * 4. / 3. * 3.1416 * 1e3; # kg/kg
Expand All @@ -38,20 +38,20 @@
# cloudiness mask - as in DYCOMS paper
cloudy_mask = np.where(nc > 20, 1, 0)

print 'nc>20 cloudy cells: ', np.sum(cloudy_mask)
print 'nc>20 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask)
print('nc>20 cloudy cells: ', np.sum(cloudy_mask))
print('nc>20 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask))

# cloudiness mask - rl > 1e-4
cloudy_mask = np.where(rl > 1e-4, 1, 0)

print 'rl>1e-4 cloudy cells: ', np.sum(cloudy_mask)
print 'rl>1e-4 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask)
print('rl>1e-4 cloudy cells: ', np.sum(cloudy_mask))
print('rl>1e-4 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask))

# cloudiness mask - as in RICO paper
cloudy_mask = np.where(rl > 1e-5, 1, 0)

print 'rl>1e-5 cloudy cells: ', np.sum(cloudy_mask)
print 'rl>1e-5 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask)
print('rl>1e-5 cloudy cells: ', np.sum(cloudy_mask))
print('rl>1e-5 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask))

# ---- adiabatic LWC ----
AF = np.zeros([nx, ny, nz])
Expand Down Expand Up @@ -88,7 +88,7 @@
parcel_th = np.mean(clb_th[clb_th>0])
parcel_rl = 0

print 'parcel init: th = ', parcel_th, ' rv = ', parcel_rv
print('parcel init: th = ', parcel_th, ' rv = ', parcel_rv)

for k in np.arange(nz):
parcel_T = parcel_th * lcmn.exner(p_e.astype(float)[k])
Expand All @@ -100,8 +100,8 @@
parcel_rl += delta_rv
adia_rl[k] = parcel_rl

print p_e
print adia_rl
print(p_e)
print(adia_rl)

#for i, j in zip(np.arange(nx), np.arange(ny)):
# parcel_rv = rv[i, j, clb_idx[i, j]]
Expand Down
66 changes: 33 additions & 33 deletions NC_vs_AF/nc_vs_af_AFperCol.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,8 @@
# adiabatic rl calculated separately for each column

from sys import argv, path, maxsize
path.insert(0, "/home/piotr/usr/local/lib/python2.7/dist-packages/")
#path.insert(0, "/home/piotr/usr/local/lib/python2.7/dist-packages/")
path.insert(0,"/home/piotr/usr/local/lib/python3/dist-packages/")

import h5py
import numpy as np
Expand All @@ -17,8 +18,8 @@

evap_lat = 2.5e6 # [J/kg] latent heat of evaporation

timesteps = [12000, 36000, 72000]
#timesteps = [36000]
#timesteps = [12000, 36000, 72000]
timesteps = [13200]

input_dir = argv[1]
outfile = argv[2]
Expand All @@ -30,7 +31,6 @@

for timestep in timesteps:
plt.clf()
print timestep
filename = input_dir + "/timestep" + str(timestep).zfill(10) + ".h5"
#rl = (h5py.File(filename, "r")["cloud_rw_mom3"][:,:,:] + h5py.File(filename, "r")["rain_rw_mom3"][:,:,:]) * 4. / 3. * 3.1416 * 1e3; # kg/kg
rl = (h5py.File(filename, "r")["cloud_rw_mom3"][:,:,:]) * 4. / 3. * 3.1416 * 1e3; # kg/kg
Expand All @@ -39,23 +39,24 @@
# cloudiness mask - as in DYCOMS paper
cloudy_mask = np.where(nc > 20, 1, 0)

print 'nc>20 cloudy cells: ', np.sum(cloudy_mask)
print 'nc>20 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask)
print('nc>20 cloudy cells: ', np.sum(cloudy_mask))
print('nc>20 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask))

# cloudiness mask - rl > 1e-4
cloudy_mask = np.where(rl > 1e-4, 1, 0)

print 'rl>1e-4 cloudy cells: ', np.sum(cloudy_mask)
print 'rl>1e-4 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask)
print('rl>1e-4 cloudy cells: ', np.sum(cloudy_mask))
print('rl>1e-4 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask))

# cloudiness mask - as in RICO paper
cloudy_mask = np.where(rl > 1e-5, 1, 0)
cloudy_mask_used = cloudy_mask

print 'rl>1e-5 cloudy cells: ', np.sum(cloudy_mask)
print 'rl>1e-5 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask)
print('rl>1e-5 cloudy cells: ', np.sum(cloudy_mask))
print('rl>1e-5 mean nc in cloudy cells: ', np.sum(nc * cloudy_mask) / np.sum(cloudy_mask))

hght_abv_clb = np.zeros([nx, ny, nz])
hght = np.arange(nz) * dz

# ---- adiabatic LWC ----
AF = np.zeros([nx, ny, nz])
Expand Down Expand Up @@ -105,20 +106,6 @@
parcel_rl += delta_rv
adia_rl[i, j, k] = parcel_rl

#for i, j in zip(np.arange(nx), np.arange(ny)):
# parcel_rv = rv[i, j, clb_idx[i, j]]
# parcel_th = th[i, j, clb_idx[i, j]]
# for k in np.arange(nz):
# if k < clb_idx[i, j]:
# adia_rl[i,j,k] = 0
# else:
# parcel_T = parcel_th * lcmn.exner(p_e.astype(float)[k])
# delta_rv = parcel_rv - lcmn.r_vs(parcel_T, p_e.astype(float)[k])
# parcel_rv -= delta_rv
# parcel_th += delta_rv * evap_lat / lcmn.c_pd / lcmn.exner(p_e.astype(float)[k])
# adia_rl[i,j,k] = rl[i,j,k] / delta_rv


#adia_rl = np.where(adia_rl > 0., adia_rl, 0)
#print adia_rl

Expand All @@ -131,8 +118,11 @@
for k in np.arange(nz):
# if rl[i,j,k] > 0 and adia_rl[k] == 0:
# print 'i: ',i, ' j: ',j, ' k: ',k, 'rl: ', rl[i,j,k], 'adia_rl: ', adia_rl[k], 'nc: ', nc[i,j,k]
AF[i, j, k] = rl[i,j,k] / adia_rl[i, j, k]
hght_abv_clb[i, j, k] = (k - clb_idx[i,j]) * 40
if adia_rl[i, j, k] == 0:
AF[i, j, k] = 0
else:
AF[i, j, k] = rl[i,j,k] / adia_rl[i, j, k]
hght_abv_clb[i, j, k] = (k - clb_idx[i,j]) * dz
# print 'i: ',i, ' j: ',j, ' k: ',k, 'rl: ', rl[i,j,k], 'adia_rl: ', adia_rl[k], 'nc: ', nc[i,j,k], 'AF: ', AF[i,j,k]

#print cloudy_mask_used[cloudy_mask_used>0]
Expand All @@ -152,9 +142,9 @@
# jet, hot
cb = plt.colorbar()
cb.set_label("Height above cloud base [m]")
plt.clim(0,1400)
plt.xlim(0,10)
plt.ylim(0,200)
# plt.clim(0,1400)
# plt.xlim(0,10)
# plt.ylim(0,200)

plt.xlabel('AF')
plt.ylabel('Nc [1/cc]')
Expand All @@ -173,7 +163,7 @@
# jet, hot
cb = plt.colorbar()
cb.set_label("Height above cloud base [m]")
plt.clim(0,1400)
# plt.clim(0,1400)
plt.gca().set_xlim(left=0)
plt.gca().set_ylim(bottom=0)
plt.ylabel('r_l')
Expand All @@ -188,9 +178,9 @@
# jet, hot
cb = plt.colorbar()
cb.set_label("Height above cloud base [m]")
plt.clim(0,1400)
plt.xlim(0,5e-3)
plt.ylim(0,200)
# plt.clim(0,1400)
# plt.xlim(0,5e-3)
# plt.ylim(0,200)
plt.xlabel('r_l')
plt.ylabel('Nc [1/cc]')
plt.savefig(outfile + '_NCvsrl_AFperCol_' + str(timestep) +'.png')
Expand All @@ -200,3 +190,13 @@
# plt.scatter(hght_abv_clb[cloudy_mask_used==1].flatten(), nc[cloudy_mask_used==1].flatten(), c = AF[cloudy_mask_used==1].flatten())
# plt.colorbar()
# plt.savefig(outfile + '_NCvsHght_AFperCol_' + str(timestep) +'.png')

# plot AF profile
plt.clf()
AF[AF==0] = np.nan
AF[cloudy_mask_used==0] = np.nan
AF_mean = np.nanmean(AF, axis=(0,1))
AF_mean[np.isnan(AF_mean)] = 0 # set AF=0 above and below cloud
# plt.xlim(0,1.5)
plt.plot(AF_mean, hght)
plt.savefig(outfile + '_AFprofile_AFperCol_' + str(timestep) +'.png')

0 comments on commit ec1dd8f

Please sign in to comment.