求Si沙氏指数的计算公式

cassienju

论坛里有两个链接都有回答啊，自己搜一下

sunflower

谢谢楼主分享，学习ing

sunflower

mofangbao 发表于 2012-3-28 13:37
善用搜索，提问前先搜索 我相信深深也是搜索了告诉你的，所以...还是自己先动手比较好

我也相信你说的

sunflower

言深深 发表于 2012-3-28 10:15
http://bbs.06climate.com/forum.php?mod=viewthread&tid=1753这个不知道和你的温和不吻合，论坛以前上传的 ...

您真是个热心的人，好人有好报哒

lucifer

还是要去攒贡献呀

一线天

你好，我看到你在气象家园发帖子问沙氏指数的grads脚本，现在有没有正确版的？我在你的基础上也一直改不出来？如果你有正确的gs文件可以发我一份么？谢谢

2005211076axu

家园里边的大咖还真是多，热心人也多，多学习来提升下自己吧

南南偶卖糕Oo

请问楼主是否计算出沙氏指数了呢？我的程序在循环部分有误，能否帮我看看呢？
我的gs如下：

'reinit'
'sdfopen G:\shuju\ECdata\T_rhum\201307_T&rhum.nc'
'sdfopen G:\shuju\ECdata\data\201307_geo.nc'

'set x 69'
'set y 68'
'set t 1 120'

'set lev 850'
'define p1=z.2*9.8' *位势乘以9.8得到位势高度
'define tk1=T.1'
'define t1=tk1-273.16'
'define rh1=r.1'
'define td850=t1-((14.55+0.114*t1)*(1-0.01*rh1) + pow((2.5+0.007*t1)*(1-0.01*rh1),3)+ (15.9+0.117*t1)*pow((1-0.01*rh1),14))'
*td850是850hPa的露点温度
'define e1=6.112*pow(10,7.5*td850/(td850+237.3))' *计算水气压
'define Tlcl11=(2840.0/(3.5*log(tk1)-log(e1)-4.805))+55.0'  *Bolton公式计算的抬升凝结高度处的温度
'define r1=622*(e1/(p1-e1))' *混合比
'define theta1=tk1*pow((1000/p1),(0.2854*(1-0.28*r1)))'
'define a1=theta1*exp((3376/tlcl11-2.54)*r1*(1+0.81*r1))' *Bolton公式计算的850hPa假相当位温


'set lev 500'
'define tk2=T.1'
'define p2=z.2*9.8' *位势乘以9.8得到位势高度
'define t2=T.1-273.16'
'define rh2=r.1'
'define td500=t2-((14.55+0.114*t2)*(1-0.01*rh2) + pow((2.5+0.007*t2)*(1-0.01*rh2),3)+ (15.9+0.117*t2)*pow((1-0.01*rh2),14))'
*td500是500hPa的露点温度
'define e2=6.112*pow(10,7.5*td500/(td500+237.3))' *计算水气压
'define Tlcl22=(2840.0/(3.5*log(tk2)-log(e2)-4.805))+55.0' *Bolton公式计算的抬升凝结高度处的温度
'define r2=622*(e2/(p2-e2))' *混合比
'define t3=tlcl22-273.16'
n=0
while(n<7)
n=n+1
'define t3=t3+0.0000001* 'n'' *假设一个气温初猜值和一个增量
'tk3=t3+273.15'
'define td3=t3-((14.55+0.114*t2)*(1-0.01*rh2) + pow((2.5+0.007*t2)*(1-0.01*rh2),3)+ (15.9+0.117*t2)*pow((1-0.01*rh2),14))'
'define e3=6.112*pow(10,7.5*td3/(td3+237.3))'
'define Tlcl33=(2840.0/(3.5*log(tk3)-log(e3)-4.805))+55.0'*Bolton公式计算的抬升凝结高度处的温度
'define r3=622*(e3/(p2-e3))'混合比
'define theta2=tk3*pow((1000/p2),(0.2854*(1-0.28*r3)))'
'define a2=theta2*exp((3376/tlcl33-2.54)*r3*(1+0.81*r3))' *Bolton公式计算的850hPa假相当位温
'define a=a2-a1'

'define JDZ=abs(a)'
if(JDZ<0.01) *满足850hPa和500hPa假相当位温之差满足绝对值小于0.01
   break;endif;
endwhile
'define si=t2-t3' *沙氏指数计算公式
'd si'

'grads2ascii si G:\shuju\ECdata\SI\txt\si_2013_daily.txt %g 1 1'  *来自兰溪哥哥的帖子，将格点资料输出为txt格式