用NCL作整层水汽通量的问题

BDRUSH

下载了ncep上u,v,shum,pres的逐日资料，然后提取某三天进行整层水汽通量的合成计算。具体脚本如下：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/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/wind_rose.ncl"
begin                    
;**************************************************************************************
;                         读U风
;**************************************************************************************
    u_file ="/cygdrive/H/Data/grd/uwind2013/uwind2013_1000to300hP.grd"                           
    a       = fbindirread(u_file,0,(/365,8,73,144/),"float")             ;0代表是第一行.  8是层次
   u        =  a(183,:,:,:)
   u1      =  a(228,:,:,:)
   u2      = a(240,:,:,:)
;**************************************************************************************
;                         读v风
;**************************************************************************************
        v_file = "/cygdrive/H/Data/grd/vwind2013/vwind2013_1000to300hP.grd"                           
    a1     = fbindirread(v_file,0,(/365,8,73,144/),"float")             ;0代表是第一行
        v      = a1(183,:,:,:)
        v1     = a1(228,:,:,:)
        v2     = a1(240,:,:,:)
;**************************************************************************************
;                         读shum比湿
;**************************************************************************************
        shum_file ="/cygdrive/H/Data/grd/shum2013/shum2013_1000to300hP.grd"                           
    a2        = fbindirread(shum_file,0,(/365,8,73,144/),"float")
        q         = a2(183,:,:,:)
        q1        = a2(228,:,:,:)
        q2        = a2(240,:,:,:)
;**************************************************************************************
;                                                  读取地面气压
;**************************************************************************************
        a3            = addfile("/cygdrive/H/Data/ncep/pres.sfc/pres.sfc.2013.nc","r")
        ps1           = short2flt(a3->pres(183,::-1,:))                              
        ps2           = short2flt(a3->pres(228,::-1,:))
        ps3           = short2flt(a3->pres(240,::-1,:))
;**************************************************************************************
;                                                  计算qu = q(比湿)乘以u(u风)
;**************************************************************************************        
        qu1 = new((/8,73,144/),float)
        qv1 = new((/8,73,144/),float)
        qu2 = new((/8,73,144/),float)
        qv2 = new((/8,73,144/),float)
        qu3 = new((/8,73,144/),float)
        qv3 = new((/8,73,144/),float)
        
        do iz=0,7
                do i=0,72
                        do j=0,143
                                qu1(iz,i,j) = q(iz,i,j)*u(iz,i,j)
                                qv1(iz,i,j) = q(iz,i,j)*v(iz,i,j)
                        end do
                end do
        end do
        
        do iz=0,7
                do i=0,72
                        do j=0,143
                                qu2(iz,i,j) = q1(iz,i,j)*u1(iz,i,j)
                                qv2(iz,i,j) = q1(iz,i,j)*v1(iz,i,j)
                        end do
                end do
        end do
        
        do iz=0,7
                do i=0,72
                        do j=0,143
                                qu3(iz,i,j) = q2(iz,i,j)*u2(iz,i,j)
                                qv3(iz,i,j) = q2(iz,i,j)*v2(iz,i,j)
                        end do
                end do
        end do
;**************************************************************************************
;                                                  给qu1,qv1赋予属性
;**************************************************************************************               
        lev =(/1000,925,850,700,600,500,400,300/)
    lat=fspan(-90,90,73)
    lon=fspan(0,357.5,144)
        
        lon!0                 = "lon"
        lon@units         = "degrees_east"
        lat!0                 = "lat"
        lat@units         = "degrees_north"
        lev!0                = "lev"
        lev@units        = "hPa"
        
        qu1!0         = "lev"
        qu1!1         = "lat"
        qu1!2         = "lon"
        qu1&lon = lon
        qu1&lat = lat
        qu1&lev = lev
        
        qu2!0         = "lev"
        qu2!1         = "lat"
        qu2!2         = "lon"
        qu2&lon = lon
        qu2&lat = lat
        qu2&lev = lev

        qu3!0         = "lev"
        qu3!1         = "lat"
        qu3!2         = "lon"
        qu3&lon = lon
        qu3&lat = lat
        qu3&lev = lev
        
        qv1!0         = "lev"
        qv1!1         = "lat"
        qv1!2         = "lon"
        qv1&lon = lon
        qv1&lat = lat
        qv1&lev = lev
        
        qv2!0         = "lev"
        qv2!1         = "lat"
        qv2!2         = "lon"
        qv2&lon = lon
        qv2&lat = lat
        qv2&lev = lev
        
        qv3!0         = "lev"
        qv3!1         = "lat"
        qv3!2         = "lon"
        qv3&lon = lon
        qv3&lat = lat
        qv3&lev = lev
