请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 27352|回复: 85

[分享资料] 做了很久沙氏指数si的计算,最后一个问题请教大神

[复制链接]

新浪微博达人勋

发表于 2016-1-28 09:39:48 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
我最近在计算各种指数,其中沙氏指数困扰我太久了,在我把其他指数都算完之后,终于也在闺蜜兼同事的帮助下把沙氏指数也搞得差不多了,但还是有问题……汗,不得不再次来请教诸位大神。
下面我主要把存在问题的循环部分写下来麻烦大家看看到底哪儿出了问题。
前面所有的计算都木有问题,把循环去掉也木有问题,可以确定问题就出在循环上了。
我用的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'
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-2-3 08:45:56 | 显示全部楼层
兰溪之水 发表于 2016-2-2 19:21
你直接d一下a看看命令窗口显示啥

是。全都d过一遍了,a可以出来的,也是对的。后来也算出来si最后的值,就是有些偏大。可能我的while和if还是存在某些拧巴的问题。再次感谢兰溪欧巴
密码修改失败请联系微信:mofangbao
回复 支持 0 反对 2

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2016-1-28 09:40:39 | 显示全部楼层
不好意思,我的主题分类可能选错,因为选择“脚本编辑”居然发不了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-2-2 09:52:10 | 显示全部楼层
有钱钱奖励哦,居然都木有人答复,肿么安心过年……
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 成长值: 19710
发表于 2016-2-2 14:00:45 | 显示全部楼层

回帖奖励 +1 金钱

你的a应该不是一个数值吧
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-2-2 14:56:19 | 显示全部楼层

回帖奖励 +1 金钱

很赞!原来grads也能计算沙氏指数
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-2-2 14:58:59 | 显示全部楼层
兰溪之水 发表于 2016-2-2 14:00
你的a应该不是一个数值吧

我的a是a2-a1的值,是一个数值,去掉循环是能出结果的。刚刚前面改成while(n<10000)的出沙氏指数的最后结果了,就是感觉不太对,1月的是20多30多,有些偏大。非常感谢兰溪,可能是if的原因,我再仔细看看。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-2-2 15:00:44 | 显示全部楼层
magua2006 发表于 2016-2-2 14:56
很赞!原来grads也能计算沙氏指数

目前结果还不太对。嗯哪,grads还是挺强大的,但做完手头的事我还是学习ncl算了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-2-2 15:13:30 | 显示全部楼层
兰溪之水 发表于 2016-2-2 14:00
你的a应该不是一个数值吧

难道if里面不能define什么的吗,木有找到类似的例子
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-2-2 15:23:42 | 显示全部楼层
兰溪之水 发表于 2016-2-2 14:00
你的a应该不是一个数值吧

兰溪哥哥救救我吧……
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 成长值: 19710
发表于 2016-2-2 19:21:06 | 显示全部楼层
你直接d一下a看看命令窗口显示啥
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表