grads计算水汽通量求助！！！！（已解决，并更新）

wet510

问题已解决，不清楚的同学可见14楼，那个是比较简单的求水汽通量的grads脚本。另外，求整层水汽通量一般采用积分的方式，而不是循环累加---------------------------------------------------------华丽的分割线---------------------------------------------------------------------
重新编辑了下帖子，关于水汽通量的单位问题，推荐去下面的帖子看看：
http://bbs.06climate.com/forum.p ... 141&page=4#lastpost
另外，尤其注意是grads中的vint函数的应用。

wet510

hustlh 发表于 2012-7-25 15:27

                               
登录/注册后可看大图

那我只计算整层（1000hpa-300hpa），是不是用Z循环就可以了？？？

其实如果算整层的话不用循环也可以，如果你要求确定时间的一个区域的整层水汽通量，可以如下设置维度，比如：

1. 'sdfopen E:\TDDOWNLOAD\nc_data\hgt.1998.nc'
   
2. 'sdfopen E:\TDDOWNLOAD\nc_data\shum.1998.nc'
   
3. 'sdfopen E:\TDDOWNLOAD\nc_data\uwnd.1998.nc'
   
4. 'sdfopen E:\TDDOWNLOAD\nc_data\vwnd.1998.nc'
   
5. 'sdfopen E:\TDDOWNLOAD\nc_data\pres.sfc.1998.nc'
   
6. 'set lat -10 80'
   
7. 'set lon 60 160'
   
8. 'set z 1 8'
   
9. 'set time 12z01jun1998 12z31aug1998'
   
10. ************vapar flux***************
    
11. 'define uq=shum.2*uwnd.3'
    
12. 'define vq=shum.2*vwnd.4'
    
13. 'define uqs=vint(pres.5(lev=0),uq,300)'
    
14. 'define vqs=vint(pres.5(lev=0),vq,300)'

复制代码

计算完，vqs和uqs实际上只有3层了，即经度、维度层、时间层，其高度层由于积分都是一层，自己可以试验下，因此，为了画图或输出数据的话，还要用
'set z 1'
'set t 1'
时间具体可以自己改，高度改掉，可以观察结果是否有变化。
然后接你要绘图或输出结果的语句

禁断2次元

研究研究啊

wet510

顶一个，高手们快来

maplesnow82

在计算多层的uq,vq时，不需要使用z循环。
set lon 115 145
set lat 35 60
set z 1 8
define es=...
define qs=...
define qv=...
define uq=..
define vq=...
按照你的外循环，出来的值肯定只有最后一个值啊。 就像 fortran的算法一样啊，因为你的uq,vq没有被定义成数组。 如果需要的话， 我建议可以把算出来的值写到另外一个文件里进行检查：
下面贴一个我以前写的算假相当位温的gs文件，把数据输到一个二进制文件里：
'reinit'
'open f:\data\wgrib.big\7-1-00.ctl'
'set gxout fwrite'
'set  fwrite  f:\data\wgrib.big\gs2005\thes\thes07.dat'

i=90
while(i<=116)
'set t  '   i
zz=1
while(zz<=26)
'set lat 0  90'
'set lon 50 180'
'set z 'zz
'define tt=TMPprs'
'define es=6.11*exp(17.26*(tt-273.16)/(tt-35.86))'
'define qs=0.622*es/(lev-0.378*es)'
'define qq=RHprs/100.*qs'

'define prs=lev*100'
'define el=log((qq/0.62197)*(prs/100.)/(1.+(qq/0.62197)))'
'define td=273.15+(243.5*el-440.8)/(19.48-el)'
'define tc=tt-1.21*(tt-td)'
'define e=qq/(qq+0.622)*prs'
'define cpm=1004*(1+0.86*qq)'
'define rm=287.04*(1+0.608*qq)'
'define pc=prs*pow(tc/tt,cpm/rm)'
'define st=tc*pow(100000/(pc-e),287.04/1004)'
'define lvc=2500.8*1000-2274.4*(tt-273.15)'
'define qse=st*exp(qq*lvc/1004/tc)'
'undefine el'
'undefine prs'
'undefine td'
'undefine tc'
'undefine es'
'undefine qs'
'undefine qq'
'undefine e'
'undefine cpm'
'undefine rm'
'undefine pc'
'undefine st'
'undefine lvc'
'd qse'
zz=zz+1
endwhile

希望有所帮助！

wet510

maplesnow82 发表于 2012-3-14 22:18

                               
登录/注册后可看大图

在计算多层的uq,vq时，不需要使用z循环。
set lon 115 145
set lat 35 60

早上刚爬起来就看到高手回帖，激动啊。我马上试试看先行谢过，不懂得地方还要请教

topmad

wet510 发表于 2012-3-15 09:21

                               
登录/注册后可看大图

早上刚爬起来就看到高手回帖，激动啊。我马上试试看先行谢过，不懂得地方还要请教

同学你好谦虚....

wet510

果然如此，不过在显示mag(uqs,vqs)前要加个环境设定，set z 1，设置高度为一层，才能画出图，但mag(uqs,vqs)每一层的值都一样，所以设置哪层无所谓。

maplesnow82

wet510 发表于 2012-3-15 09:54

                               
登录/注册后可看大图

果然如此，不过在显示mag(uqs,vqs)前要加个环境设定，set z 1，设置高度为一层，才能画出图，但mag(uqs,vqs ...

Congratulations! Problem is solved!

wutongyeyu

顶一个

hustlh

呵呵，我也遇到这个问题了，幸好有高手指点啊，真是太感谢了！

我再试试。

buddy83664

测试 重新登录情况