求水汽通量散度时间剖面图的gs，帮忙修改下面的也可以

kaidada

下面的gs，水汽通量散度和实况差太多，其他都还好，请教是不是那个水汽通量散度的公式有问题，哪位大神可以忙修改下，不要纠结我*的2.5哈，那是小问题，万分感谢。
'reinit'
'sdfopen e:\grads\2011.1.nc\uwnd.2011.nc'
'sdfopen e:\grads\2011.1.nc\vwnd.2011.nc'
'sdfopen e:\grads\2011.1.nc\rhum.2011.nc'
'sdfopen e:\grads\2011.1.nc\shum.2011.nc'
'sdfopen e:\grads\2011.1.nc\omega.2011.nc'
'sdfopen e:\grads\2011.1.nc\air.2011.nc'
'set t 803 815'
*-----水平散度------
'set lev 1000 100'
'define div=hdivg(uwnd.1*2.5,vwnd.2*2.5)*1e5'
*-----相对涡度------
'set lev 1000 100'
'define vor=hcurl(uwnd.1*2.5,vwnd.2*2.5)*1e5'
*-----水汽通量散度------
'set lev 1000 300'
'define qdiv=hdivg(uwnd.1*2.5*shum.4*2.5*1e3/9.8,vwnd.2*2.5*shum.4*2.5*1e3)*(1e5)'
*---------水汽通量--------
'define qu=(uwnd.1*shum.4*2.5*1e3/9.8)'
'define qv=(vwnd.2*shum.4*2.5*1e3/9.8)'
*--------计算假相当位温----------
'define br=(0.1158*log(rhum.3)-1.5332787)*air.6-789.92816'                              
'define cr=789.92816*air.6'                                                         
'define tc2=(-br-sqrt(br*br-4*cr))/2.0'                                               
'define pc2=lev*pow(tc2/air.6,3.5005574)'                                                        
'define etc2=6.1078*pow(273.16/tc2,5.1347779)*exp(3148.0973*(tc2-273.16)/(126.06334*tc2))'   
'define wc2=0.622*etc2/(pc2-etc2)'                                                        
'define ed2=tc2*pow(1000.0/(pc2-etc2),0.28765924)'                                          
'define lc2=2500.79-2.3697*(tc2-273.16)'                                                  
'define eqt=ed2*exp(wc2*lc2/1.0048/tc2)'
*--------GRADS设置------------
'set grads off'
'set grid off'
'set csmooth on'
'set map 15 1 9'
'set parea 1 10.2 1.5 8'
'set xlopts 1 6 0.14'
'set ylopts 1 6 0.14'
'set font 0'

*---------输出图形1------------
'enable print e:\grads\20110721.2\qdiv.gmf'
'set lon 108.5'
'set lat 31.75'
'set zlog on'
*'set xlint 4'
'set ylevs 1000 925 850 700 600 500 400 300'
'set xlabs 20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20'
'set gxout shaded'
'set cint 0.5'
'set cmax 0'
'd qdiv'

'set ylpos 0 r'
'set gxout contour'
'set cint 0.5'
'set cthick 11'
*'set clskip 1'
'd qdiv'
'cbarn 1 0.5'
'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 11'
'writehz 1 8.1 (a)水汽通量散度（经度108.5纬度31.75） 1 2.5 1 1.3 0.6 0 11'
'writehz 0.7 1.0 7月20日 1 2 1 1.3 0.6 0 7'
'writehz 2.8 1.0 7月21日 1 2 1 1.3 0.6 0 7'
'writehz 5.9 1.0 7月22日 1 2 1 1.3 0.6 0 7'
'writehz 9.1 1.0 7月23日 1 2 1 1.3 0.6 0 7'
'print'
'disable print'
'printim e:/grads/20110721.2/qdiv.gif white'
'c'
*---------输出图形2------------
'enable print e:\grads\20110721.2\jia.gmf'
'set grads off'
'set lon 108.5'
'set lat 31.75'
'set zlog on'
'set ylevs 1000 925 850 700 600 500 400 300'
'set xlabs 20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20'
'set gxout contour'
'set cint 3'
'set cthick 6'
'set clskip 1'
'd eqt'
'set ylpos 0 r'
'set gxout contour'
'set cint 3'
'set cthick 11'
*'set clskip 1'
'd eqt'
'cbarn 1 0.5'
'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 7'
'writehz 1 8.1 (a)假相当位温（经度108.5纬度31.75） 1 2.5 1 1.3 0.6 0 7'
'writehz 0.7 1.0 7月20日 1 2 1 1.3 0.6 0 7'
'writehz 2.8 1.0 7月21日 1 2 1 1.3 0.6 0 7'
'writehz 5.9 1.0 7月22日 1 2 1 1.3 0.6 0 7'
'writehz 9.1 1.0 7月23日 1 2 1 1.3 0.6 0 7'
'print'
'disable print'
'printim e:/grads/20110721.2/jia.gif white'
'c'
*---------输出图形3------------
'enable print e:\grads\20110721.2\omega.gmf'
'set lev 1000 100'
'set grads off'
'set lon 108.5'
'set lat 31.75'
'set zlog on'
'set ylevs 1000 925 850 700 600 500 400 300 250 200 150 100'
'set xlabs 20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20'
'set gxout shaded'
'set cint 5'
'set cmax 0'
'd omega.5*100'

