[已解决]GraDS画理查森数Ri

hdswy

想利用WRF的模拟输出资料画理查森数Ri，但是画出的等值线图1几乎没有显示出多少，想请教大家什么地方出了问题。
图2 为在格点上显示出的值（'d skip(Ri,2) '）


图1      

grid

图2


gs文件如下：
'reinit'
' open d:/arw/testd01.ctl '
'set mpdset cnriver'
'set grid off'
'set lon 102 114 '
'set lat 31 41 '
'set t 1 '

' set grads off '
'set lev 850 500'
'define thse=theta*exp(2.5e6*QVAPOR/1004/tk) '
'define rou=pressure*1.e2/(9.8*height*1.e3)'

'set lev 850'
'define ro=(rou(lev=850)+rou(lev=500))/2 '
'define dth=thse(lev=850)-thse(lev=500) '
'define the=(thse(lev=850)+thse(lev=500))/2 '
'define du=u(lev=850)-u(lev=500)'
'define dv=v(lev=850)-v(lev=500)'
'define Ri=-dth*3.5e5/(the*(du*du+dv*dv)*ro)'
'd Ri '

gs中第二段资料中的变量为
theta         19  0  Potential Temperature (K)
QVAPOR        19  0  Water vapor mixing ratio (kg kg-1)
tk            19  0  Temperature (K)
height        19  0  Model height (km)
pressure      19  0  Model pressure (hPa)
U             19  0  x-wind component (m s-1)
V             19  0  y-wind component (m s-1)

定义的thse为假相当位温，rou是密度，ro是850到500hPa的平均密度，用的公式是Ri =-g/θse*(əθse/əz)/(əVh/əz)2,差分形式Ri=-Δθse / θse平 * Δp / （ρ*(ΔVh)2）



------------------------------------------------------------------------------

经过别人提醒，gs中第三段 Δp 是350hPa，应该是3.5e4，另外用maskout函数隐去了大于30的值。另外，由于850hPa上缺值太多，所以上下层分别取500hPa和700hPa，作出图如图3。
图3

mofangbao

楼主先把网格点的值标志网格上看看，看能不能看出点问题，比如有个特别大的值或者特别小的值影响了等值线插值时的计算呢

hdswy

mofangbao 发表于 2012-3-25 11:27

                               
登录/注册后可看大图

楼主先把网格点的值标志网格上看看，看能不能看出点问题，比如有个特别大的值或者特别小的值影响了等值线插 ...

嗯 图2右下角附近的数太大。一般Ri算出来也就一位数到两位数，不知道我哪里算错了

hdswy

mofangbao 发表于 2012-3-25 11:27

                               
登录/注册后可看大图

楼主先把网格点的值标志网格上看看，看能不能看出点问题，比如有个特别大的值或者特别小的值影响了等值线插 ...

额……其实是我自己把量级算大了一位……谢谢诶

HZK

空气密度和假相当位温是计算出来的么？求计算公式。。

HZK

能不能用这个公式，用z坐标算，高度取等压面上的位势高度乘以g，这样可以不用算位温和密度

hdswy

HZK 发表于 2013-3-1 20:49

                               
登录/注册后可看大图

空气密度和假相当位温是计算出来的么？求计算公式。。

密度就是根据p=rou*R*T算的，数据结果中有位温，所以假相当位温thse=theta*exp(2.5e6*Q/1004/tk)

徐徐向前

能不能换个色标呀

涛声依旧否

顶一个先！

nzy718d

请教楼主：是哪个地方的量级错了，能指导下吗