计算水汽通量用到垂直积分遇到了问题

gercibee

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/popRemap.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
    sp1 = addfile("rhum.2013.nc","r")
    v1 = addfile("vwnd.2013.nc","r")
    u1 = addfile("uwnd.2013.nc","r")
    p1 = addfile("pres.sfc.2013.nc","r")

    v3 = short2flt(v1->vwnd(:,0:7,:,:))
    u3 = short2flt(u1->uwnd(:,0:7,:,:))
    sp3 = short2flt(sp1->rhum(:,0:7,:,:))
    p3 = short2flt(p1->pres)

    time = v1->time
;-------------计算最大地面气压-----------------------------
    p4 = p3/100
    pmax = max(p4)
;----------------------------------------------------------
    uq1 = u3*sp3
    vq1 = v3*sp3


    lev = (/1000,925,850,700,600,500,400,300/)
    lat = fspan(-90,90,73)
    lon = fspan(0,357.5,144)

    time!0 = "time"
    lon!0 = "lon"
    lon@units = "degrees_east"
    lat!0 = "lat"
    lat@units = "degrees_north"
    lev!0 = "lev"
    lev@unite = "hpa"

    uq1!0 = "time"
    uq1!1 = "lev"
    uq1!2 = "lat"
    uq1!3 = "lon"
    uq1&lon = lon
    uq1&lat = lat
    uq1&lev = lev
    uq1&time = time

    vq1!0 = "time"
    vq1!1 = "lev"
    vq1!2 = "lat"
    vq1!3 = "lon"

    p4!0 = "time"
    p4!1 = "lat"
    p4!2 = "lon"
    p4&lon = lon
    p4&lat = lat
    p4&time = time

    uq2 = uq1(time|:,lat|:,lon|:,lev|:)
    vq2 = vq1(time|:,lat|:,lon|:,lev|:)

    p = (/1000.,925.,850.,700.,600.,500.,400.,300./)
    linlog = 1  
    pbot = 1000.
    ptop = 300.
    uq3 = (vibeta(p,uq2,linlog,p4,pbot,ptop))/9.8
    vq3 = (vibeta(p,vq2,linlog,p4,pbot,ptop))/9.8

    uq3!0 = "time"
    uq3!1 = "lat"
    uq3!2 = "lon"
    uq3&lon = lon
    uq3&lat = lat
    uq3&time = time

    vq3!0 = "time"
    vq3!1 = "lat"
    vq3!2 = "lon"
    vq3&lon = lon
    vq3&lat = lat
    vq3&time = time

    printVarSummary(uq3)
end
算到这一步的时候总出现warning:vibeta: there must be at least three levels with data above the surface这样的提示，uq3没有值，请大家指正一下代码有哪里出现错误了吗？

longlivehj

水汽通量计算要用比湿，看代码楼主好像用的是相对湿度，这个要注意一下。
警告是由于地形以上层数不足3层造成的，可以考虑人为添加几层，数值为0，试试看。

xiaohezhou

来学习学习

排山倒海

我按照你的脚本，用我自己的数据算了一下，也是这个结果，请问你的解决了吗？

gercibee

排山倒海 发表于 2015-3-10 12:07
我按照你的脚本，用我自己的数据算了一下，也是这个结果，请问你的解决了吗？

我后来用的是ERA的数据，就没有这个问题，结果还不错。如果非要用ncep数据，我知道的是有人把缺测值改成了别的数。这样就可以运算了。

排山倒海

gercibee 发表于 2015-3-10 19:48
我后来用的是ERA的数据，就没有这个问题，结果还不错。如果非要用ncep数据，我知道的是有人把缺测值改成 ...

好的，我再去试试，谢啦~

奚子惠

gercibee 发表于 2015-3-10 19:48
我后来用的是ERA的数据，就没有这个问题，结果还不错。如果非要用ncep数据，我知道的是有人把缺测值改成 ...

楼主你好~我也遇到了这样的情况，我的资料是ncep月平均资料，但是我记得这个数据里面是没有缺省值的吧~所以“有人把缺测值改成别的数”，这个咋理解呢

奚子惠

gercibee 发表于 2015-3-10 19:48
我后来用的是ERA的数据，就没有这个问题，结果还不错。如果非要用ncep数据，我知道的是有人把缺测值改成 ...

其实我还有一个想法就是可不可以把比湿的资料（原来是8层）自己补充扩展到跟u，v原来的层数相同（17层），然后 p       = (/1000.,925.,850.,700.,600.,500.,400.,300.,250.,200.,150.,100.,70.,50.,30.,20.,10./) ，不过我还没有试，不清楚可不可以，因为我在算其他物理量的垂直积分时 在层次上就是原本资料的层数没有改动 ，程序就没有报错

sunyue19901203

楼主，请问一下积分的时候一定要把地面气压的单位换算成hPa么？还有你最后计算出来的水汽通量是多少量级呢？单位是什么？

timevigor

请问计算地面最大气压要除以100？，不可以直接用p1的数据吗