整层水汽通量的讨论

日不不不落 · 发表于 2015-12-28 09:26:07

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

关于水汽通量的帖子查了好多，根据大家的帮助也使用NCL画出来了，但是感觉出来的图有点怪怪的，希望大神们能帮忙看看，可能是哪里出错了，关于积分上限和下限值的设置，是要遵循什么规律吗

日不不不落 · 发表于 2015-12-28 10:07:57

andrewsoong 发表于 2015-12-28 10:02
的确是有点分辨率不高哇，你可以贴上脚本让大家给你看看吧

;2014年台风季——整层水汽输送
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"
begin
latS = 0
latN = 60
lonL = 60
lonR = 180
r1 = "/gpfsES/geo/zywang/get_flux_myan/vertical/"
r2 = "/gpfsES/geo/zywang/data_CESM/"
n1 = "U_on_p.nc"
n2 = "V_on_p.nc"
n3 = "SHUM_on_p.nc"
n4 = "MASK_SV/T_SV.nc"
f1 = addfile(r1+n1,"r")
f2 = addfile(r1+n2,"r")
f3 = addfile(r1+n3,"r")
f4 = addfile(r2+n4,"r")
u = f1->U
lev_p = f3->lev_p
lat = f4->lat
lon = f4->lon
v = f2->V
shum = f3->SHUM
lev_p = f3->lev_p
pres = f4->PS ;suface pressure(Pa)
;***********************
summer_u = new((/2000,10,48,96/),float)
summer_v = new((/2000,10,48,96/),float)
summer_shum = new((/2000,10,48,96/),float)
summer_pres = new((/2000,48,96/),float)
jj = 0;
do ii = 0,23988,12
summer_u(jj,:,:,:) = dim_avg_n_Wrap(u(ii+4:ii+7,:,:,:),0)
summer_v(jj,:,:,:) = dim_avg_n_Wrap(v(ii+4:ii+7,:,:,:),0)
summer_shum(jj,:,:,:) = dim_avg_n_Wrap(shum(ii+4:ii+7,:,:,:),0)
summer_pres(jj,:,:) = dim_avg_n_Wrap(pres(ii+4:ii+7,:,:),0)
jj = jj+1;
end do
summer_u11 = summer_u(:,:,23:40,16:48)
summer_v11 = summer_v(:,:,23:40,16:48)
summer_shum11 = summer_shum(:,:,23:40,16:48)
summer_pres11 = summer_pres(:,23:40,16:48)
; add attributes
year = ispan(1,2000,1)
summer_u11!0 = "time"
summer_u11&time = year
summer_u11!1 = "lev_p"
summer_u11&lev_p = lev_p
summer_u11!2 = "lat"
summer_u11&lat = lat(23:40)
summer_u11!3 = "lon"
summer_u11&lon = lon(16:48)
summer_v11!0 = "time"
summer_v11&time = year
summer_v11!1 = "lev_p"
summer_v11&lev_p = lev_p
summer_v11!2 = "lat"
summer_v11&lat = lat(23:40)
summer_v11!3 = "lon"
summer_v11&lon = lon(16:48)
summer_shum11!0 = "time"
summer_shum11&time = year
summer_shum11!1 = "lev_p"
summer_shum11&lev_p = lev_p
summer_shum11!2 = "lat"
summer_shum11&lat = lat(23:40)
summer_shum11!3 = "lon"
summer_shum11&lon = lon(16:48)
summer_pres11!0 = "time"
summer_pres11&time = year
summer_pres11!1 = "lat"
summer_pres11&lat = lat(23:40)
summer_pres11!2 = "lon"
summer_pres11&lon = lon(16:48)
print(max(summer_pres11))
print(min(summer_pres11))
printVarSummary(summer_pres11)
printVarSummary(summer_pres11)
;calculate qu = shum(specific humidity)*u(wind)
qu = new(dimsizes(summer_shum11),"float")
qu = summer_shum11*summer_u11
copy_VarMeta(summer_shum11,qu)
; calculate qv = shum(specific humidity)*v(wind)
qv = new(dimsizes(summer_shum11),"float")
qv = summer_shum11*summer_v11
copy_VarMeta(summer_shum11,qv)
printVarSummary(qu)
printVarSummary(qv)
; vertically integrated water vapor transport(kg/(m.s))
linlog=1
pbot=1000 (right?)
ptop=100 (right?)
g=9.8
vint_qu=vibeta(lev_p,qu(time|:,lat|:,lon|:,lev_p|:),linlog,summer_pres11,pbot,ptop)/g
copy_VarMeta(qu(:,0,:,:),vint_qu)
vint_qv=vibeta(lev_p,qv(time|:,lat|:,lon|:,lev_p|:),linlog,summer_pres11,pbot,ptop)/g
copy_VarMeta(qv(:,0,:,:),vint_qv)
vint_qu_W1 = vint_qu(270:369,:,:)
vint_qu_W2 = vint_qu(780:879,:,:)
vint_qu_W3 = vint_qu(1069:1168,:,:)
vint_qu_W4 = vint_qu(1900:1999,:,:)
vint_qu_C1 = vint_qu(610:709,:,:)
vint_qu_C2 = vint_qu(1440:1539,:,:)
vint_qu_C3 = vint_qu(1639:1738,:,:)
vint_qu_C4 = vint_qu(1760:1859,:,:)
vint_qu_W = (vint_qu_W1+vint_qu_W2+vint_qu_W3+vint_qu_W4)/4.0
vint_qu_C = (vint_qu_C1+vint_qu_C2+vint_qu_C3+vint_qu_C4)/4.0
avg_vint_qu = dim_avg_n_Wrap(vint_qu,0)
vint_qv_W1 = vint_qv(270:369,:,:)
vint_qv_W2 = vint_qv(780:879,:,:)
vint_qv_W3 = vint_qv(1069:1168,:,:)
vint_qv_W4 = vint_qv(1900:1999,:,:)
vint_qv_C1 = vint_qv(610:709,:,:)
vint_qv_C2 = vint_qv(1440:1539,:,:)
vint_qv_C3 = vint_qv(1639:1738,:,:)
vint_qv_C4 = vint_qv(1760:1859,:,:)
vint_qv_W = (vint_qv_W1+vint_qv_W2+vint_qv_W3+vint_qv_W4)/4.0
vint_qv_C = (vint_qv_C1+vint_qv_C2+vint_qv_C3+vint_qv_C4)/4.0
avg_vint_qv = dim_avg_n_Wrap(vint_qv,0)
do pp = 0,99,1
vint_qu_W(pp,:,:) = vint_qu_W(pp,:,:)-avg_vint_qu(:,:);
vint_qu_C(pp,:,:) = vint_qu_C(pp,:,:)-avg_vint_qu(:,:);
vint_qv_W(pp,:,:) = vint_qv_W(pp,:,:)-avg_vint_qv(:,:);
vint_qv_C(pp,:,:) = vint_qv_C(pp,:,:)-avg_vint_qv(:,:);
end do
printVarSummary(vint_qu_W)
printVarSummary(vint_qu_C)
printVarSummary(vint_qv_W)
printVarSummary(vint_qv_C)
qu_W = dim_avg_n_Wrap(vint_qu_W,0)
qu_C = dim_avg_n_Wrap(vint_qu_C,0)
qv_W = dim_avg_n_Wrap(vint_qv_W,0)
qv_C = dim_avg_n_Wrap(vint_qv_C,0)
print(min(qu_W))
print(max(qv_W))
print(min(qu_C))
print(max(qv_C))
copy_VarMeta(vint_qu_W(0,:,:),qu_W)
copy_VarMeta(vint_qu_C(0,:,:),qu_C)
copy_VarMeta(vint_qv_W(0,:,:),qv_W)
copy_VarMeta(vint_qv_C(0,:,:),qv_C)
mm_W = sqrt(qu_W^2+qv_W^2)
mm_C = sqrt(qu_C^2+qv_C^2)
copy_VarMeta(qu_W,mm_W)
copy_VarMeta(qu_C,mm_C)
print(max(mm_W))
print(min(mm_W))
print(max(mm_C))
print(min(mm_C))
wks =gsn_open_wks("eps","vapour_flux_W")
gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
vcres = True
vcres@gsnAddCyclic = False
vcres@gsnDraw = False ; don't draw yet
vcres@gsnFrame = False ; don't advance frame yet
vcres@mpDataSetName = "Earth..3" ; This new database contains
vcres@mpDataBaseVersion = "MediumRes" ; Medium resolution database
vcres@mpOutlineOn = True ; Turn on map outlines
vcres@mpProjection = "CylindricalEquidistant"
vcres@mpGeophysicalLineThicknessF = 0.7 ; double the thickness of geophysical boundaries
vcres@mpNationalLineThicknessF = 1.0 ; double the thickness of national boundaries
vcres@pmTickMarkDisplayMode = "Always"
vcres@mpMinLatF = latS
vcres@mpMaxLatF = latN
vcres@mpMinLonF = lonL
vcres@mpMaxLonF = lonR
vcres@lbLabelBarOn = True
vcres@vcRefMagnitudeF = 0.2 ; make vectors larger
vcres@vcRefLengthF = 0.05 ; ref vec length
vcres@vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector
vcres@vcRefAnnoArrowLineColor = "black" ; change ref vector color
vcres@vcRefAnnoArrowUseVecColor = False ; don't use vec color for ref
vcres@vcGlyphStyle = "CurlyVector" ; turn on curly vectors
;vcres@vcMinDistanceF = 0.017 ; thin out vectors
;vcres@vcRefAnnoOn = True
;vcres@vcMonoFillArrowFillColor = True
vcres@vcLineArrowThicknessF = 2.0
vcres@vcVectorDrawOrder= "PostDraw"
vcres@tiMainString ="vaper_flux : kg/(m*s)"
vcres@gsnRightString ="Rammasun(1409)"
plot=gsn_csm_vector_map(wks,qu_W,qv_W,vcres) ; create plot
cnres = True
cnres@china = True ;draw china map or not
cnres@river = False ;draw changjiang&huanghe or not
cnres@province = False ;draw province boundary or not
cnres@nanhai = False ;draw nanhai or not
cnres@diqu = False ; draw diqujie or not
chinamap = add_china_map(wks,plot,cnres)
txres = True
txres@txFontHeightF = 0.020 ; Set the font height
txres@txPerimOn = True
txres@txBackgroundFillColor =0
txres@txPerimColor = 0
res = True ; plot mods desired
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@gsnLeftString=""
res@gsnRightString=""
res@cnLevelSelectionMode="ManualLevels"
res@cnMinLevelValF = 0
res@cnMaxLevelValF = 0.14
res@cnLevelSpacingF = 0.014
;res@cnLevels=(/0,0.014,0.028,0.042,0.056,0.07,0.084,0.098,0.112,0.126,0.14/)
res@lbLabelBarOn = True
res@lbLabelStride = 2 ; plot every other colar bar label
res@lbOrientation = "vertical" ; vertical label bars
;res@lbLeftMarginF =-0.3
res@lbAutoManage = True
res@lbLabelFontHeightF = 0.009
;mm_W = smth9(mm_W,0.5,0.25,False)
plot1 = gsn_csm_contour(wks,mm_W,res)
overlay(plot,plot1)
mstring = "p"
fontnum = 37
xoffset = 0.0
yoffset = 0.0
ratio = 1.0
size = 1.0
angle = 0.0
new_index = NhlNewMarker(wks, mstring, fontnum, xoffset, yoffset, ratio, size, angle)
mkres = True
mkres@gsMarkerSizeF = 0.02
mkres@gsMarkerColor = "red"
mkres@gsMarkerThicknessF=5.0
mkres@gsMarkerIndex = new_index
mkres@tfPolyDrawOrder = "PostDraw"
;dum=gsn_add_polymarker(wks, plot, x(i), y(i), mkres)
maximize_output(wks,False)
end

