How can we calculate potential vorticity (PV) along the trajectory by using MeteoInfoLab?

I want to know the PV values to examine the impact of stratosphere.

I calculated back trajectories by using GDAS (1 degree grids) data.

I tried to write a script based on this example (http://meteothink.org/examples/meteoinf ... _data.html), but I haven't got PV values.

Can someone please share the script for calculation of PV along the trajectory?

## How to calculate potential vorticity (PV) along the trajectory

### Re: How to calculate potential vorticity (PV) along the trajectory

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 (http://meteothink.org/examples/meteoinf ... icity.html). 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...'
```

### Re: How to calculate potential vorticity (PV) along the trajectory

Thank you for answering my question.

The calculation went really well.

The calculation went really well.