[已解决]grads对fnl资料绘制湿位涡，出现点问题，求高手指点

良辰

1. *湿位涡剖面图
   
2. 'reinit'
   
3. 'open F:\taifeng\blw_fnl\fnl.ctl'
   
4. 'enable printF:\taifeng\blw_fnl\xxxpv1pm.gmf'
   
5. 'set z 2 20'
   
6. 'set zlog on'
   
7. 'set lon 100 130'
   
8. 'set lat 20 40'
   
9. 'set t 35'
   
10. 
    
11. *假相当位温
    
12. 'define tc=tmpprs-273.15'
    
13. 'define rh=rhprs'
    
14. 'define prs=lev'  
    
15. 'define es=(6.112*exp((17.67*tc)/(tc+243.5)))'
    
16. 'define qs=0.62197*es/(prs-0.378*es)'
    
17. 'define q=rh*qs/100'
    
18. 'define e=prs*q/(0.62197+q)+1e-10'
    
19. 'define tlcl=55.0+2840.0/(3.5*log(tc+273.16)-log(e)-4.805)'
    
20. 'define theta=(tc+273.16)*pow((1000/prs),(0.2854*(1.0-0.28*q)))'
    
21. 'define eqt=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'
    
22. 
    
23. 
    
24. *湿位涡3
    
25. 'define u=UGRDprs'
    
26. 'define v=VGRDprs'
    
27. 'define vo=hcurl(u,v)'
    
28. 'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
    
29. 'define g=9.8'
    
30. 'define dp=100*(lev(z-1)-lev(z+1))'
    
31. 'define deqt=eqt(z-1)-eqt(z+1)'
    
32. 'define du=u(z-1)-u(z+1)'
    
33. 'define dv=u(z-1)-u(z+1)'
    
34. 'define dx=2.0*6370949.0*cos(lat*3.14159/180.0)*3.14159/180.0'
    
35. 'define dy=2.0*6370949.0*3.14159/180.0'
    
36. 'define dtx=cdiff(eqt,x)'
    
37. 'define dty=cdiff(eqt,y)'
    
38. 'define pv1=-g*(vo+f)*deqt/dp'  
    
39. 'define pv2=g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))'
    
40. 'define pv=pv1+pv2'
    
41. 
    
42. 'set lon 113'
    
43. 'set lat 22 35'
    
44. 'set grid off'
    
45. 'set grads off'
    
46. 'set gxout contour'
    
47. 'set csmooth on'
    
48. 'set cthick 7'
    
49. *'set clopts -1 -1 0.10'
    
50. 'd pv1*10e5'
    
51. *'d pv2*10e6'
    
52. *'d pv*10e5'
    
53. 'print'
    
54. 'disable print'
    
55. ;

复制代码

良辰

问题：1.grads中提示的错误是什么意思？
      2.画出的图上边总少一块，设置过ylev没有效果，怎么解决？
      3.这个算法合理不？是否正确？

mofangbao

错误提示我相信你即使谷歌一下那些英文也能看懂的，就是说enable的地方出现了语法上的错误。很明显，是你enable后面的print和f之间少了一个空格。
少一块应该是你范围设置不对或者数据本身的问题
另外，楼主第一个问题实在是应该自己解决的。先自己仔细思考后再提问比较合适哈。

良辰

mofangbao 发表于 2012-9-28 08:58

                               
登录/注册后可看大图

错误提示我相信你即使谷歌一下那些英文也能看懂的，就是说enable的地方出现了语法上的错误。很明显，是你en ...

那好吧，第一个问题我承认我错了，修改后出现 no hardcopy metafile open ，对grads没有深入了解，希望高手指点！
第二个问题一直还没有解决！

mikado

良辰 发表于 2012-9-28 10:01

                               
登录/注册后可看大图

那好吧，第一个问题我承认我错了，修改后出现 no hardcopy metafile open ，对grads没有深入了解，希望高 ...

第二个问题你应该检查一下FNL资料。我记得RH还是哪个物理量了在FNL里面100hPa以上经常是没有数据的

良辰

mikado 说的有道理，但是不是主要解决办法，主要是由于插值造成的，先进行范围较大的插值计算，绘图时候选择小范围绘图就ok了！谢谢各位！

shuihu88888

为什么我用你的程序画不出图来呢
'reinit'
'open E:\grads\doc\fnl_20110707_00_00.ctl'
'enable print E:\20110708\weiwo.jpg'
'set z 2 20'
'set zlog on'
'set lon 100 130'
'set lat 20 40'
'set t 35'
'define tc=tmpprs-273.15'
'define rh=rhprs'
'define prs=lev'  
'define es=(6.112*exp((17.67*tc)/(tc+243.5)))'
'define qs=0.62197*es/(prs-0.378*es)'
'define q=rh*qs/100'
'define e=prs*q/(0.62197+q)+1e-10'
'define tlcl=55.0+2840.0/(3.5*log(tc+273.16)-log(e)-4.805)'
'define theta=(tc+273.16)*pow((1000/prs),(0.2854*(1.0-0.28*q)))'
'define eqt=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'

'define u=UGRDprs'
'define v=VGRDprs'
'define vo=hcurl(u,v)'
'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
'define g=9.8'
'define dp=100*(lev(z-1)-lev(z+1))'
'define deqt=eqt(z-1)-eqt(z+1)'
'define du=u(z-1)-u(z+1)'
'define dv=u(z-1)-u(z+1)'
'define dx=2.0*6370949.0*cos(lat*3.14159/180.0)*3.14159/180.0'
'define dy=2.0*6370949.0*3.14159/180.0'
'define dtx=cdiff(eqt,x)'
'define dty=cdiff(eqt,y)'
'define pv1=-g*(vo+f)*deqt/dp'  
'define pv2=g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))'
'define pv=pv1+pv2'

'set lon 113'
'set lat 22 35'
'set grid off'
'set grads off'
'set gxout contour'
'set csmooth on'
'set cthick 7'
*'set clopts -1 -1 0.10'
'd pv1*10e5'
*'d pv2*10e6'
*'d pv*10e5'
'print'
'printim E:\20110708\weiwo.jpg white'
'disable print'

river

shuihu88888 发表于 2012-10-20 10:18

                               
登录/注册后可看大图

为什么我用你的程序画不出图来呢
'reinit'
'open E:\grads\doc\fnl_20110707_00_00.ctl'

具体的错误提示之类的呢

良辰

shuihu88888 发表于 2012-10-20 10:18

                               
登录/注册后可看大图

为什么我用你的程序画不出图来呢
'reinit'
'open E:\grads\doc\fnl_20110707_00_00.ctl'

希望这位同学好好看看我的gs,不能盲目照搬的，我的是批量的fnl资料，所以设置了时间变量t，您把那行去掉就ok了！祝你好运！

ppya

应该是画不来的，需要两步操作，我是这样的，pv1可以，pv2不行，不知道有没有高手可以不两步实现。