MeteoInfoLab脚本示例：计算假相当位温

MeteoInfo

计算假相当位温的例子，算法参考了此贴：http://bbs.06climate.com/forum.php?mod=viewthread&tid=22546

1. f_air = addfile('D:/Temp/nc/air.2011.nc')
   
2. f_rhum = addfile('D:/Temp/nc/rhum.2011.nc')
   
3. levels = f_air['level'][:]
   
4. z = 0
   
5. t0 = f_air['air'][0,z,:,:]
   
6. rh = f_rhum['rhum'][0,z,:,:]
   
7. prs = levels[z]
   
8. es = 6.1078*exp(17.2693882*(t0-273.16)/(t0-35.86))
   
9. qq = rh*(0.62197*es/(prs-0.378*es))/100.
   
10. e = prs*qq/(0.62197+qq)+1e-10
    
11. tlcl = 55.0+2840.0/(3.5*log(t0)-log(e)-4.805)
    
12. theta = t0*pow((1000/prs),(0.2854*(1.0-0.28*qq)))
    
13. eqt = theta*exp(((3376./tlcl)-2.54)*qq*(1.0+0.81*qq))
    
14. thse = eqt-273.15
    
15. #Plot
    
16. axesm()
    
17. lworld = shaperead('D:/Temp/Map/country1.shp')
    
18. geoshow(lworld, edgecolor='k')
    
19. layer = contourm(eqt, 12)
    
20. clabel(layer)
    
21. colorbar(layer)
    
22. title('pseudo-equivalent potential temperature')

f_air = addfile('D:/Temp/nc/air.2011.nc')<br /> f_rhum = addfile('D:/Temp/nc/rhum.2011.nc')<br /> levels = f_air['level'][:]<br /> z = 0<br /> t0 = f_air['air'][0,z,:,:]<br /> rh = f_rhum['rhum'][0,z,:,:]<br /> prs = levels[z]<br /> es = 6.1078*exp(17.2693882*(t0-273.16)/(t0-35.86))<br /> qq = rh*(0.62197*es/(prs-0.378*es))/100.<br /> e = prs*qq/(0.62197+qq)+1e-10<br /> tlcl = 55.0+2840.0/(3.5*log(t0)-log(e)-4.805)<br /> theta = t0*pow((1000/prs),(0.2854*(1.0-0.28*qq)))<br /> eqt = theta*exp(((3376./tlcl)-2.54)*qq*(1.0+0.81*qq))<br /> thse = eqt-273.15<br /> #Plot<br /> axesm()<br /> lworld = shaperead('D:/Temp/Map/country1.shp')<br /> geoshow(lworld, edgecolor='k')<br /> layer = contourm(eqt, 12)<br /> clabel(layer)<br /> colorbar(layer)<br /> title('pseudo-equivalent potential temperature')

kmdqlb

王老师威武

baggio

{:eb302:}{:eb302:}{:eb302:}{:eb302:}

有事秦天

王老师您好，我想绘制价相当位温剖面图，怎样才能每层计算的时候都用到当前层次高度，在叠加，这个问题该怎么解决？
f = addfile('D:/M07/a201507.nc')
levels = f['level'][:]
time = 106
level = '1000:200'
lat = '36:44'
lon = '112:118'
t0 = f['t'][time,level,lat,lon][::-1,:,:]
rh = f['r'][time,level,lat,lon][::-1,:,:]
prs = levels[level]
es = 6.1078*exp(17.2693882*(t0-273.16)/(t0-35.86))
qq = rh*(0.62197*es/(prs-0.378*es))/100.
e = prs*qq/(0.62197+qq)+1e-10
tlcl = 55.0+2840.0/(3.5*log(t0)-log(e)-4.805)
theta = t0*pow((1000/prs),(0.2854*(1.0-0.28*qq)))
eqt = theta*exp(((3376./tlcl)-2.54)*qq*(1.0+0.81*qq))
thse = eqt-273.15

wd1 = eqt[:,'40',:]
lev1 = wd1.dimvalue(0)
lev1 = lev1[::-1]
lev2 = meteo.p2h(lev1)
levels = []
for j in range(0, len(lev1)):
    levels.append('%i' % lev1[j])
wd1.setdimvalue(0, lev2)

axesm()
lworld = shaperead('D:/MeteoInfo/map/bou2_4p.shp')
geoshow(lworld, edgecolor='k')
layer = contourm(eqt)
clabel(layer)
#colorbar(layer)
title('pseudo-equivalent potential temperature')
xlim(115, 119)
ylim(36, 44)

瑛雩

不知道为什么，一模一样的数据和资料，我计算700hPa上的假相当位温，MeteoInfoLab计算得到287.84，grads计算得到318.65搞不懂什么原因啊

瑛雩

已经搞清楚了，气压不能是整数，后面要加个点变成小数折腾了好久

Eegle

瑛雩 发表于 2019-9-15 19:43
已经搞清楚了，气压不能是整数，后面要加个点变成小数折腾了好久

结果是一样的吗？我算出的850hPa大概是312K左右，可是同一时刻的初始场应该是350K左右，不懂错在哪了

Eegle

王老师您好，我用您的算法计算prs=850.0，T=293.16（20℃），rh=85%带入计算，eqt才308K。好像不太准确呀，至少也要330K左右，不知道是哪里错了

sunflower

感谢楼主分享

猫猫2003

瑛雩 发表于 2019-9-15 19:18
不知道为什么，一模一样的数据和资料，我计算700hPa上的假相当位温，MeteoInfoLab计算得到287.84，grads计 ...

我也遇到同样的问题，你找到原因了吗？