MeteoInfoLab脚本示例：计算水平螺旋度

尝试了用MeteoInfoLab编写计算水平螺旋度的脚本，结果未经验证。

脚本程序：

print 'Open data files...'

f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc')

f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc')

print 'Calculate average wind field from 850 to 600 hpa...'

tidx = 173    # Jun 23, 2011

t = f_uwnd.gettime(tidx)

level = [850,600]    # 850 - 600 hPa

lat = [15,55]

lon = [70,135]

uvar = f_uwnd['uwnd']

vvar = f_vwnd['vwnd']

uwnd = uvar[tidx,level,lat,lon]

vwnd = vvar[tidx,level,lat,lon]

uc = uwnd.ave()

vc = vwnd.ave()

sp = magnitude(uwnd, vwnd)

speed = sp.ave()

#speed = magnitude(uc, vc)

direc = atan2(vc, uc)

cdirec = direc-40./180.*3.14159

if cdirec < -3.14159:

    cdirec=cdirec+3.14159*2

cuc = speed*cos(cdirec)

cvc = speed*sin(cdirec)

print 'Calculate horizontal helicity...'

Hrs = 0.

zHrs = 0.

hh = 6

while hh <= 10:

    zHrs = (uvar[tidx,hh+1,lat,lon]-cuc)*(vvar[tidx,hh,lat,lon]-cvc)- \

        (uvar[tidx,hh,lat,lon]-cuc)*(vvar[tidx,hh+1,lat,lon]-cvc)

    Hrs = Hrs+zHrs

    hh += 1

Hrs = Hrs[::-1,:]

print 'Plot...'

axesm()

mlayer = shaperead('D:/Temp/map/country1.shp')

geoshow(mlayer, edgecolor='black')

layer = contourfm(Hrs, 20)

title('Horizontal helicity (' + t.strftime('%Y-%m-%d') + ')')

colorbar(layer)