分享：水汽通量散度

1649518749

;图形和grads、MI对比了一下，量级都是1e-5，grads大小要比ncl、MI偏大4倍左右，估计是里面函数风速单位的问题。
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/cnmap/cnmap.ncl"

begin

   a=addfile("/mnt/d/wrfmodel/2014/air.2014.nc","r")
   b=addfile("/mnt/d/wrfmodel/2014/uwnd.2014.nc","r")
   c=addfile("/mnt/d/wrfmodel/2014/vwnd.2014.nc","r")
   d=addfile("/mnt/d/wrfmodel/2014/rhum.2014.nc","r")
   air= short2flt(a->air(0,{700},:,:))
   u=short2flt(b->uwnd(0,{700},:,:))
   v=short2flt(c->vwnd(0,{700},:,:))
   rhum=short2flt(d->rhum(0,{700},:,:))
   lat=a->lat
   lon=a->lon
   es = 6.112*exp(17.67*(air-273.16)/(air-29.65))
   qs = 0.62197*es*1000/(700-0.378*es)
   q = qs*rhum/100
q = smth9(q, 0.50, -0.25, False)
copy_VarCoords(air,es)
  copy_VarCoords(air,qs)
   copy_VarCoords(air,q)
   qu = new((/73,144/),float)
   qv = new((/73,144/),float)
                do i=0,72
                        do j=0,143
                                qu(i,j) = q(i,j)*u(i,j)/9.8
                                qv(i,j) = q(i,j)*v(i,j)/9.8
                        end do
                end do
   copy_VarCoords(air,qu)
   copy_VarCoords(air,qv)
    ;    qu = smth9(qu, 0.50, -0.25, False)
    ;    qv = smth9(qv, 0.50, -0.25, False)
        vapord = new((/73,144/),float)
   vapord=uv2dv_cfd(qu,qv,lat,lon,2)               ; u,v ==> divergence
   vapord = vapord*100000
copy_VarCoords(air,vapord)
printVarSummary(vapord)
;>=========================================================<
wks       = gsn_open_wks("pdf","vaporsandu")
gsn_define_colormap( wks ,"BlWhRe")  


res                         = True            
res@tiMainString            =""
res@gsnMaximize             = False
res@gsnDraw                 = False
res@gsnFrame                = False

;>--------------------------------------------<
;            set for the map
;>--------------------------------------------<
res@mpMinLatF               = 15.                        
res@mpMaxLatF               = 55.
res@mpMinLonF               = 72.
res@mpMaxLonF               = 136.
res@pmTickMarkDisplayMode = "Always"
res@mpFillOn                = True
res@mpOutlineOn             = False  ; Use outlines from shapefile

res@mpDataBaseVersion       = "MediumRes"
res@mpDataSetName           = "Earth..4"
res@mpAreaMaskingOn         = True
res@mpMaskAreaSpecifiers    = (/"China","Taiwan","Disputed area between India and China","India:Arunachal Pradesh"/)
res@mpLandFillColor         = "white"
res@mpInlandWaterFillColor  = "white"
res@mpOceanFillColor        = "white"
res@mpOutlineBoundarySets   = "NoBoundaries"
;>--------------------------------------------<
; set for the plot
res@cnInfoLabelOn = False
res@gsnLeftString  = ""
res@gsnRightString =""
res@cnFillOn                = True               
res@cnLinesOn               = False
res@cnLineLabelsOn         = False   
res@cnLineColor = "blue"  
res@cnLevelSelectionMode = "ManualLevels"
res@cnMinLevelValF   = -2.5 ; set min contour level
res@cnMaxLevelValF=  3. ; set max contour level            
res@gsnSpreadColors         = True         
res@lbLabelAutoStride       = True      
;res@lbBoxEndCapStyle = "TriangleBothEnds"   
;res@cnFillDrawOrder         = "PreDraw"     
res@cnLevelSpacingF= 0.25
;res@gsnContourPosLineDashPattern=0  ;>0
;res@gsnContourNegLineDashPattern=1 ;<0线的种类
;res@gsnContourZeroLineThicknessF=2.
map=gsn_csm_contour_map(wks,vapord,res)

;>============================================================<
;                      add China map
;>------------------------------------------------------------<;
cnres           = True
cnres@china     = True       ;draw china map or not
cnres@river     = False       ;draw changjiang&huanghe or not
cnres@province  = True       ;draw province boundary or not
cnres@nanhai    = False       ;draw nanhai or not
cnres@diqu      = False       ; draw diqujie or not

chinamap = add_china_map(wks,map,cnres)
;>============================================================<

draw(map)
frame(wks)
end  

1649518749

大南瓜 发表于 2017-10-15 12:34
亲，想问一下你这个水汽通量散度的单位是什么呀？

g*cm-2*hpa-1*s-1。

鱼鱼鱼鱼

请问大神，我使用wrfout结果进行画图直接读取（时间，高度层，纬度，经度）数据，但是wrfout输出的是格点数，不是按照经纬度进行区分的，在 vapord=uv2dv_cfd(qu,qv,lat,lon,2) 这一步应该如何设置呢？

1649518749

鱼鱼鱼鱼 发表于 2017-6-8 11:02
请问大神，我使用wrfout结果进行画图直接读取（时间，高度层，纬度，经度）数据，但是wrfout输出的是格点数 ...

用wrf_user_getvar 读取就行了。

鱼鱼鱼鱼

1649518749 发表于 2017-6-8 12:42
用wrf_user_getvar 读取就行了。

wrf_user_getVar好像是在高度上不分层的，没办法求某一层的水汽通量散度

1649518749

鱼鱼鱼鱼 发表于 2017-6-8 15:43
wrf_user_getVar好像是在高度上不分层的，没办法求某一层的水汽通量散度

wrf_user_intrp3d  再配合这个函数，就可以读出某个气压层的数据

Soaring

奉献精神很强嘛

1649518749

Soaring 发表于 2017-6-9 10:26
奉献精神很强嘛

赚积分{:lol:}

鱼鱼鱼鱼

1649518749 发表于 2017-6-8 17:14
wrf_user_intrp3d  再配合这个函数，就可以读出某个气压层的数据

谢谢！最后我自己定义的经度纬度，也是可以出图的。

初翠翠

感谢分享，受益匪浅

tcbyle

{:eb502:}