Meteoinfolab Script Example: calculating the horizontal helix degree

An attempt was made to write a script that calculates the horizontal helix in Meteoinfolab, with the result being unverified.

Script Program:

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 hpa ...'Tidx= 173#June,t =f_uwnd.gettime (TIDX) level= [850,600]#850-600 HPalat = [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.14159ifCdirec <-3.14159: Cdirec=cdirec+3.14159*2CUC= speed*cos (cdirec) CVC= speed*sin (cdirec)Print 'Calculate Horizontal helicity ...'Hrs=0.zHrs=0.hh= 6 whileHH <= 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+ = 1Hrs= 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)

Meteoinfolab Script Example: calculating the horizontal helix degree