【求助】两个时间段整层水汽通量的t检验出问题了

WOT

这个脚本是可以运行的，但是和文献中的阴影区域差别很大很大，有没有哪位大神解答一下呀！
*************1988-2001与1973-1986七、八月整层水汽通量差**************
'reinit'
'sdfopen E:\GrADS\testwww\uwnd.mon.mean.nc'
'sdfopen E:\GrADS\testwww\vwnd.mon.mean.nc'
'sdfopen E:\GrADS\testwww\shum.mon.mean.nc'
'sdfopen E:\GrADS\testwww\pres.mon.mean.nc'

'set parea 1.5 10.0 1.5 8.0'
'set map 1 1 3'
'set mpdset cnworld cnriver'

'set xlopts -10'
'set ylopts -10'
'set xlint 10'
'set ylint 10'

*****************1988-2001年7、8月整层水汽通量***************
************七月*************************
'set t 487 643'
'set z 1 8 '
'define aqu7=uwnd*shum.3/1000'
'define aqv7=vwnd.2*shum.3/1000'
'set t 415'
'set z 1 8 '
'define aqud7=ave(aqu7,t+0,t=643,12)'
'define aqvd7=ave(aqv7,t+0,t=643,12)'
**************八月***********************
'set t 488 644'
'set z 1 8 '
'define aqu8=uwnd*shum.3/1000'
'define aqv8=vwnd.2*shum.3/1000'
'set t 488'
'set z 1 8 '
'define aqud8=ave(aqu8,t+0,t=644,12)'
'define aqvd8=ave(aqv8,t+0,t=644,12)'
************七、八月平均*********
'set z 1 8 '
'define aqud=(aqud7+aqud8)/2'
'define aqvd=(aqvd7+aqvd8)/2'
******************水汽通量***************
'set z 1'
'define aqus=vint(pres.4(z=1),aqud,500)'
'define aqvs=vint(pres.4(z=1),aqvd,500)'

'define ashuiqi=mag(aqus,aqvs)'


***********************1973-1986年7、8月整层水汽通量**********************************
************七月*************************
'set t 307 463'
'set z 1 8 '
'define bqu7=uwnd*shum.3/1000'
'define bqv7=vwnd.2*shum.3/1000'
'set t 307'
'set z 1 8 '
'define bqud7=ave(bqu7,t+0,t=463,12)'
'define bqvd7=ave(bqv7,t+0,t=463,12)'
**************八月***********************
'set t 308 464'
'set z 1 8 '
'define bqu8=uwnd*shum.3/1000'
'define bqv8=vwnd.2*shum.3/1000'
'set t 308'
'set z 1 8 '
'define bqud8=ave(bqu8,t+0,t=464,12)'
'define bqvd8=ave(bqv8,t+0,t=464,12)'
************七、八月平均*********
'set z 1 8 '
'define bqud=(bqud7+bqud8)/2'
'define bqvd=(bqvd7+bqvd8)/2'
******************水汽通量***************
'set z 1'
'define bqus=vint(pres.4(z=1),bqud,500)'
'define bqvs=vint(pres.4(z=1),bqvd,500)'
'define bshuiqi=mag(bqus,bqvs)'


************************差值***************************************
'define shuiqi=ashuiqi-bshuiqi'
'define qus=aqus-bqus'
'define qvs=aqvs-bqvs'
***********************水汽通量方差**********************
**********************1988-2001年水汽通量方差******
i=487
'set z 1 8 '
'define aaqu7=uwnd*shum.3/1000'
'define aaqv7=vwnd.2*shum.3/1000'
'define assum=aaqu7-aaqu7'
while (i<=643)
'set t ' i
'define aaqu7=uwnd*shum.3/1000'
'define aaqv7=vwnd.2*shum.3/1000'

'set t ' i+1
'set z 1 8 '
'define aaqu8=uwnd*shum.3/1000'
'define aaqv8=vwnd.2*shum.3/1000'

'set z 1 8 '
'define aaqud=(aaqu7+aaqu8)/2'
'define aaqvd=(aaqv7+aaqv8)/2'
'set z 1'
'define aaqus=vint(pres.4(z=1),aaqud,500)'
'define aaqvs=vint(pres.4(z=1),aaqvd,500)'

'define aashuiqi=mag(aaqus,aaqvs)'
'define assum=pow(aashuiqi-ashuiqi,2)+assum'
i=i+12
endwhile
'define avari=assum/14'

*******************1973-1986年水汽通量方差**********
i=307
'set z 1 8 '
'define bbqu7=uwnd*shum.3/1000'
'define bbqv7=vwnd.2*shum.3/1000'
'define bssum=bbqu7-bbqu7'
while (i<=463)
'set t ' i
'define bbqu7=uwnd*shum.3/1000'
'define bbqv7=vwnd.2*shum.3/1000'

'set t ' i+1
'set z 1 8 '
'define bbqu8=uwnd*shum.3/1000'
'define bbqv8=vwnd.2*shum.3/1000'

'set z 1 8 '
'define bbqud=(bbqu7+bbqu8)/2'
'define bbqvd=(bbqv7+bbqv8)/2'
'set z 1'
'define bbqus=vint(pres.4(z=1),bbqud,500)'
'define bbqvs=vint(pres.4(z=1),bbqv7,500)'
'define bbshuiqi=mag(bbqus,bbqvs)'

'define bssum=pow(bbshuiqi-bshuiqi,2)+bssum'

i=i+12
endwhile
'define bvari=bssum/14'

*************************前后两段标准差乘积**************
'define prodvari=sqrt(avari*bvari+shuiqi*(shuiqi-shuiqi))'

***************************t值求算**********************
'define stutest=shuiqi/sqrt((avari+bvari-2*0.53241*prodvari)/13)'

*******************************************************************

'set lon 20 120'
'set lat -20 60'

'set gxout shaded'
*'E:\GrADS\test\wo3color.gs'
*'d shuiqi'
*'E:\GrADS\test\cbar_matlab 1 1 0'
'set rgb 20 255 255 255'
'set rgb 21 104 106 107'
'set rgb 22 66 68 68'

'set rbcols 20 21 22 '
'set clevs 2.16 3.012'

'd abs(stutest)'
*'run  cbarn.gs'

'set gxout barb'
'set ccolor 1'
'set gxout vector'
'set arrscl 0.5 40'
'set arrowhead 0.06'

'd qus;qvs'
*'cbarn 1.2 0'
'run C:\OpenGrADS\Classic\scripts\axis.gs -type b -interval 20 -position o -suffix `3.'
'run C:\OpenGrADS\Classic\scripts\axis.gs -type L -interval 10 -position o -suffix `3.'

'draw title 1000-500hpa vapor flux difference'

'printim E:\GrADS\test\456.png  white x2500 y2000'
'disable print'