求助帖：利用GRADS函数hcurl()求水平涡度值

sweet33

grads新手在线等~~~各路大神，求救~~~
利用hcurl()函数求解水平x方向的涡度值，需要计算的范围：12.5~25N，80~90E，1000~600hpa的空间范围，原意是计算该范围66年4月的水平涡度值，以此来代表环流圈的强度。但是在计算过程发现，1000hpa边界层的涡度值一并计算出来了（后来才改为925hpa，但是计算结果没什么区别），这似乎不合理，缺乏真实性。结果也是有数值变化的，附件有图。用这个函数是否能求得水平涡度有待考证，在这里贴上gs。各路大神，请问有没有谁做过水平涡度的，请指教一二，小女子在此拜谢啦~~~
file:///C:/Users/Administrator/AppData/Roaming/Tencent/Users/634463954/QQ/WinTemp/RichOle/($Z1)L%7DH9PM~9O~U0D_1XVC.png

'reinit'
'sdfopen F:\data\NCEP\vwnd.mon.mean.nc'
'sdfopen F:\data\NCEP\omega.mon.mean.nc'
'enable print F:\reviseFig\shuipingwodu\w.gmf'

t1=4
while(t1<=792)
'set t 't1''
say t1
'set lat 10 27.5'
'set lon 77.5 92.5'
'set lev 925 600'
'set gxout line'
'set grads off'
'set grid off'
'set zlog on'
'define vort=hcurl(vwnd,omega.2)*1e6'
'define vu=aave(vort,lon=80,lon=90,lat=12.5,lat=25)'
'define zw=ave(vu,lev=925,lev=600)'
'set lat 13'
*'set lon 80 90'
*'set lev 600'
'd zw'
'print'
'c'
t1=t1+12
endwhile
'disable print'
'reinit'
;

river

你看一下grads求涡度的原理，去掉1000hPa应该有很大变化吗，我也不清楚，很久没弄过这些了。
define dv = cdiff(v,x)
define dx = cdiff(lon,x)*3.1416/180
define du = cdiff(u*cos(lat*3.1416/180),y)
define dy = cdiff(lat,y)*3.1416/180
display (dv/dx-du/dy)/(6.37e6*cos(lat*3.1416/180))
The above example assumes an X-Y varying dimension environment. Note that the intrinsic variables lat and lon give results in degrees and must be converted to radians in the calaculation. Also note the radius of the earth is assumed to be 6.37e6 meters thus the U and V winds are assumed to have units of m/sec.

sweet33

river 发表于 2015-3-24 13:24
你看一下grads求涡度的原理，去掉1000hPa应该有很大变化吗，我也不清楚，很久没弄过这些了。
define dv =  ...

嗯嗯，O(∩_∩)O谢谢~~~弄了几天，终于弄清楚了。由量纲分析将涡度方程简化为ξ=- 9v / 9z，计算时直接去掉第一项。因为我求的是z方向的涡度，不能直接用cdiff来计算。根据4月的风场v资料，根据中央差分的原理公式用Grads手动计算偏导，在变微分为差分时，对于925hPa——600hPa各层，计算ξ时取向上差分与向下差分的平均（中央差分），1000hPa按照向上差分进行计算；对80～90°E，12.5～25°N，1000～600hPa范围内各点计算ξ，得到一个类似温度场的涡度空间场。依次计算1948-2013年逐年区域平均，得到66年的时间序列，再对其标准化，得到标准化强度指数，即为环流强度指数。

river

sweet33 发表于 2015-3-27 10:17
嗯嗯，O(∩_∩)O谢谢~~~弄了几天，终于弄清楚了。由量纲分析将涡度方程简化为ξ=- 9v / 9z，计算时直接去 ...

好复杂的说，楼主弄好了可以来论坛分享一下

sweet33

river 发表于 2015-3-27 12:16
好复杂的说，楼主弄好了可以来论坛分享一下

分享下我求出的x方向的水平涡度：
gs：
'reinit'
'open F:\reviseFig\shuipingwodu\v.wodu\April.v.66.ctl'
'set fwrite F:\reviseFig\shuipingwodu\v.wodu\v.wodu.grd'
'set gxout fwrite'

'set lon 80 90'
'set lat 12.5 25'
i=1
while(i<=66)
'set t 'i''

'set z 1'
'define dv1=vwnd(z=2)-vwnd(z=1)'
'define dp1=100*75'
'define dvdp1=-(dv1/dp1)'
'd dvdp1'

'set z 2'
'define dv2=vwnd(z=3)-vwnd(z=1)'
'define dp2=100*150'
'define dvdp2=-(dv2/dp2)'
'd dvdp2'

'set z 3'
'define dv3=vwnd(z=4)-vwnd(z=2)'
'define dp3=100*125'
'define dvdp3=-(dv3/dp3)'
'd dvdp3'

'set z 4'
'define dv4=vwnd(z=5)-vwnd(z=3)'
'define dp4=100*150'
'define dvdp4=-(dv4/dp4)'
'd dvdp4'

'set z 5'
'define dv5=vwnd(z=6)-vwnd(z=4)'
'define dp5=100*200'
'define dvdp5=-(dv5/dp5)'
'd dvdp5'

i=i+1
endwhile
'disable fwrite'
'reinit'
;

此处的数据资料为ncep月资料，经过处理后只留下4月的资料。

river

sweet33 发表于 2015-3-30 12:43
分享下我求出的x方向的水平涡度：
gs：
'reinit'

感谢分享

xxxzsn

sweet33 发表于 2015-3-30 12:43
分享下我求出的x方向的水平涡度：
gs：
'reinit'

这样写出来的数据，后面的不会把前面的覆盖掉吗？

sweet33

xxxzsn 发表于 2015-4-8 15:51
这样写出来的数据，后面的不会把前面的覆盖掉吗？

不会的哈，出来的数据是80～90°E，12.5～25°N，1000～600hPa范围的涡度空间场。

木易

sweet33 发表于 2015-3-30 12:43
分享下我求出的x方向的水平涡度：
gs：
'reinit'

楼主我想请教您问一下每个层次中的dp=100*125或者100*200表示的是什么意思呢？？我看每个层次dp数值都不一样

SHIXIAOLAN

谢谢楼主分享{:eb301:}