利用NCEP资料画剖面图的一点儿经验

哈小婧

用这个gs画图，为什么总是会出错呢？显示：”enable error:invalid operand“

river

哈小婧 发表于 2013-6-6 11:40

                               
登录/注册后可看大图

用这个gs画图，为什么总是会出错呢？显示：”enable error:invalid operand“

把你的gs贴出来看看吧

哈小婧

'reinit'
'sdfopen c:\nc\uwnd.2013.nc'
'sdfopen  c:\nc\vwnd.2013.nc'
'sdfopen  c:\nc\rhum.2013.nc'
'sdfopen c:\nc\shum.2013.nc'
'sdfopen  c:\nc\omega.2013.nc'
'sdfopen c:\nc\air.2013.nc'
'set t 493 505'
*-----水平散度------
'set lev 1000 100'
'define div=hdivg(uwnd.1,vwnd.2)*1e5'
*-----水汽通量散度------
'set lev 1000 300'
'define qdiv=hdivg(uwnd.1*shum.4*1e3/9.8,vwnd.2*shum.4*1e3/9.8)*(1e5)'
*---------水汽通量--------
'define qu=(uwnd.1*shum.4*1e3/9.8)'
'define qv=(vwnd.2*shum.4*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 c:\nc\qdiv.gmf'
'set lon 105'
'set lat 37'
'set zlog on'
*'set xlint 4'
'set ylevs 1000 925 850 700 600 500 400 300 250 200 150 100'
'set xlabs 08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08'
'set gxout shaded'
'set cint 0.3'
'set cmax 0'
'd qdiv'
'run cbar.gs'
'set ylpos 0 r'
'set gxout contour'
'set cint 0.3'
'set cthick 11'
*'set clskip 1'

'd qdiv'
'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 11'
'writehz 1 8.1 水汽通量散度（经度105纬度37） 1 2.5 1 1.3 0.6 0 11'
'writehz 0.7 1.0 5月4日 1 2 1 1.3 0.6 0 7'
'writehz 2.8 1.0 5月5日 1 2 1 1.3 0.6 0 7'
'writehz 5.9 1.0 5月6日 1 2 1 1.3 0.6 0 7'
'writehz 9.1 1.0 5月7日 1 2 1 1.3 0.6 0 7'

'print'
'disable print'
'printim c:/nc/qdiv.gif white'
'c'
*---------输出图形2------------
'enable print c:\nc\jia.gmf'
'set grads off'
'set lon 105'
'set lat 37'
'set zlog on'
'set ylevs 1000 925 850 700 600 500 400 300 250 200 150 100'
'set xlabs 08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08'
'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'

'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 7'
'writehz 1 8.1 假相当位温（经度105纬度37） 1 2.5 1 1.3 0.6 0 7'
'writehz 0.7 1.0 5月4日 1 2 1 1.3 0.6 0 7'
'writehz 2.8 1.0 5月5日 1 2 1 1.3 0.6 0 7'
'writehz 5.9 1.0 5月6日 1 2 1 1.3 0.6 0 7'
'writehz 9.1 1.0 5月7日 1 2 1 1.3 0.6 0 7'

'print'
'disable print'
'printim c:/nc/jia.gif white'
'c'
*---------输出图形3------------
'enable printc:\nc\omega.gmf'
'set lev 1000 100'
'set grads off'
'set lon 105'
'set lat 37'
'set zlog on'
'set ylevs 1000 925 850 700 600 500 400 300 250 200 150 100'
'set xlabs 08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08'
'set gxout shaded'
'set cint 0.05'
'set cmax 0'
'd omega.5'
'run cbar.gs'
'set ylpos 0 r'
'set gxout contour'
'set cint 0.05'
'set cthick 11'
'd omega.5'
'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 7'
'writehz 1 8.1 垂直速度（经度105纬度37） 1 2.5 1 1.3 0.6 0 7'
'writehz 0.7 1.0 5月4日 1 2 1 1.3 0.6 0 7'
'writehz 2.8 1.0 5月5日 1 2 1 1.3 0.6 0 7'
'writehz 5.9 1.0 5月6日 1 2 1 1.3 0.6 0 7'
'writehz 9.1 1.0 5月7日 1 2 1 1.3 0.6 0 7'
'print'
'disable print'
'printim c:/nc/omega.gif white'
'c'
*---------输出图形4------------
'enable c:\nc\div.gmf'
'set grads off'
'set lon 105'
'set lat 37'
'set zlog on'
'set ylevs 1000 925 850 700 600 500 400 300 250 200 150 100'
'set xlabs 08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08'
'set gxout shaded'
'set cint 0.5'
'set cmax 0'
'd div'
'run cbar.gs'
'set ylpos 0 r'
'set gxout contour'
'set cint 0.5'
'set cthick 11'
*'set clskip 1'
'd div'
'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 7'
'writehz 1 8.1 散度（经度105纬度37） 1 2.5 1 1.3 0.6 0 7'
'writehz 0.7 1.0 5月4日 1 2 1 1.3 0.6 0 7'
'writehz 2.8 1.0 5月5日 1 2 1 1.3 0.6 0 7'
'writehz 5.9 1.0 5月6日 1 2 1 1.3 0.6 0 7'
'writehz 9.1 1.0 5月7日 1 2 1 1.3 0.6 0 7'
'print'
'disable print'
'printim c:/nc/div.gif white'
'c'


我就是用楼主你给的这个  画着看看~~  是不是哪里改错了呢。。。。

哈小婧

river 发表于 2013-6-6 12:44

                               
登录/注册后可看大图

把你的gs贴出来看看吧

你看一下  我贴出来了~~

river

哈小婧 发表于 2013-6-6 21:45

                               
登录/注册后可看大图

'reinit'
'sdfopen c:\nc\uwnd.2013.nc'
'sdfopen  c:\nc\vwnd.2013.nc'

你这个只是改了路径吧，别的感觉都没动呢。你把你出图的路径改到其他盘试试，有可能是c盘有保护。

哈小婧

river 发表于 2013-6-6 22:14

                               
登录/注册后可看大图

你这个只是改了路径吧，别的感觉都没动呢。你把你出图的路径改到其他盘试试，有可能是c盘有保护。

是啊 就是想用这个画着试一下~~~~好的 我改改路径看~！~谢谢哈~~~

油麦面疙瘩

楼主，我想问一下，nc文件里的比湿的单位是kg/kg么？为什么感觉我用nc资料画出来的水汽通量散度图和用ncep的grib1的资料画出来的差了一个量级。。。

river

油麦面疙瘩 发表于 2013-6-7 20:55

                               
登录/注册后可看大图

楼主，我想问一下，nc文件里的比湿的单位是kg/kg么？为什么感觉我用nc资料画出来的水汽通量散度图和用ncep的 ...

是0.00001 kg/kg，你去官网上看看说明吧，我也好长时间了，记不清了

油麦面疙瘩

river 发表于 2013-6-7 21:31

                               
登录/注册后可看大图

是0.00001 kg/kg，你去官网上看看说明吧，我也好长时间了，记不清了

嗯，我看过说明了，就是0.00001kg/kg，所以才搞不懂，那这个gs文件里的比湿最后得出的结果，单位是什么？如果要让单位是g/kg的话，不是得乘以1e8么？这样感觉量级好像总不对。。。

river

油麦面疙瘩 发表于 2013-6-8 09:37

                               
登录/注册后可看大图

嗯，我看过说明了，就是0.00001kg/kg，所以才搞不懂，那这个gs文件里的比湿最后得出的结果，单位是什么？ ...

是呀，gs里面不是先乘了1e3，最后又乘的1e5么