关于fnl资料 set fwrite 的问题

云淡ぐ风轻

我想利用fnl资料绘制湿位涡，想先通过grads 计算位涡MPV,‘set fwrite’写到二进制文件中，最后配上ctl，再绘图。
发现存在问题，第一：写文件时，提示“data request warning：reguest is completely outside file limits”
第二：出的图（见附件）在700-500hpa存在问题。
有请高手解答。（其中资料就是29个时次的fnl资料）

这个是set fwrite的gs文件
'reinit'
'open D:\grads\data\20130525\fnl_20130523_00_00.ctl'
'set lon 50 160'
'set lat 0 80'
'set z 1 21'
'set t 1 29'
***********************************
  'define a=6.371*1e6'
  'define pi=3.1416/180'
  'define g=9.8'
  'define f=2*7.292*0.00001*sin(lat*3.14159/180)'
  'define vor=hcurl(UGRDprs,VGRDprs)'  
  'define es=6.112*exp(17.67*(TMPprs-273.16)/(TMPprs-29.65))'
  'define r=RHprs*(0.622*es/(lev-0.378*es))/100.0'
  'define e=lev*r/(0.62197+r)+1e-10'
  'define tk=55.0+2840.0/(3.5*log(TMPprs)-log(e)-4.805)'
  'define theta=TMPprs*pow(1000.0/lev,0.286)'
  'define thetase=theta*exp(((3376./tk)-2.54)*r*(1.0+0.81*r))'
  'define dtheZ=thetase(z+1)-thetase(z-1)'
  'define dtheX=cdiff(thetase,x)'
  'define dtheY=cdiff(thetase,y)'
  'define du=UGRDprs(z+1)-UGRDprs(z-1)'
  'define dv=VGRDprs(z+1)-VGRDprs(z-1)'
  'define dp=100*(lev(z+1)-lev(z-1))'
  'define dx=cdiff(lon,x)*3.1416/180'
  'define dy=cdiff(lat,y)*3.1416/180'
  'define mpv0=-g*vor*dtheZ/dp*1e6'
  'define mpv1=-g*(vor+f)*dtheZ/dp*1e6'
  'define mpv2=g*((dv/dp)*(dtheX/dx)-(du/dp)*(dtheY/dy))/(a*cos(lat*pi))*1e6'
  'define mpv=(mpv1+mpv2)'
  'define mmpv=(mpv0+mpv2)'
**********湿位涡*******************
'set gxout fwrite'
'set fwrite D:\grads\data\20130525\pv.dat'
it=1
while(it<=29)
  'set t 'it  
  
iz=1
while(iz<=21)
'set z 'iz
'set lon 60 150'
'set lat 10 70'
'd mpv'
iz=iz+1
endwhile
iz=1
while(iz<=21)
'set z 'iz
'set lon 60 150'
'set lat 10 70'
'd mpv1'
iz=iz+1
endwhile
iz=1
while(iz<=21)
'set z 'iz
'set lon 60 150'
'set lat 10 70'
'd mpv2'
iz=iz+1
endwhile
iz=1
while(iz<=21)
'set z 'iz
'set lon 60 150'
'set lat 10 70'
'd mmpv'
iz=iz+1
endwhile
iz=1
while(iz<=21)
'set z 'iz
'set lon 60 150'
'set lat 10 70'
'd thetase'

iz=iz+1
endwhile

it=it+1  
endwhile
'disable fwrite'
'reinit'
;


输出后配的PV.CTL
dset D:\grads\data\20130525\pv.dat
undef 9.999E+20
title d:/grads/data/20130525/fnl_20130523_00_00
ydef 61 linear 10.000000 1
xdef 91 linear 60.000000 1.000000
tdef 29 linear 00Z23may2013 6hr
zdef 21 levels 1000 975 950 925 900 850 800 775 500 700 650 600 550 500 450 400 350 300 250 200 150 100
vars 5
mpv 21 99 mpvall
mpv1 21 99 mpvall1
mpv2 21 99 mpvall2
mmpv 21 99 mmpvall
thetase 21 99 thetase all
ENDVARS

gs绘图
'reinit'
'open D:\grads\data\20130525\pv.ctl'
  'set lon 100 130'
  'set lat 35'
'set lev 975 250'
'set t 12'
'set gxout contour'
*'set csmooth on'
*'set cterp on'
'd smth9(mpv)'
'printim D:\grads\data\20130525\gmf\pvtest\122mpv.bmp white '
  'print'  
  'c'
  'disable print'
;


贴上图片

传说中的谁

貌似只能计算单个时次单一层次的，你试试看。多时次多层次的要用到数组，具体怎么用我也不知道