'set ylpos 0 r'
'set gxout contour'
'set cint 5'
'set cthick 11'
'd omega.5*100'
'cbarn 1 0.5'
'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 7'
'writehz 1 8.1 (a)垂直速度（经度108.5纬度31.75） 1 2.5 1 1.3 0.6 0 7'
'writehz 0.7 1.0 7月20日 1 2 1 1.3 0.6 0 7'
'writehz 2.8 1.0 7月21日 1 2 1 1.3 0.6 0 7'
'writehz 5.9 1.0 7月22日 1 2 1 1.3 0.6 0 7'
'writehz 9.1 1.0 7月23日 1 2 1 1.3 0.6 0 7'
'print'
'disable print'
'printim e:/grads/20110721.2/omega.gif white'
'c'
*---------输出图形4------------
'enable print e:\grads\20110721.2\div.gmf'
'set grads off'
'set lon 108.5'
'set lat 31.75'
'set zlog on'
'set ylevs 1000 925 850 700 600 500 400 300 250 200 150 100'
'set xlabs 20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20'
'set gxout shaded'
'set cint 0.5'
'set cmax 0'
'd div'

'set ylpos 0 r'
'set gxout contour'
'set cint 0.5'
'set cthick 11'
*'set clskip 1'
'd div'
'cbarn 1 0.5'
'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 7'
'writehz 1 8.1 (a)散度（经度108.5纬度31.75） 1 2.5 1 1.3 0.6 0 7'
'writehz 0.7 1.0 7月20日 1 2 1 1.3 0.6 0 7'
'writehz 2.8 1.0 7月21日 1 2 1 1.3 0.6 0 7'
'writehz 5.9 1.0 7月22日 1 2 1 1.3 0.6 0 7'
'writehz 9.1 1.0 7月23日 1 2 1 1.3 0.6 0 7'
'print'
'disable print'
'printim e:/grads/20110721.2/div.gif white'
'c'
*---------输出图形5------------
'enable print e:\grads\20110721.2\vor.gmf'
'set grads off'
'set lon 108.5'
'set lat 31.75'
'set zlog on'
'set ylevs 1000 925 850 700 600 500 400 300 250 200 150 100'
'set xlabs 20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20'
'set gxout shaded'
'set cint 1.5'
'set cmin 0'
'd vor'

'set ylpos 0 r'
'set gxout contour'
'set cint 1.5'
'set cthick 11'
*'set clskip 1'
'd vor'
'cbarn 1 0.5'
'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 7'
'writehz 1 8.1 (a)涡度（经度108.5纬度31.75） 1 2.5 1 1.3 0.6 0 7'
'writehz 0.7 1.0 7月20日 1 2 1 1.3 0.6 0 7'
'writehz 2.8 1.0 7月21日 1 2 1 1.3 0.6 0 7'
'writehz 5.9 1.0 7月22日 1 2 1 1.3 0.6 0 7'
'writehz 9.1 1.0 7月23日 1 2 1 1.3 0.6 0 7'
'print'
'disable print'
'printim e:/grads/20110721.2/vor.gif white'
;

river

和实况是怎么差太多的？有没有图来对比一下。
另外那个2.5是绝对不需要乘的

kaidada

river 发表于 2014-9-17 19:24
和实况是怎么差太多的？有没有图来对比一下。
另外那个2.5是绝对不需要乘的

画了风场，不成严重不符合，乘2.5和实况几乎一样

river

kaidada 发表于 2014-9-17 19:25
画了风场，不成严重不符合，乘2.5和实况几乎一样

风场在画图的时候才需要乘以2.5，这个和grads软件本身有关。计算过程中不需要，否则会引起不必要的误会。

独行侠

shum是相对湿度，不是比湿吧

简君

飘过飘过。我是小小白。

kaidada

river 发表于 2014-9-17 19:43
风场在画图的时候才需要乘以2.5，这个和grads软件本身有关。计算过程中不需要，否则会引起不必要的误会。

嗯嗯，已经发现这个问题了，谢谢提醒

kaidada

独行侠 发表于 2014-11-19 19:49
shum是相对湿度，不是比湿吧

我确实觉得公式有问题，可是又不确定，天气学学的太渣

xiaohezhou

好的很呐。

笑哈乐呵

试试，下载下来