1957年-2001年多年夏季平均水汽通量，我这样写不知道对不对

洗雨

刚开始学习画图，很多都忘记了不知道我的做法对不对！？
这个是读取比湿的gs 思路是读取全球的夏季平均资料，数据共8层，45年。
'reinit'
'sdfopen g:\mon.mean\shum.mon.mean.nc'
'set x 1 144'
'set y 1 73'
'set gxout fwrite'
'set fwrite g:\shsunor.dat'
z=1,
mt=109;while(mt<=648)
'set z 'z
'define aa=ave(shum,t='6+mt',t=8+mt')'
'd aa'
mt=mt+12;endwhile
z=z+1;endwhile
'disable fwrite'
'reinit'
提取地表气压
'reinit'
'sdfopen G:\mon.mean\pres.sfc.mon.mean.nc'
'set x 1 144'
'set y 1 73'
'set gxout fwrite'
'set fwrite g:\presunor.dat'
mt=109;while(mt<=648)
'define aa=ave(pres,t='6+mt',t='8+mt')'
'd aa'
mt=mt+12;endwhile
'disable fwrite'
'reinit'
提取经向风场
'reinit'
'sdfopen g:\mon.mean\vwnd.mon.mean.nc'
'set x 1 144'
'set y 1 73'
'set gxout fwrite'
'set fwrite g:\vwsunor.dat'

mt=109;while(mt<=648)
z=1;while(z<=17)
'set z 'z
'define se=ave(vwnd,t='6+mt',t='8+mt')'
'd se'
z=z+1;endwhile
mt=mt+12;endwhile
'disable fwrite'
'reinit'
提取纬向风场
'reinit'
'sdfopen g:\mon.mean\uwnd.mon.mean.nc'
'set x 1 144'
'set y 1 73'
'set gxout fwrite'
'set fwrite g:\uwsunor.dat'

mt=109;while(mt<=648)
z=1;while(z<=17)
'set z 'z
'define b=ave(uwnd,t='6+mt',t='8+mt')'
'd b'
z=z+1;endwhile
mt=mt+12;endwhile
'disable fwrite'
'reinit'
以下为比湿的ctl，这里比湿变量我用a表示
DSET g:\shsunor.dat
TITLE  Summer shum  ave data
UNDEF -9.99E33
XDEF 144 LINEAR 0 2.5
YDEF 73 LINEAR -90 2.5
ZDEF 8 LEVELS 1000 925 850 700 600 500 400 300
TDEF 45 LINEAR 1Sep1957 1yr
vars 1
a 8 99 shum at summer
endvars
地表气压ctl
DSET g:\presunor.dat
TITLE pres long term mean
UNDEF -9.99E33
XDEF 144 LINEAR 0 2.5
YDEF 73 LINEAR -90 2.5
ZDEF 1 LEVELS 1000
TDEF 45 LINEAR 1jun1957 1yr
vars 1
pres 0 99 shum at summer
endvars

uwnd的ctl
DSET g:\uwsunor.dat
TITLE  Summer shum  ave data
UNDEF -9.99E33
XDEF 144 LINEAR 0 2.5
YDEF 73 LINEAR -90 2.5
ZDEF 17 LEVELS 1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10
TDEF 45 LINEAR 1jun1957 1yr
vars 1
b 17 99 uwnd at summer
endvars

vwnd的ctl
DSET g:\vwsunor.dat
TITLE  Summer shum  ave data
UNDEF -9.99E33
XDEF 144 LINEAR 0 2.5
YDEF 73 LINEAR -90 2.5
ZDEF 17 LEVELS 1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10
TDEF 45 LINEAR 1jun1957 1yr
vars 1
se 17 99 uwnd at summer
endvars

画图的gs如下
思路是先算一个时次的的水汽通量，然后再在时间上ave，和别人画的图不合
'reinit'
'open g:\gs\shum.ctl'
'open g:\gs\uwnd.ctl'
'open g:\gs\vwnd.ctl'
‘open g:\gs\pres.ctl'
'enable print g:\picture\1957-2001Aavequv.gmf'
'set lon 60 160'
'set lat -20 50'
'set grads off'
'set grid off'
'set t 1'
'set xlint 5'
'set ylint 5'
'set mpdset cnworld'
'set map 1 1 1'
'draw map'
'set xlopts 1 1 0.15'
'set ylopts 1 6 0.15'
'set clopts 1 1 0.22'
'run define_colors.gs'

'set gxout shaded'
*'set ccols 45 44 43 42 0'
*'set csmooth on'
*'set clevs  -6 -4 -2 0'
'define u=vint(lpres.4,b.2*a.1,300)/1000'
'define v=vint(pres.4,se.3*a.1,300)/1000'

'set gxout shaded'
'set lon 65 160'
'set lat -15 45'
'set strsiz 0.28'

'set gxout vector'
'set arrowhead 0.1'
'define qu=ave(u,t=1,t=45)'
'define qv=ave(v,t=1,t=45)'
'd mag(qu,qv)'
*'cbarn'
'd qu;qv'
'print'
'disable print'

修改后，特别注意算整层水汽通量，要有地表气压参与。ctl书写变量时的个数要与z层数相同。还可以再除以一个9.8

洗雨

这个图片多发了一次，彩图是自己画的，黑白图是在论文上找的，论文图所用资料为NCEP/NCAR 1957-2001年日平均的再分析资料包括水平u纬向风，V经向风,比湿q，精度2.5*2.5的

river

你看看你的gs里面高度有循环吗？

quezhiping

图有的地方没有值是因为你“'define q=hdivg(u,v)*100000'”用到了求散度的函数，需要将定义范围扩大一点然后画图范围稍小一点就可以了

洗雨

river 发表于 2015-4-15 08:17
你看看你的gs里面高度有循环吗？

计算水汽通量不是从第一层积分到第八层了吗！我在读取数据时就有高度循环了。后面不需要吧？

river

洗雨 发表于 2015-4-15 13:06
计算水汽通量不是从第一层积分到第八层了吗！我在读取数据时就有高度循环了。后面不需要吧？

我是说你比湿那个gs

洗雨

quezhiping 发表于 2015-4-15 09:30
图有的地方没有值是因为你“'define q=hdivg(u,v)*100000'”用到了求散度的函数，需要将定义范围扩大一点然 ...

谢谢，这个倒是改进了

洗雨

river 发表于 2015-4-15 16:29
我是说你比湿那个gs

没注意看，后面几次修改所以不小心把那个while语句去掉了。图倒是画出来的是夏季平均的，但是和人家用日平均画出来的不同

洗雨

这个问题又不知道怎么解决了

乖乖孙迪

学习啦！！！