计算沙氏指数si过程中遇到的问题~

南南偶卖糕Oo

各位前辈们好，我最近需要计算出沙氏指数的值，参考了许多做过相关内容的前辈的帖子，也读了文献，根据这些帖子和文献写了一个gs。问题就是循环部分，while（n<7）这里，选择7是因为文献中说，一般迭代循环6次左右就差不多可以保证精确度，但后来发现无论n取什么值，得出的结果都是一样的，而且循环出来的值都偏大（有十几二十几甚至三十几），且几乎为正（很明显是错误的）。而且if语句似乎也没有起到作用。算了好几天了，很急又很崩溃，希望各位大佬帮我研究研究，找出问题所在，不胜感激！参考的文献如下：

[2015气象]沙氏指数的计算.pdf (861.33 KB, 下载次数: 24)

2019-4-10 18:52 上传

点击文件名下载附件

我的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格式
;

南南偶卖糕Oo

给自己顶一个防沉

astiny

所以大佬你现在成了吗？

南南偶卖糕Oo

astiny 发表于 2019-5-24 11:31
所以大佬你现在成了吗？

很抱歉唉，我放弃了，做出来全是正值，而且数据很大，如果你做出来一定要告诉我哦！

astiny

南南偶卖糕Oo 发表于 2019-5-31 22:40
很抱歉唉，我放弃了，做出来全是正值，而且数据很大，如果你做出来一定要告诉我哦！

http://bbs.06climate.com/forum.p ... p;extra=#pid1013574
瞎写的。你将就着看吧。欢迎试用&指正

南南偶卖糕Oo

astiny 发表于 2019-6-3 11:00
http://bbs.06climate.com/forum.php?mod=viewthread&tid=90981&page=1&extra=#pid1013574
瞎写的。你将 ...

感谢感谢啦！

ExpectantSam

我用这篇文章的Bolton公式计算的850hPa假相当位温，和怀俄明网站上的对不上，差了40-50左右，反而和相当位温差别不大……