关于湿位涡的画图，请大神帮忙！

蓝岚

脚本我是在论坛中找到的，按照需要更改了，但图有错，请帮忙？
***假相当位温***
'reinit'
'set gxout fwrite'
'set fwrite h:\thse.dat'
'sdfopen h:\air.2016.nc'
'sdfopen h:\rhum.2016.nc'
'sdfopen h:\uwnd.2016.nc'
'sdfopen h:\vwnd.2016.nc'
i=80
while(i<=110)
'set lon 40 160'
'set lat 20 90'
'set t 'i''
'set lev 925'
'define prs=lev*100'
'define ms=17.67*(air-273.15)/(air-29.65)'
'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 thse=theta*exp(((3376./tlcl)-2.54)*qv*(1.0+0.81*qv))'
'd thse'
'set lev 850'
'define prs=lev*100'
'define ms=17.67*(air-273.15)/(air-29.65)'
'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 thse=theta*exp(((3376./tlcl)-2.54)*qv*(1.0+0.81*qv))'
'd thse'
'set lev 700'
'define prs=lev*100'
'define ms=17.67*(air-273.15)/(air-29.65)'
'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 thse=theta*exp(((3376./tlcl)-2.54)*qv*(1.0+0.81*qv))'
'd thse'
'set lev 600'
'define prs=lev*100'
'define ms=17.67*(air-273.15)/(air-29.65)'
'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 thse=theta*exp(((3376./tlcl)-2.54)*qv*(1.0+0.81*qv))'
'd thse'
'set lev 500'
'define prs=lev*100'
'define ms=17.67*(air-273.15)/(air-29.65)'
'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 thse=theta*exp(((3376./tlcl)-2.54)*qv*(1.0+0.81*qv))'
'd thse'
'set lev 400'
'define prs=lev*100'
'define ms=17.67*(air-273.15)/(air-29.65)'
'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 thse=theta*exp(((3376./tlcl)-2.54)*qv*(1.0+0.81*qv))'
'd thse'
   j=2
   while(j<=7)
   'set lon 40 160'
   'set lat 20 90'
   'set t 'i''
   'set z 'j''
   'd uwnd.3'
    j=j+1
    endwhile
    j=2
   while(j<=7)
   'set lon 40 160'
   'set lat 20 90'
   'set t 'i''
   'set z 'j''
   'd vwnd.4'
    j=j+1
    endwhile
i=i+1
endwhile
'disable fwrite'
'reinit'
****配合ctl文件****
DSET h:\thse.dat
TITLE Upper Data
undef -9.99E33
xdef 49 linear 40 2.5
ydef 29 linear 20 2.5
zdef 6 levels 925 850 700 600 500 400
tdef 31 linear 18z20JAN2016 6hr
vars 3
thse 6 99 temps
uwnd 6 99 u-wind
vwnd 6 99 v-wind
Endvars
***计算湿位涡****
'reinit'
'set gxout fwrite'
'set fwrite h:\mpv.dat'
'open h:\thse.ctl'
i=1
while(i<=31)
'set lon 40 160'
'set lat 20 90'
'set t 'i''
'set lev 850'
'define vo=hcurl(uwnd,vwnd)'
'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
'define g=9.8'
'define dp=100*(925-700)'
'define dthse=thse(lev=925)-thse(lev=700)'
'define du=uwnd(lev=925)-uwnd(lev=700)'
'define dv=vwnd(lev=925)-vwnd(lev=700)'
'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 mpv1*1000000'
'set lev 700'
'define vo=hcurl(uwnd,vwnd)'
'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
'define g=9.8'
'define dp=100*(850-600)'
'define dthse=thse(lev=850)-thse(lev=600)'
'define du=uwnd(lev=850)-uwnd(lev=600)'
'define dv=vwnd(lev=850)-vwnd(lev=600)'
'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 mpv1*1000000'
'set lev 600'
'define vo=hcurl(uwnd,vwnd)'
'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
'define g=9.8'
'define dp=100*(700-500)'
'define dthse=thse(lev=700)-thse(lev=500)'
'define du=uwnd(lev=700)-uwnd(lev=500)'
'define dv=vwnd(lev=700)-vwnd(lev=500)'
'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 mpv1*1000000'
'set lev 500'
'define vo=hcurl(uwnd,vwnd)'
'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
'define g=9.8'
'define dp=100*(600-400)'
'define dthse=thse(lev=600)-thse(lev=400)'
'define du=uwnd(lev=600)-uwnd(lev=400)'
'define dv=vwnd(lev=600)-vwnd(lev=400)'
'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 mpv1*1000000'

'set lev 850'
'define vo=hcurl(uwnd,vwnd)'
'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
'define g=9.8'
'define dp=100*(925-700)'
'define dthse=thse(lev=925)-thse(lev=700)'
'define du=uwnd(lev=925)-uwnd(lev=700)'
'define dv=vwnd(lev=925)-vwnd(lev=700)'
'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 mpv2*1000000'
'set lev 700'
'define vo=hcurl(uwnd,vwnd)'
'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
'define g=9.8'
'define dp=100*(850-600)'
'define dthse=thse(lev=850)-thse(lev=600)'
'define du=uwnd(lev=850)-uwnd(lev=600)'
'define dv=vwnd(lev=850)-vwnd(lev=600)'
'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 mpv2*1000000'
'set lev 600'
'define vo=hcurl(uwnd,vwnd)'
'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
'define g=9.8'
'define dp=100*(700-500)'
'define dthse=thse(lev=700)-thse(lev=500)'
'define du=uwnd(lev=700)-uwnd(lev=500)'
'define dv=vwnd(lev=700)-vwnd(lev=500)'
'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 mpv2*1000000'
'set lev 500'
'define vo=hcurl(uwnd,vwnd)'
'define f=2*7.292*sin(lat*3.14159/180.0)*0.00001'
'define g=9.8'
'define dp=100*(600-400)'
'define dthse=thse(lev=600)-thse(lev=400)'
'define du=uwnd(lev=600)-uwnd(lev=400)'
'define dv=vwnd(lev=600)-vwnd(lev=400)'
'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 mpv2*1000000'
i=i+1
endwhile
'disable fwrite'
'reinit'
***配合ctl文件****
DSET h:\mpv.dat
TITLE Upper Data
undef -9.99E33
xdef 49 linear 40 2.5
ydef 29 linear 20 2.5
zdef 4 levels 850 700 600 500  
tdef 31 linear 18z20JAN2016 6hr
vars 2
mpv1 4 99 mpv-1
mpv2 4 99 mpv-2
Endvars
****进行画图****
'reinit'
'open h:\mpv.ctl'
i=1
while(i<=31)
'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 lon 40 160'
'set lat 20 90'
'set lev 850'
'set gxout contour'
'set csmooth on'
'set cthick 7'
'set clopts -1 -1 0.10'
'set grads off'
'd mpv1'
'q time'
res=subwrd(result,3)
'draw title 'res'    'mpv1
'printim h:\pv1.'month'.'day'.'hour'.400.gif white x1000 y600'
'c'
i=i+1
endwhile
'reinit'
****图形如下***

river

可能缺测值有问题吧，改改缺测值试试。最好用其他软件把数据读出来看看缺测值是哪个

初翠翠

哇哦，简直棒棒的，好全