grads假相当位温500-850的gs

ykldmy

本人新手，有没有大神帮我看看，这个gs问题出在哪里？500-850的假相当位温
'reinit'
'sdfopen e:\2010\uwnd.2010.nc'
'sdfopen e:\2010\vwnd.2010.nc'
'sdfopen e:\2010\rhum.2010.nc'
'sdfopen e:\2010\shum.2010.nc'
'sdfopen e:\2010\omega.2010.nc'
'sdfopen e:\2010\air.2010.nc'
'set grid off'
'set grads off'
'set xlopts 1 5 0.18'
'set ylopts 1 5 0.18'
'set lon 80 140'
'set lat 10 50'
'set t 663 675'
'set lev 1000 300'
*--------计算假相当位温----------
'define br=(0.1158*log(rhum.3)-1.5332787)*air.6-789.92816'            


'define cr=789.92816*air.6'                                            


'define tc2=(-br-sqrt(br*br-4*cr))/2.0'                                


'define pc2=lev*pow(tc2/air.6,3.5005574)'                              


'define etc2=6.1078*pow(273.16/tc2,5.1347779)*exp(3148.0973*(tc2-

273.16)/(126.06334*tc2))'   
'define wc2=0.622*etc2/(pc2-etc2)'                                    


'define ed2=tc2*pow(1000.0/(pc2-etc2),0.28765924)'                     


'define lc2=2500.79-2.3697*(tc2-273.16)'                              


'define eqt=ed2*exp(wc2*lc2/1.0048/tc2)'
'define eqt1=eqt(lev500)-eqt(lev850)'

*---------输出图形2------------
'enable print e:\jia2010.gmf'
i=763
while(i<=775)
'set t 'i''
'set mpdset cnworld'
'set gxout contour'
'set cthick 11'
'set clopts -1 -1 0.14'
'set cint 2'
'd eqt1'
'draw title 'i'i'
'print'
'c'
i=i+1
endwhile
'disable print'

river

好歹给个错误提示什么的吧，就这么看我还觉得是对的呢

qixiangren34

恰好我正在弄，发现个问题，我原来是用gs分时分层循环计算假相当位温，输出二进制文件后写ctl,现在用本网站提供的gs计算的结果与原来的图不太一致，位置有偏差，但我更倾向于原来的算法，因为能较好的说明问题。仅提供了个思路，不知能否帮助你。

liutaoz21

前来学习～～～

mofangbao

此种方式提问不靠谱啊~
别人也就只能给你找找语法错误了，甚至都懒得看完的。
楼主先自己把错误位置定好再说

ykldmy

mofangbao 发表于 2013-7-3 08:06

                               
登录/注册后可看大图

此种方式提问不靠谱啊~
别人也就只能给你找找语法错误了，甚至都懒得看完的。
楼主先自己把错误位置定好再 ...

e:\QQ图片

ykldmy

[img]file:///C:\Users\Administrator\AppData\Roaming\Tencent\Users\1123259272\QQ\WinTemp\RichOle\~2IJ4R9F[L`A)8VYBAKNYME.jpg[/img]

ykldmy

river 发表于 2013-7-2 16:17

                               
登录/注册后可看大图

好歹给个错误提示什么的吧，就这么看我还觉得是对的呢

river

ykldmy 发表于 2013-7-15 16:50

                               
登录/注册后可看大图

错误提示说的很明显了，变量错误：需要加减乘除符号的时候却出现一个5，错误的变量名是eqt，那就看看你eqt那一句的等式有没有问题啊，看看是不是符号是中文状态下输入进去的什么的

ykldmy

river 发表于 2013-7-15 21:36

                               
登录/注册后可看大图

错误提示说的很明显了，变量错误：需要加减乘除符号的时候却出现一个5，错误的变量名是eqt，那就看看你eq ...

'reinit'
'sdfopen e:\2013\uwnd.2013.nc'
'sdfopen e:\2013\vwnd.2013.nc'
'sdfopen e:\2013\rhum.2013.nc'
'sdfopen e:\2013\shum.2013.nc'
'sdfopen e:\2013\omega.2013.nc'
'sdfopen e:\2013\air.2013.nc'
'set grid off'
'set grads off'
'set xlopts 1 5 0.18'
'set ylopts 1 5 0.18'
'set lon 60 150'
'set lat 0 60'
'set t 45 62'
'set lev 500'
*--------计算假相当位温----------
'define br=(0.1158*log(rhum.3)-1.5332787)*air.6-789.92816'                              
'define cr=789.92816*air.6'                                                         
'define tc2=(-br-sqrt(br*br-4*cr))/2.0'                                               
'define pc2=lev*pow(tc2/air.6,3.5005574)'                                                        
'define etc2=6.1078*pow(273.16/tc2,5.1347779)*exp(3148.0973*(tc2-273.16)/(126.06334*tc2))'   
'define wc2=0.622*etc2/(pc2-etc2)'                                                        
'define ed2=tc2*pow(1000.0/(pc2-etc2),0.28765924)'                                          
'define lc2=2500.79-2.3697*(tc2-273.16)'                                                  
'define eqt1=ed2*exp(wc2*lc2/1.0048/tc2)'
*------------------
'set lev 850'
*--------计算假相当位温----------
'define br=(0.1158*log(rhum.3)-1.5332787)*air.6-789.92816'                              
'define cr=789.92816*air.6'                                                         
'define tc2=(-br-sqrt(br*br-4*cr))/2.0'                                               
'define pc2=lev*pow(tc2/air.6,3.5005574)'                                                        
'define etc2=6.1078*pow(273.16/tc2,5.1347779)*exp(3148.0973*(tc2-273.16)/(126.06334*tc2))'   
'define wc2=0.622*etc2/(pc2-etc2)'                                                        
'define ed2=tc2*pow(1000.0/(pc2-etc2),0.28765924)'                                          
'define lc2=2500.79-2.3697*(tc2-273.16)'                                                  
'define eqt2=ed2*exp(wc2*lc2/1.0048/tc2)'
*---------输出图形2------------
'enable print e:\2013\jia3.gmf'
i=45
while(i<=62)
'set t 'i''
'set lon 80 140'
'set lat 10 50'
'set gxout contour'
'set cint 5'
'set cthick 11'
'draw title 'i'500-850'
'd eqt1-eqt2'
'print'
'c'
i=i+1
endwhile
'disable print'
;
最后我改成这样就OK了，但是原先的那个我还是不知道哪里有问题，改不好，郁闷死了