;**************************************************************************************
        qu11 = qu1(lat|:, lon|:, lev|:)
        delete(qu1)
        qv11 = qv1(lat|:, lon|:, lev|:)
        delete(qv1)
        
        qu22 = qu2(lat|:, lon|:, lev|:)
        delete(qu2)
        qv22 = qv2(lat|:, lon|:, lev|:)
        delete(qv2)
        
        qu33 = qu3(lat|:, lon|:, lev|:)
        delete(qu3)
        qv33 = qv3(lat|:, lon|:, lev|:)
        delete(qv3)
;**************************************************************************************        
        p = (/1000,925,850,700,600,500,400,300/)
        linlog = 2                                                                                                         ; 1为线性积分  2为对数积分
        data11 = (vibeta(p,qu11,linlog,ps1,1000,300))/9.8           ; 返回的数值是2维的(lat,lon)
        data12 = (vibeta(p,qv11,linlog,ps1,1000,300))/9.8
        data21 = (vibeta(p,qu22,linlog,ps2,1000,300))/9.8           ; 返回的数值是2维的(lat,lon)
        data22 = (vibeta(p,qv22,linlog,ps2,1000,300))/9.8
        data31 = (vibeta(p,qu33,linlog,ps3,1000,300))/9.8           ; 返回的数值是2维的(lat,lon)
        data32 = (vibeta(p,qv33,linlog,ps3,1000,300))/9.8

        qu  = new((/73,144/),float)
        qv  = new((/73,144/),float)
        mm  = new((/73,144/),float)
        mm1 = new((/73,144/),float)
        mm2 = new((/73,144/),float)
        mm3 = new((/73,144/),float)

        mm1 = sqrt(data11^2+data12^2)
        mm2 = sqrt(data21^2+data22^2)
        mm3 = sqrt(data31^2+data32^2)
        
        do i=0,72
                do j=0,143
                        mm(i,j) = (mm1(i,j)+mm2(i,j)+mm3(i,j))/3
                end do
        end do

        do i=0,72
                do j=0,143
                        qu(i,j) = (data11(i,j)+data21(i,j)+data31(i,j))/3
                end do
        end do

        do i=0,72
                do j=0,143
                        qv(i,j) = (data12(i,j)+data22(i,j)+data32(i,j))/3
                end do
        end do        

        mm!0 = "lat"
        mm!1 = "lon"
        mm&lat = lat
        mm&lon = lon
        lon!0 = "lon"
        lon@units = "degrees_east"
        lat!0 = "lat"
        lat@units = "degrees_north"        
        qu!0 = "lat"
        qu!1 = "lon"
        qu&lat = lat
        qu&lon = lon
        lon!0 = "lon"
        lon@units = "degrees_east"
        lat!0 = "lat"
        lat@units = "degrees_north"
        qv!0 = "lat"
        qv!1 = "lon"
        qv&lat = lat
        qv&lon = lon
        lon!0 = "lon"
        lon@units = "degrees_east"
        lat!0 = "lat"
        lat@untis = "degrees_north"
;************************************************************************************************
        wks =gsn_open_wks("pdf","111")
        gsn_define_colormap(wks,"BlAqGrYeOrReVi200")

        map_res                                  = True
        map_res@gsnFrame      = False
        map_res@gsnDraw       = False
        map_res@mpDataSetName         = "Earth..4"   ; This new database contains
                                           ; divisions for other countries.
        map_res@mpDataBaseVersion     = "MediumRes"  ; Medium resolution database
        map_res@mpFillOn=False
        map_res@mpGridMaskMode           = "MaskLand"  ; Don't draw grid over land.
        map_res@mpGeophysicalLineColor =   "black"
        map_res@mpGeophysicalLineThicknessF =  1.0
        map_res@mpOutlineOn           = True         ; Turn on map outlines

        map_res@mpLimitMode="LatLon"
        map_res@mpCenterLonF=120   
        map_res@mpMinLonF=60        
        map_res@mpMaxLonF=180         
        map_res@mpMinLatF=0
        map_res@mpMaxLatF=80
        map_res@mpPerimOn=True
        
        map = gsn_csm_map(wks,map_res)
