关于计算沙氏指数的gs

sunflower

'set lev 850'
'define p1=hgtprs'
'define t1=tmpprs-273.16'
'define rh1=rhprs'
'define td1=t-((14.55+0.114*t1)*(1-0.01*rh1) + pow((2.5+0.007*t1)*(1-0.01*rh1),3)+ (15.9+0.37*t1)*pow((1-0.01*rh),14))'
'define tk1=tmpprs'
'define est1=6.112*exp((17.67*t1)/(t1+243.5))'%计算饱和水汽压
'define etd1=6.112*exp((17.67*td1)/(td1+243.5))'%计算实际水汽压
'define e1=etd1'
'define u1=(etd1/est1)*100'%计算相对湿度
'define tl1=1/(1/(tk1-55)-log(u1/100)/2840)+55'%计算抬升凝结高度的绝对温度（K）
'define r1=622*(e1/(p1-e1))'%计算混合比
'set lev 500'
'define p2=hgtprs'
'define t2=tmpprs-273.16'
'define rh2=rhprs'
'define td2=t-((14.55+0.114*t2)*(1-0.01*rh2) + pow((2.5+0.007*t2)*(1-0.01*rh2),3)+ (15.9+0.37*t2)*pow((1-0.01*rh2),14))'
'define tk2=tmpprs'
'define est2=6.112*exp((17.67*t2)/(t2+243.5))'%计算饱和水汽压
'define etd2=6.112*exp((17.67*td2)/(td2+243.5))'%计算实际水汽压
'define e2=etd2'
'define u2=(etd2/est2)*100'%计算相对湿度
'define tl2=1/(1/(tk2-55)-log(u2/100)/2840)+55'%计算抬升凝结高度的绝对温度（K）
'define r2=622*(e2/(p2-e2))'%计算混合比
'define T3=tl1-273.15'%绝对温度转换为摄氏温度设为气块初始温度
'define a1=(tk1.*((1000/p1)^(0.2854*(1-0.00028*r1)))).*exp((3.376/tl1-0.00254).*r1.*(1+0.00081*r1))' %计算850hpa的相当位温（K）
n=0
while 1
    n=n+1;
    'define T3=T3-0.0000001*n' %用迭代法来得到850hpa气块先沿干绝热线上升，到达抬升凝结高度后再沿湿绝热线上升至500hpa时所具有的气块温度（摄氏度）
    'define 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  e=etd'
    'define  u=(etd/est)*100'
    'define  tl=1/(1/(tk-55)-log(u/100)/2840)+55'
    'define  r=622*(e/(p2-e))'
    'define  a2=(tk.*((1000/p2)^(0.2854*(1-0.00028*r)))).*exp((3.376/tl-0.00254).*r*(1+0.00081*r));%计算500hpa相当位温
    'define  a=a2-a1'
    if abs(a)<0.002 %根据相对位温守恒，当两处的相当位温差小于0.001℃时迭代结束
    *'define  a2=a2'
    'define T=T3'
    break
    endif
  endwhile
'd si=T2-T' %计算沙瓦特指数
        

sunflower

恳求大神指教，实在是折腾太久了……

四叶草

楼主仔细检查下变量名，比如你gs中'define td1=t-((14.55+0.114*t1)*(1-0.01*rh1) + pow((2.5+0.007*t1)*(1-0.01*rh1),3)+ (15.9+0.37*t1)*pow((1-0.01*rh),14))'的t从何来？

river

楼上所说就是要害

sunflower

四叶草 发表于 2015-9-24 17:08
楼主仔细检查下变量名，比如你gs中'define td1=t-((14.55+0.114*t1)*(1-0.01*rh1) + pow((2.5+0.007*t1)*(1 ...

非常感谢指教，我真是眼睛有问题呀，再改改看，再次感谢，中秋快乐哈

sunflower

river 发表于 2015-9-24 20:32
楼上所说就是要害

是嘞，我还真是很粗心呢

sunflower

四叶草 发表于 2015-9-24 17:08
楼主仔细检查下变量名，比如你gs中'define td1=t-((14.55+0.114*t1)*(1-0.01*rh1) + pow((2.5+0.007*t1)*(1 ...

都改回来了，还是不对，也不说错误原因了，累觉不爱

river

sunflower 发表于 2015-9-25 17:29
是嘞，我还真是很粗心呢

'define  a2=(tk.*((1000/p2)^(0.2854*(1-0.00028*r)))).*exp((3.376/tl-0.00254).*r*(1+0.00081*r));%计算500hpa相当位温
还有850假相当位温那一句，里面那些点是干嘛的？
if abs(a)<0.002 %根据相对位温守恒，当两处的相当位温差小于0.001℃时迭代结束  
不是小于0.001吗，为什么写0.002呢

四叶草

sunflower 发表于 2015-9-25 17:28
非常感谢指教，我真是眼睛有问题呀，再改改看，再次感谢，中秋快乐哈

不谢不谢，你也中秋快乐！算法我也不懂，因为之前没接触过，还是得自己多看看相关知识。

sunflower

river 发表于 2015-9-25 20:26
'define  a2=(tk.*((1000/p2)^(0.2854*(1-0.00028*r)))).*exp((3.376/tl-0.00254).*r*(1+0.00081*r));% ...

多谢大神指点，其实我自己也是从一个m文件里改过来的，自身能力不足又疏忽大意，实在是不好意思，我好好改改推敲一下哈，非常感谢