求助：整层水汽通量散度积分问题

红色的兔子

最近用grads计算了整层水汽通量散度积分
我查过相关帖子，发现在计算整层水汽通量散度时有个问题
一般步骤是
1.首先计算各层的水汽通量qu和qv
但是第二步计算整层水汽通量散度的时候有两种方法：
有人先对qu、qv分别积分 得quall和qvall
然后再hdivg（quall,qvall）得到整层的水汽通量散度

还有人先求出各层水汽通量散度q=hdivg（qu,qv）
然后再对q进行积分 得到整层的水汽通量散度

在 http://bbs.06climate.com/forum.p ... 21104&extra=&page=1 帖子中，回帖的人多数认为两种方法都可以。
但是最近我计算中发现这两种方法计算结果有明显差异，我现在将两个gs和两个图都放上来，请问两个结果的正确性。

详见下网址
http://pan.baidu.com/s/1eS7DN8Y

两个gs：
第一个：
'reinit'
'open H:\js\shum.ctl'
'open H:\js\uave.ctl'
'open H:\js\vave.ctl'
'open H:\js\pres.ctl'

****************East area**************************
'set x 1 144'
'set y 1 73'
'set z 1 8'
'define uq=ave(uwnd.2*q.1/9.8,t=1,t=60)'
'define vq=ave(vwnd.3*q.1/9.8,t=1,t=60)'
'define sandu=hdivg(uq,vq)'
'set z 1'
'define ps=ave(0.01*pres.4,t=1,t=60)'
'define zhcsandu=vint(ps,sandu,300)'

'set grads off';'set grid off';;'set map 1 1 3'
'set parea 0.7 10.5 0.8 8'
'set lat 15 50'
'set lon 95 140'
'set xlopts 1 4 0.20';'set ylopts 1 4 0.20';
'set xlint 10';'set ylint 10 '
'set grads off'

'set gxout contour'
'set ccols 1'
'set csmooth on'
'set cint 2e-6'
*'set clab off'
'd zhcsandu'

'printim H:\js\huabei1951-2010zhcsandu.png white'
'printim H:\js\huabei1951-2010zhcsandu.eps white'
;

第1个gs对应的图，先计算每层的散度，再进行积分

第二个：
'reinit'
'open H:\js\shum.ctl'
'open H:\js\uave.ctl'
'open H:\js\vave.ctl'
'open H:\js\pres.ctl'

****************East area**************************
'set x 1 144'
'set y 1 73'
'set z 1 8'
'define uq=ave(uwnd.2*q.1/9.8,t=1,t=60)'
'define vq=ave(vwnd.3*q.1/9.8,t=1,t=60)'
'set z 1'
'define ps=ave(0.01*pres.4,t=1,t=60)'
'define uqx=vint(ps,uq,300)'
'define vqx=vint(ps,vq,300)'
'define sandu=hdivg(uqx,vqx)'

'set grads off';'set grid off';;'set map 1 1 3'
'set parea 0.7 10.5 0.8 8'
'set lat 15 50'
'set lon 95 140'
'set xlopts 1 4 0.20';'set ylopts 1 4 0.20';
'set xlint 10';'set ylint 10 '
'set grads off'

'set gxout contour'
'set ccols 1'
'set cint 2e-6'
'set csmooth on'
*'set clab off'
'd sandu'

'printim H:\js\huabei1951-2010zhcsandu2.png white'
'printim H:\js\huabei1951-2010zhcsandu2.eps white'
;

第2个gs对应的图，先计算了整层水汽通量，再计算散度

帅气小葛

我看到的所有文献中，垂直积分的符号是在水汽散度的外面，这就意味着应该先求各层散度，然后再垂直积分，我个人的理解

男紫汗

我之前在版上看到有人说散度是周围几个格点求出来的，先求散还是先积分在边界上会不同。对uq作平均也可以？我只知道u和q分别作平均再相乘，唉，这里的问题搞不清楚，你uq平均后出来的水汽通量形势和欧洲中心或者文献中一样么？

红色的兔子

男紫汗 发表于 2016-7-1 12:13
我之前在版上看到有人说散度是周围几个格点求出来的，先求散还是先积分在边界上会不同。对uq作平均也可以？ ...

uq的平均只是时间上的平均，会对结果产生影响吗？

nunu18

楼主好细心，这个要估计数学很牛的人才能解答吧，我看着两种方法差不多，不知道为啥结果差那么多，帮你顶吧，希望帖子不要沉下来，有个高手出来解惑

river

没太注意过，楼主看下官方说明吧
vint

vint(psexpr,expr,top)

This function performs a mass-weighted vertical integral in mb pressure coordinates. The three arguments to vint are:

psexpr   a GrADS expression for the surface pressure, in mb, which bounds the integral on the bottom.
expr     a GrADS expression representing the quantity to be integrated.
top      the bounding top pressure, in mb. This value must be a constant and cannot be provided as an expression.

The calculation is a sum of the mass-weighted layers:

f/g * sum(expr * Delta(level))
The bounds of the integration are the surface pressure and the indicated top value. The scale factors are f=100 and g=9.8. The summation is done for each layer present that is between the bounds. The layers are determined by the Z levels of the default file. Each layer is considered to be from the midpoints between the levels actually present, and is assumed to have the same value throughout the layer, namely the value of the gridpoint at the middle of the layer.

Usage Notes

The summation is done using the Z levels from the default file, so it is important that the default file have the same Z dimension coordinates as expr.
Data levels below and above the bounds of the summation are ignored.
The Z dimension in world-coordinate units is assumed to be pressure values given in millibars (mb). The units of g are such that when the expression integrated is specific humidity (q) in units of g/g, the result is kg of water per square meter, or precipitable water in mm.
It is usually a good idea to make the top pressure value to be at the top of a layer, which is midway between grid points. For example, if the default file (and the data) have pressure levels of ...,500,400,300,250,... then a good value for top might be 275, the value at the top of the layer that extends from 350 to 275 mb.
The vint function operates only in an X-Y varying dimension environment.
Be sure the units of the surface pressure (argument 1) are in millibars (mb).
Examples

1. This expression will integrate specific humidity to obtain precipitable water, in mm:

vint(ps,q,275)
2. This is an artificial example that demonstrates a vertical integration from a fixed lower bound of 1000mb to the top of the atmosphere, and integrating a field of all 1's. This gives an answer of 10204.1 (or 100000/9.8) which is the mass of air (in kg) of a 1 meter squared column when the surface pressure is 1 bar and the accelleration due to gravity is assumed to be exactly 9.8m/sec**2 over the entire column.

vint(const(ps,1000),const(t,1),0)

nunu18

这个帖子终究没有人解答么

大扫把

小半年过去了。这么重要的问题还没有答案。我也想明白

may6117

几个版本，也是有点晕，继续研究

Alexander

{:5_235:}{:5_275:}{:5_275:}为什么啊

缥缈

理论上第二种是对的吧，一孔之见