MeteoInfoLab脚本示例：计算K指数

MeteoInfo

利用FNL数据计算K指数的例子，计算方法参考了此贴：http://bbs.06climate.com/forum.php?mod=viewthread&tid=1854 。

1. #Add file
   
2. f = addfile('D:/Temp/grib/fnl_20110416_00_00')
   
3. #Calculate K index
   
4. #850hPa
   
5. T850 = f['Temperature_isobaric'][0,'85000.0','10:35','90:125'] - 273.16
   
6. rh = f['Relative_humidity_isobaric'][0,'85000.0','10:35','90:125']
   
7. Td850 = T850-((14.55+0.114*T850)*(1-0.01*rh) + pow((2.5+0.007*T850)*(1-0.01*rh),3) + (15.9+0.117*T850)*pow((1-0.01*rh),14))
   
8. #700hPa
   
9. T700 = f['Temperature_isobaric'][0,'70000.0','10:35','90:125'] - 273.16
   
10. rh = f['Relative_humidity_isobaric'][0,'70000.0','10:35','90:125']
    
11. Td700 = T700-((14.55+0.114*T700)*(1-0.01*rh) + pow((2.5+0.007*T700)*(1-0.01*rh),3) + (15.9+0.117*T700)*pow((1-0.01*rh),14))
    
12. #500hPa
    
13. T500 = f['Temperature_isobaric'][0,'50000.0','10:35','90:125'] - 273.16
    
14. K = T850-T500+Td850-(T700-Td700)
    
15. #Plot
    
16. axesm()
    
17. lworld = shaperead('D:/Temp/Map/country1.shp')
    
18. geoshow(lworld, edgecolor='k')
    
19. levs = arange(-40,36,2.5)
    
20. layer = contourm(K, levs)
    
21. clabel(layer)
    
22. colorbar(layer)
    
23. t = f.gettime(0)
    
24. title('K index (' + t.strftime('%Y-%m-%d %H:00') + ')')

#Add file<br /> f = addfile('D:/Temp/grib/fnl_20110416_00_00')<br /> #Calculate K index<br /> #850hPa<br /> T850 = f['Temperature_isobaric'][0,'85000.0','10:35','90:125'] - 273.16<br /> rh = f['Relative_humidity_isobaric'][0,'85000.0','10:35','90:125']<br /> Td850 = T850-((14.55+0.114*T850)*(1-0.01*rh) + pow((2.5+0.007*T850)*(1-0.01*rh),3) + (15.9+0.117*T850)*pow((1-0.01*rh),14))<br /> #700hPa<br /> T700 = f['Temperature_isobaric'][0,'70000.0','10:35','90:125'] - 273.16<br /> rh = f['Relative_humidity_isobaric'][0,'70000.0','10:35','90:125']<br /> Td700 = T700-((14.55+0.114*T700)*(1-0.01*rh) + pow((2.5+0.007*T700)*(1-0.01*rh),3) + (15.9+0.117*T700)*pow((1-0.01*rh),14))<br /> #500hPa<br /> T500 = f['Temperature_isobaric'][0,'50000.0','10:35','90:125'] - 273.16<br /> K = T850-T500+Td850-(T700-Td700)<br /> #Plot<br /> axesm()<br /> lworld = shaperead('D:/Temp/Map/country1.shp')<br /> geoshow(lworld, edgecolor='k')<br /> levs = arange(-40,36,2.5)<br /> layer = contourm(K, levs)<br /> clabel(layer)<br /> colorbar(layer)<br /> t = f.gettime(0)<br /> title('K index (' + t.strftime('%Y-%m-%d %H:00') + ')')

好多鱼

很好的帖子，留下来了。

398891490

感谢分享