How to calculate potential vorticity (PV) along the trajectory

Post Reply
sachio3
Posts: 2
Joined: October 31st, 2018, 3:32 am
Registered HYSPLIT User: Yes

How to calculate potential vorticity (PV) along the trajectory

Post by sachio3 » October 31st, 2018, 4:31 am

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?

yaqiang
Posts: 40
Joined: October 5th, 2015, 10:20 pm
Registered HYSPLIT User: Yes

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

Post by yaqiang » November 1st, 2018, 4:06 am

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...'

sachio3
Posts: 2
Joined: October 31st, 2018, 3:32 am
Registered HYSPLIT User: Yes

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

Post by sachio3 » November 7th, 2018, 12:59 am

Thank you for answering my question.
The calculation went really well.

Post Reply