You can calculate PV array from the existing variables and than interpolate to the end point positions of the trajectories. There is a MPV example online (
). Replace pseudo-equivalent potential temperature to potential temperature will get PV result.
Code: Select all
# Calculate potential vorticity
# Set working directory
trajDir = 'D:/Temp/HYSPLIT'
meteoDir = 'D:/Temp/ARL'
# Open meteorological data file
print 'Open meteorological data file...'
meteofn = os.path.join(meteoDir, 'gdas1.mar15.w5')
print 'Meteorological file: ' + meteofn
meteof = addfile(meteofn)
# Read data
print 'Read data...'
latlim = '10:60'
lonlim = '60:140'
rh = meteof['RELH'][:,:,latlim,lonlim]
nt,nz,ny,nx = rh.shape
lat = rh.dimvalue(2)
lev = rh.dimvalue(1)
t0 = meteof['TEMP'][:,:nz,latlim,lonlim]
uwnd = meteof['UWND'][:,:nz,latlim,lonlim]
vwnd = meteof['VWND'][:,:nz,latlim,lonlim]
vort = hcurl(uwnd, vwnd)
prs = zeros([nt,nz,ny,nx])
prs = dim_array(prs, rh.dims)
for i in range(nz):
prs[:,i,:,:] = lev[i]
# Calculate potential temperature
print 'Clalulate potential temperature...'
pt = meteo.potential_temperature(prs, t0)
# Calculate potential vorticity
print 'Calculate potential vorticity...'
pv = zeros([nt,nz,ny,nx], dtype='double')
pv = dim_array(pv, rh.dims)
pv.setdimvalue(1, lev[1:nz-1])
for t in range(nt):
tt = meteof.gettime(t)
print tt.strftime('%Y-%m-%d %H:00')
for z in range(1, nz-1):
#print '\tLevel: %i' % z
f = zeros([ny,nx])
f1 = 2*7.292*sin(lat*3.14159/180.0)*0.00001
for i in range(nx):
f[:,i] = f1
g = 9.8
dp = 100*(lev[z-1]-lev[z+1])
dpt = pt[t,z-1,:,:]-pt[t,z+1,:,:]
du = uwnd[t,z-1,:,:]-uwnd[t,z+1,:,:]
dv = vwnd[t,z-1,:,:]-vwnd[t,z+1,:,:]
dx1 = 6370949.0*cos(lat*3.14159/180.0)*3.14159/180.0
dx = zeros([ny,nx])
for i in range(nx):
dx[:,i] = dx1
dy = 6370949.0*3.14159/180.0
dtx = cdiff(pt[t,z,:,:], 1)
dty = cdiff(pt[t,z,:,:], 0)
pv1 = -g*(vort[t,z,:,:]+f)*dpt/dp
pv2 = g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))
ppv = pv1+pv2
pv[t,z-1,:,:] = ppv
#Plot test
axesm()
geoshow('country', edgecolor='k')
t = 0
tt = meteof.gettime(t)
z = 5
clevs = arange(-2,2.1,0.2)
layer = contourfm(pv[t,z,:,:]*1e6, clevs)
colorbar(layer)
title('Potential vorticity (%i hPa)\n' % lev[z] + \
tt.strftime('%Y-%m-%d %H:00'))
print 'Finish...'