做了很久沙氏指数si的计算，最后一个问题请教大神

sunflower

我最近在计算各种指数，其中沙氏指数困扰我太久了，在我把其他指数都算完之后，终于也在闺蜜兼同事的帮助下把沙氏指数也搞得差不多了，但还是有问题……汗，不得不再次来请教诸位大神。
下面我主要把存在问题的循环部分写下来麻烦大家看看到底哪儿出了问题。
前面所有的计算都木有问题，把循环去掉也木有问题，可以确定问题就出在循环上了。
我用的ncep下的grib2资料算的，先算抬升凝结高度处的温度，再根据位温守恒通过迭代计算干绝热上升到500hPa时的温度。
...............
'define t3=tl1-273.15'
'define a1=tk1*pow((1000/p1,(0.2854*(1.0.00028*r1)))*exp((3.376/tl1-0.00254)*r1*(1+0.00081*r1))'
n=0
while(1)
n=n+1
'define t3=t3-0.0000001* 'n''
'tk=t3+273.15'
'define est=6.112*exp((17.67*t3)/(t3+243.5))'
'define etd=6.112*exp((17.67*td2)/(td2+243.5))'
'define t=etd'
'define u=(etd/ets)*100'
'define tl=1/(1/(tk-55)-log(u/100)/2840)+55'
'define r=622*(e/(p2-e))'
'define a2=tk1*pow((1000/p2,(0.2854*(1.0.00028*r)))*exp((3.376/tl-0.00254)*r*(1+0.00081*r))'
'define a=a2-a1'
if(abs(a)<0.001)
'define t=t3'
break
endif
endwhile
'define si=t2-t'
'd si'

sunflower

兰溪之水 发表于 2016-2-2 19:21
你直接d一下a看看命令窗口显示啥

是。全都d过一遍了，a可以出来的，也是对的。后来也算出来si最后的值，就是有些偏大。可能我的while和if还是存在某些拧巴的问题。再次感谢兰溪欧巴

sunflower

不好意思，我的主题分类可能选错，因为选择“脚本编辑”居然发不了

sunflower

有钱钱奖励哦，居然都木有人答复，肿么安心过年……

兰溪之水

你的a应该不是一个数值吧

magua2006

很赞！原来grads也能计算沙氏指数

sunflower

兰溪之水 发表于 2016-2-2 14:00
你的a应该不是一个数值吧

我的a是a2-a1的值，是一个数值，去掉循环是能出结果的。刚刚前面改成while（n<10000）的出沙氏指数的最后结果了，就是感觉不太对，1月的是20多30多，有些偏大。非常感谢兰溪，可能是if的原因，我再仔细看看。

sunflower

magua2006 发表于 2016-2-2 14:56
很赞！原来grads也能计算沙氏指数

目前结果还不太对。嗯哪，grads还是挺强大的，但做完手头的事我还是学习ncl算了

sunflower

兰溪之水 发表于 2016-2-2 14:00
你的a应该不是一个数值吧

难道if里面不能define什么的吗，木有找到类似的例子

sunflower

兰溪之水 发表于 2016-2-2 14:00
你的a应该不是一个数值吧

兰溪哥哥救救我吧……

兰溪之水

你直接d一下a看看命令窗口显示啥