编辑沙氏指数脚本出错，同求绘制沙氏指数脚本

大力

要绘制沙氏指数平面图，先通过计算单点沙氏指数，但总是不对。同时，求计算沙氏指数脚本脚本如下：
'reinit'
'open  D:\data\analysis\20170802\2017080100.ctl'

'set lon 116'
'set lat 40'

'set t 5'

'set lev 850'
'define t850=TMPprs'
************************求850hPa露点温度**************
'define a=(7.5*t850)/(273.3+t850)'
'define es=6.112*pow(10,a)'
'define e=es*RHprs/100'
'define td850=273.3*log10(e/6.112)/(7.5-log10(e/6.112))'

************************求抬升凝结高度和温度**************
'define tc=t850-(t850-td850)*0.976/(0.976-0.000833*pow((237.3+td850),2)/273+td850)'
'define b=(273+tc)/(273+td850)'
'define pc=850*pow(b,3.5)'
*if(tc>=273.15)
'define esc=6.112*pow(10,7.5*tc/(273.3+tc))'
*else
*esc=6.112*pow(10,9.5*tc/(265.5+tc))
*endif

'define ec=esc*RHprs(lev='pc')/100'

if(tc-0.1>=273.15)
esc2=6.112*pow(10,(7.5*(tc-0.1)/(273.3+(tc-0.1))))
else
esc2=6.112*pow(10,(9.5*(tc-0.1)/(265.5+(tc-0.1))))
endif
'define ec2='esc2'*RHprs(lev='pc')/100'

**********求湿绝热过程至500hPa的温度，利用假相当位温守恒，采用二分法*****
'define rc=0.622*ec/(pc-ec)'
'define rc1=0.622*ec/((pc-10)-ec)'
'define rc2=0.622*ec2/((pc-10)-ec2)'
'define thse=tc*pow(1000/(pc-ec),0.286)*exp(2500000*rc/(1005*tc)'
'define thse1=tc*pow(1000/((pc-10)-ec),0.286)*exp(2500000*rc1/(1005*tc)'
'define thse2=(tc-0.1)*pow(1000/((pc-10)-ec2),0.286)*exp(2500000*rc2/(1005*(tc-0.1))'
'define t500=tc-(thse1-thse)*0.1/(thse1-thse2)'

'set lev 500'
'define tmp500=TMPprs'

'define si=tmp500-t500'

'd si'

红尘滚滚1

非常感谢楼主的内容，受教了！！！！！！！

大力

红尘滚滚1 发表于 2017-11-3 10:44
非常感谢楼主的内容，受教了！！！！！！！

但我的不对啊，没有算出来

南南偶卖糕Oo

请问楼主有做出来吗？可以分享一下gs嘛？QAQ谢谢楼主救命！