;***************************************************************************
;
;***************************************************************************        
        vcres                                         = True                    ; plot mods desired
        vcres@gsnAddCyclic                            = False
        vcres@gsnDraw                                 = False                   ; don't draw yet
        vcres@gsnFrame                                = False                   ; don't advance frame yet

        vcres@pmTickMarkDisplayMode         = "Always"        
        vcres@vcMinFracLengthF                        = 0.33        ;设置了矢量对于参考矢量的长度的最小数量级(default=0)0.33
        vcres@vcRefMagnitudeF                        = 10          ;20     风速矢量长度参考
        vcres@vcRefLengthF                                = 0.022           ;0.045  0.02
        vcres@vcGlyphStyle                                 = "CurlyVector" ; turn on curly vectors
        vcres@vcMonoLineArrowColor                 = False
        vcres@vcMinDistanceF                     = 0.013
        vcres@vcLineArrowHeadMaxSizeF         = 0.008   ;default = 0.012

        vcres@lbLabelBarOn                                 = False



        vcres@vcMonoFillArrowFillColor  = True
        vcres@vcLineArrowThicknessF                  = 1.0
        vcres@vcRefAnnoOn                                  =False

        vcres@gsnLeftString                                 =""
        vcres@gsnRightString                                 =""
        vcres@tiMainString                                         =""
        qu = qu*100
        qv = qv*100
        vector = gsn_csm_vector(wks,qu,qv,vcres)  ; create plot
;***************************************************************************        
;
;***************************************************************************        
        res                                         = True                    ; plot mods desired
        res@gsnAddCyclic                            = False
        res@gsnDraw                                 = False                   ; don't draw yet
        res@gsnFrame                                = False                   ; don't advance frame yet
        
        res@cnFillOn                                = True                    ; turn on color
        res@gsnSpreadColors                         = True                    ; use full colormap
        res@gsnSpreadColorStart                 = 0
        res@gsnSpreadColorEnd                   = 199
        res@cnLinesOn                                      = False                   ; turn off contour lines
        res@cnLineLabelsOn                          = False                   ; tuen off line labels
        res@cnInfoLabelOn                                 = False

        res@cnLevelSelectionMode                = "ExplicitLevels"
        res@cnLevels                                        = (/2,4,6,8,10,12,14,16,18,20/)                 ;(/14,18,22,26,30,34,38,42,46,50,52/)

        res@gsnLeftString                                 =""
        res@gsnRightString                                 =""
        res@tiMainString                                 =""

        res@lbLabelBarOn                                 = True

        res@lbOrientation                        = "Vertical"         ; vertical label bars
        res@lbAutoManage                            = True
        res@lbLabelFontHeightF                  = 0.009              ; make labels smaller

        mm = smth9(mm,0.5,0.25,False)
        mm = mm*100
        contour = gsn_csm_contour(wks,mm,res)
        overlay(map,contour)
        overlay(map,vector)

        draw(map)        
        frame(wks)
end

出来的图效果如下：
对于出来的图但是总感觉有问题，我想有可能是对vibeta这个积分函数不够了解，但具体的原因我也无从得之，苦思冥想些时日未果之后，希望大神来帮忙解答。

BDRUSH

搜索了论坛里面，我发现没有NCL整层积分方面的贴子，很是苦恼，希望有志之士一起来研究到底用NCL该如何画出整层水汽通量。

xiahou90

mark 我也正在做

freekiller

好长啊，有没有和grads画的比一下？

BDRUSH

freekiller 发表于 2014-9-25 20:51
好长啊，有没有和grads画的比一下？

grads不太会啊，不然也不用这么纠结了，关于NCL积分的vibeta函数与grads的vint不知道有何区别，造成结果的误差

sun_shine_Xia

没在NCL里面算过~都是用grads算出来，在NCL里面画图的~

BDRUSH

sun_shine_Xia 发表于 2014-9-25 22:00
没在NCL里面算过~都是用grads算出来，在NCL里面画图的~

所以我想试试在NCL里面能算出来不~开了这个帖子

longlivehj

水汽通量、散度啥的，如果觉得不对，可以考虑检查下输入数据单位。
另外，你这程序运行时间应该不短吧！程序里的循环都是没有必要的，尽量用数组直接运算。

BDRUSH

longlivehj 发表于 2014-9-25 23:13
水汽通量、散度啥的，如果觉得不对，可以考虑检查下输入数据单位。
另外，你这程序运行时间应该不短吧！程 ...

数据应该没有问题，都进行过了验证，另外程序也运行得挺快的，就纳闷是不是vibeta函数理解的问题

阿吉木柯

不知道necp里面有没有 total water vapor..你可以用vebeta函数算整层水汽 然后和twv对比一下