画湿位涡搞不定了，请哪位高手帮我看看是哪里出了问题？

xiaoyuyly

我本来是想画湿位涡的垂直剖面，可是一直出不来，参考了论坛里很多人的帖子，结果都运行不通，好吧，那我先画一层的，我画700hpa的湿位涡，可还是出不来，我调试了假相当位温，是可以出来的（垂直层），因为我用的是FNL 1*1的资料，所以层数为26层，我要画700hpa的话是第9层，所以我对z取8-10层，为了是用上下两层来做差分，不知道是哪里的问题，还是出不来，请高手帮我指点指点吧，为了这个我都快发狂了，错误提示是“cannot contour grid-all undefined values”，请高手一定给以指点，太感谢了，{:soso_e109:}{:soso_e109:}{:soso_e109:}{:soso_e109:}{:soso_e109:}{:soso_e109:}{:soso_e109:}{:soso_e109:}{:soso_e109:}
'reinit'
'open e:\2012\ncep\20130823\2013082412.ctl'
'set grads off'
'set mpdset yunnan'
'set lon 70 150'
'set lat 15 55'
'set ylevs 1000 975 950 925 900 850 800 750 700 650 600 550 500 450 400 350 300 250 200 150 100 70 50 30 20 10'
'set z 8 10'
'set zlog on'
*计算广义位温*
'define prs=lev'
'define aa=TMPprs-273.15'
'define es=(6.112*exp((17.67*aa)/(aa+243.5)))'
'define qs=0.62197*es/(prs-0.378*es)'
'define q=RHprs*qs/100'
'define e=prs*q/(0.62197+q)+1e-10'
'define tlcl=55.0+2840.0/(3.5*log(aa+273.16)-log(e)-4.805)'
'define theta=(aa+273.16)*pow((1000/prs),(0.2854*(1.0-0.28*q)))'
'define theta1=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'
*****计算700hPa的广义湿位涡：
'set z 9'
'define Rd=287'
'define P=70000'
'define a=Rd/P*theta*pow(P/100000,0.286)'
***算好比容
'define vor=hcurl(ugrdprs,vgrdprs)'
'define rr=6.37e6'
***地球半径
'define dx=cdiff(lon,x)*3.14/180*cos(lat*3.14/180)*rr'
'define dy=cdiff(lat,y)*3.14/180*rr'
***差分时，x和y方向的实际距离
     'define g=9.8'
     'define zizhuan=7.292e-5'
***地球自转角速度
     'define ee=2*zizhuan*cos(lat*3.14/180)'
     'define f=2*zizhuan*sin(lat*3.14/180)'
****差分时，z方向的距离
***算700hPa广义湿位涡
     'define wz=-a*vvelprs/g'
***z坐标下的垂直速度**********此处进行了修改，dz需要放在固定层的后面
'define dz=(hgtprs(z+1)-hgtprs(z-1))'
'define pvx=a*(cdiff(wz,y)/dy-(vgrdprs(z+1)-vgrdprs(z-1))/dz)*cdiff(theta1,x)/dx'
'define pvy=a*(ee+(ugrdprs(z+1)-ugrdprs(z-1))/dz-cdiff(wz,x)/dx)*cdiff(theta1,y)/dy'
'define pvz=a*(f+cdiff(vgrdprs,x)/dx-cdiff(ugrdprs,y)/dy)*(theta1(z+1)-theta1(z-1))/dz'
'define pv=pvx+pvy+pvz'
****上边是广义湿位涡的x、y、z上的分量，算好相加即为广义湿位涡
'set grid off'
'set grads off'
'set xlopts 1 6 0.25'
'set ylopts 1 6 0.25'
'set xlint 2'
'set ylint 100'
'set cthick 5'
'set cstyle 3'
'set cmark 4'
'set clopts -1 -1 0.2'
'set lat 18 32'
'set lon 94 110'
*'set lev  700'
'set gxout shaded'
*'set ccolor rainbow'
*'d theta1'
*'d skip(u,2);skip(v,2)'
'd pv*1e6'        
*'cnbasemap (pvx+pvy+pvz)*1e6'
'C:\OpenGrADS\Contents\Resources\Scripts\cbarn.gs'
'cbarn.gs'
'printim e:\2012\ncep\20130823\2013082412swo1.gif gif x800 y600 white'
'enable print  e:\2012\ncep\20130823\2013082412swo1.gmf'
'disable print'

river

你自己不是知道要做差分么，但是你看看你的差分放在哪了？
放在set z 9的后面，那你自己把层次定死了，后面的差分不就错了么

xiaoyuyly

river 发表于 2013-10-24 16:26

                               
登录/注册后可看大图

你自己不是知道要做差分么，但是你看看你的差分放在哪了？
放在set z 9的后面，那你自己把层次定死了，后面 ...

可是我改了以后，把SET Z 9 放在 'define wz=-a*vvelprs/g' 的前面，还是一样的错误出现啊，真的是奇怪了，我要发疯了！！

xiaoyuyly

求哪个高手帮我指点解决一下吧，我为了画这个湿位涡花了我好几天的功夫，网上搜到的所有帖子的gs我都基本试过了，可是在我这里就是调不通，真的是见鬼了！！
哪个帮我解决我万分万分感谢！

river

xiaoyuyly 发表于 2013-10-24 16:34

                               
登录/注册后可看大图

可是我改了以后，把SET Z 9 放在 'define wz=-a*vvelprs/g' 的前面，还是一样的错误出现啊，真的是奇怪了 ...

'define wz=-a*vvelprs/g'这句后面没有差分了？你自己的脚本怎么自己都不熟悉到这个地步(⊙o⊙)？

xiaoyuyly

river 发表于 2013-10-24 16:44

                               
登录/注册后可看大图

'define wz=-a*vvelprs/g'这句后面没有差分了？你自己的脚本怎么自己都不熟悉到这个地步(⊙o⊙)？

我是菜鸟，请不要见怪，是还有差分，呵呵，我刚学习grads不久，这种涉及到差分的作图和计算还是第一次做，所以还真的不太知道，请大侠您告诉我怎么改可以吗？万分万分万分感谢！

river

xiaoyuyly 发表于 2013-10-24 16:49

                               
登录/注册后可看大图

我是菜鸟，请不要见怪，是还有差分，呵呵，我刚学习grads不久，这种涉及到差分的作图和计算还是第一次做， ...

这个不能说怪谁，grads的编程论坛里也就那么几个大神比较精通，比如清风，比如兰溪。而且差分确实平时用的也比较少。所以一般很少有人可以直接看出来你的脚本运行错误根本原因，都是一步一步排除的。首先你要把一些显而易见的问题改正，然后再慢慢调试。你这个你自己也说了为了差分，设置了z取8-10层，但是你的脚本里所有的差分都是在set Z 9的后面，也就是说你固定了Z ，那差分的运算就无法进行，所以报错了

north23

不知搞定没有啊  

hxyj

帖子里提到的广义位温，是高守亭老师提出的吗？如果是，感觉计算公式有问题。

就要唯一

相同的问题，困扰我好几天了