Moisture potential vorticity

The example to calcluate moisture potential vorticity.

# Calculate moisture potential vorticity
# Set working directory
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'
tidx = 0
rh = meteof['RELH'][tidx,:,latlim,lonlim]
nz,ny,nx = rh.shape
lon = rh.dimvalue(2)
lat = rh.dimvalue(1)
lev = rh.dimvalue(0)
t0 = meteof['TEMP'][tidx,:nz,latlim,lonlim]
uwnd = meteof['UWND'][tidx,:nz,latlim,lonlim]
vwnd = meteof['VWND'][tidx,:nz,latlim,lonlim]

dx, dy = meteolib.lat_lon_grid_deltas(lon, lat)
vort = meteolib.vorticity(uwnd, vwnd, dx, dy)
prs = zeros([nz,ny,nx])
prs = dim_array(prs, rh.dims)
for i in range(nz):
    prs[i] = lev[i]

# Calculate pseudo-equivalent potential temperature
print('Clalulate pseudo-equivalent potential temperature...')
dewp = meteolib.dewpoint_from_relative_humidity(t0, rh / 100)
eqt = meteolib.equivalent_potential_temperature(prs, t0, dewp)

# Calculate moisture potential vorticity
print('Calculate moisture potential vorticity...')
mpv_3d = zeros([nz,ny,nx], dtype='double')
mpv_3d = dim_array(mpv_3d, rh.dims)
mpv_3d.setdimvalue(0, lev[1:nz-1])

tt = meteof.gettime(tidx)
print(tt.strftime('%Y-%m-%d %H:00'))
f1 = 2*7.292*sin(lat*3.14159/180.0)*0.00001
f = f1.reshape(ny, 1).repeat(nx, axis=1)
g = 9.8
dx1 = 6370949.0*cos(lat*3.14159/180.0)*3.14159/180.0
dx = dx1.reshape(ny, 1).repeat(nx, axis=1)
dy = 6370949.0*3.14159/180.0
for z in range(1, nz-1):
    dp = 100*(lev[z-1]-lev[z+1])
    deqt = eqt[z-1]-eqt[z+1]
    du = uwnd[z-1]-uwnd[z+1]
    dv = vwnd[z-1]-vwnd[z+1]
    dtx = meteolib.cdiff(eqt[z], 1)
    dty = meteolib.cdiff(eqt[z], 0)
    mpv1 = -g*(vort[z]+f)*deqt/dp
    mpv2 = g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))
    mpv = mpv1+mpv2
    mpv_3d[z-1] = mpv

#Plot test
axesm()
geoshow('coastline', color='k')
z = 5
clevs = arange(-2, 2.1, 0.5)
layer = contourf(mpv_3d[z]*1e6, clevs, cmap='matlab_jet')
colorbar(layer)
title('Moisture potential vorticity (%i hPa)\n' % lev[z] +
    tt.strftime('%Y-%m-%d %H:00'))