复制代码

日不不不落 · 发表于 2015-12-28 09:27:28

画出来的矢量稀稀疏疏的不知道可能是哪里出错了呢望大神解救谢谢~

andrewsoong · 发表于 2015-12-28 09:57:59

你资料分辨率多少？

日不不不落 · 发表于 2015-12-28 10:00:32

andrewsoong 发表于 2015-12-28 09:57
你资料分辨率多少？

3.75*3.75的用这个资料单画风场的时候不是这么稀疏的

andrewsoong · 发表于 2015-12-28 10:02:56

日不不不落发表于 2015-12-28 10:00
3.75*3.75的用这个资料单画风场的时候不是这么稀疏的

的确是有点分辨率不高哇，你可以贴上脚本让大家给你看看吧

日不不不落 · 发表于 2015-12-28 10:15:05

上面是代码有点长跪求大神们解答谢谢~

日不不不落 · 发表于 2015-12-28 15:41:24

本帖最后由日不不不落于 2015-12-28 15:55 编辑

andrewsoong 发表于 2015-12-28 10:02
的确是有点分辨率不高哇，你可以贴上脚本让大家给你看看吧

解决了谢谢你，谢谢大家的帮助首先，气压的单位最好换算成Pa，其次，画图之前，变量的经纬度属性要设置！

姨妈血红 · 发表于 2016-1-8 10:31:21

日不不不落发表于 2015-12-28 10:07

请问下你的脚本里面的各个数据的单位是多少，最后得出的结果是什么量级的啊？

姨妈血红 · 发表于 2016-1-8 10:39:44

还有使用vibeta函数的时候linlog应为1还是2 啊？

		自动登录	找回密码
密码			立即注册

[作图] 整层水汽通量的讨论

登录后查看更多精彩内容~

浏览过的版块