ncep资料求湿位涡GS程序运行错误，求高手指点·

林子

'reinit'
'sdfopen d:\ncep\air.2011.nc'
'sdfopen d:\ncep\rhum.2011.nc'
'sdfopen d:\ncep\uwnd.2011.nc'
'sdfopen d:\ncep\vwnd.2011.nc'
'set mpdset cnworld'
'set lat 15 45'
'set lon 95 145'
*'set mproj lambert'
'set ylevs 995 900 800 700 600 500 400 300'
*选lev和z都会出现cannot contour grid-all undefined values的错误
*'set lev 995 300'
'set z 2 20'
'set zlog on'
i=925
while(i<=979)
'set t 'i''
     'q time'
     res = subwrd(result,3)
     year= substr(res,1,12)
     'q time'
     res = subwrd(result,3)
     hour = substr(res,1,2)
     say hour
     'q time'
     res = subwrd(result,3)
     day=substr(res,4,2)
     say day
     'q time'
     res = subwrd(result,3)
     month=substr(res,6,3)
     say month

*****计算假相当位温*****
'set grads off'
'define prs=lev*100'
'define es=100*(6.112*exp(17.67*（air-273.15)/(air-29.65)))'
'define qs=0.62197*es/(prs-es)'
'define qv=rhum.2*qs/100'
'define e=(prs/100.)*qv/(0.62197+qv)+1e-10'
'define tlcl=55.0+2840.0/(3.5*log(air)-log(e)-4.775)'
'undefine e'
'define theta=air*pow((100000./prs),(0.2854*(1.0-0.28*qv)))'
'define eqt=theta*exp(((3376./tlcl)-2.54)*qv*(1.0+0.81*qv))'

*****计算湿位涡3*****
'define vo=hcurl(uwnd.3,vwnd.4)'
'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=uwnd.3(z-1)-uwnd.3(z+1)'
'define dv=vwnd.4(z-1)-vwnd.4(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 121'
'set lat 21 37'
'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'
*'disable print'
'q time'
res=subwrd(result,3)
'draw title 'res'    'pv
'printim d:\ncep\pv\pv1.'month'.'day'.'hour'.500.gif white x1000 y600'
'c'
i=i+1
endwhile

mofangbao

错误就是前面注释的那里么，有没有不保密的数据测试下呢，楼主最好把ctl和出错的详细信息说一下啦

shenyuweisheny

即使你 'set z 2 20'，也不能这样做垂直差分
'define dp=100*(lev(z-1)-lev(z+1))'
'define deqt=eqt(z-1)-eqt(z+1)'
'define du=uwnd.3(z-1)-uwnd.3(z+1)'
'define dv=vwnd.4(z-1)-vwnd.4(z+1)'

必须对垂直层循环，一层一层的计算，不要想一次搞定，如下：
'set gxout fwrite'
'set fwrite mpv_ver_'year'.grd'

zz=1
while(zz<=3)
level=1000-zz*25
lev1=1000-(zz-1)*25
lev2=1000-(zz+1)*25
'set lev 'level''
'define vo=hcurl(UGRD.3,VGRD.4)'
'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
'define g=9.8'
'define dp=5000.0'
'define dthse=thse(lev='lev1')-thse(lev='lev2')'
'define du=UGRD.3(lev='lev1')-UGRD.3(lev='lev2')'
'define dv=VGRD.4(lev='lev1')-VGRD.4(lev='lev2')'
'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(thse,x)'
'define dty=cdiff(thse,y)'
'define mpv1=-g*(vo+f)*dthse/dp'  
'define mpv2=g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))'
'define mpv=mpv1+mpv2'
'd mpv'
zz=zz+1
endwhile

zz=2
while(zz<=17)
level=1000-zz*50
lev1=1000-(zz-1)*50
lev2=1000-(zz+1)*50
'set lev 'level''
'define vo=hcurl(UGRD.3,VGRD.4)'
'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
'define g=9.8'
'define dp=10000.0'
'define dthse=thse(lev='lev1')-thse(lev='lev2')'
'define du=UGRD.3(lev='lev1')-UGRD.3(lev='lev2')'
'define dv=VGRD.4(lev='lev1')-VGRD.4(lev='lev2')'
'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(thse,x)'
'define dty=cdiff(thse,y)'
'define mpv1=-g*(vo+f)*dthse/dp'  
'define mpv2=g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))'
'define mpv=mpv1+mpv2'
'd mpv'
zz=zz+1
endwhile

林子

数据是ncep2.5*2.5资料！

林子

mofangbao 发表于 2012-3-7 23:29

                               
登录/注册后可看大图

错误就是前面注释的那里么，有没有不保密的数据测试下呢，楼主最好把ctl和出错的详细信息说一下啦

ncep 2.5*2.5资料！

abd

请问楼主这个问题解决了吗？麻烦贴出来共享一下好吗？我也在纠结这个湿位涡的计算了....谢谢！

林子

abd 发表于 2012-3-12 09:31

                               
登录/注册后可看大图

请问楼主这个问题解决了吗？麻烦贴出来共享一下好吗？我也在纠结这个湿位涡的计算了....谢谢！

没有呢!不知哪位高人能指点一下!!!

林子

abd 发表于 2012-3-12 09:31

                               
登录/注册后可看大图

请问楼主这个问题解决了吗？麻烦贴出来共享一下好吗？我也在纠结这个湿位涡的计算了....谢谢！

没有呢!不知哪位高人能指点一下!!!

林子

abd 发表于 2012-3-12 09:31

                               
登录/注册后可看大图

请问楼主这个问题解决了吗？麻烦贴出来共享一下好吗？我也在纠结这个湿位涡的计算了....谢谢！

没有呢!不知哪位高人能指点一下!!!

林子

abd 发表于 2012-3-12 09:31

                               
登录/注册后可看大图

请问楼主这个问题解决了吗？麻烦贴出来共享一下好吗？我也在纠结这个湿位涡的计算了....谢谢！

没有呢!不知哪位高人能指点一下!!!