- 积分
- 4112
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2021-1-14
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 guoguohh 于 2021-11-20 19:49 编辑
看到家园里几乎都是用vibeta求垂直积分的
而era5用vibeta函数求垂直积分函数一直报错，师姐在混合坐标系下使用vibeta就没有问题，
我按照一样的步骤用就一直报错“p没有从bottom to top”，但是实际上是符合要求的，
改了很多次也并没有解决，所以放弃用vibeta。
Dennis在ncl-talk里面回复说“"vibeta" was one of NCL's original functions. I would like
to see its use discontinued.”https://www.ncl.ucar.edu/Support/talk_archives/2010/0787.html不过也没说清楚原因，可能大家以后尽量就不要用vibeta了。
最后修改了官网给的画水汽通量里的方法画出了图，官网脚本链接https://www.ncl.ucar.edu/Applications/Scripts/mfc_div_1.ncl
附我的脚本和图：
begin
fu = addfile("ujja1979-2018.nc","r");m/s
fv = addfile("vjja1979-2018.nc","r");m/s
fq = addfile("shum.nc","r");kg/kg
fps = addfile("pres.nc", "r");Pa=(kg/(m*s2))
pp=fv->level({1000:300});hPa
p=pp*100;Pa
u =fu->uwnd(:,{1000:300},:,:);m/s eofts_high
v =fv->vwnd(:,{1000:300},:,:);m/s
q =short2flt(fq->q(:,{1000:300},:,:));kg/kg
ps =short2flt(fps->sp(:,:,:));Pa=(kg/(m*s2))
;--------------计算水汽通量 moisture flux=(dp(Pa)*u(m/s)*q(kg/kg))/g(m/s^2) ----------------------
ptop=30000
g= 9.80665 ; m/s2
dp = dpres_plevel_Wrap(p, ps, ptop, 0);这里变量单位必须全部统一成pa或者hpa
dpg = dp/g
dpg@long_name = "Layer Mass Weighting"
dpg@units = "kg/m2"
copy_VarCoords(u,dpg)
uq = u*q ; (:,:,:,:)
uq@long_name = "Zonal Moisture Flux [uq]"
uq@units = "["+u@units+"]["+q@units+"]"
copy_VarCoords(u,uq)
vq = v*q ; (:,:,:,:)
vq@long_name = "Meridional Moisture Flux [vq]"
vq@units = "["+v@units+"]["+q@units+"]"
copy_VarCoords(v,vq)
uq_dpg = uq*dpg
iuq = dim_avg_n_Wrap(dim_sum_n(uq_dpg, 1),0);kg/m*s
vq_dpg = vq*dpg
ivq = dim_avg_n_Wrap(dim_sum_n(vq_dpg, 1),0);kg/m*s
copy_VarCoords(u(0,0,:,:), iuq(:,:))
copy_VarCoords(u(0,0,:,:), ivq(:,:))
printMinMax(iuq,0)
printMinMax(ivq,0)
wks = gsn_open_wks("x11","flux") ; send graphics to PNG file
res = True ; plot mods desired
res@gsnLeftString = "Wind"
res@gsnDraw = False ; don't draw yet
res@gsnFrame =False ; don't advance frame yet
res@gsnAddCyclic = False ; plotted dataa are not cyclic
res@gsnLeftString=" "
res@gsnRightString=" "
res@pmTickMarkDisplayMode="Always";在坐标标签上添加°符号
res@mpFillOn = True ; turn off map fill
res@mpMinLatF = 0. ; zoom in on map
res@mpMaxLatF =60.
res@mpMinLonF =30.
res@mpMaxLonF = 145.
res@mpDataSetName = "Earth..4" ; This new database contains; divisions for other countries.
res@mpDataBaseVersion = "MediumRes" ; Medium resolution database
res@mpOutlineOn = True ; Turn on map outlines
res@mpAreaMaskingOn = True
res@mpMaskAreaSpecifiers = (/"land"/)
res@mpOutlineSpecifiers = (/"China"/)
res@mpLandFillColor = "white"
res@mpInlandWaterFillColor = "white"
res@mpOceanFillColor = "white"
res@vcRefMagnitudeF =100. ; define vector ref mag
res@vcRefLengthF = 0.04 ; define length of vec ref
res@vcGlyphStyle = "CurlyVector" ; turn on curly vectors
res@vcMinDistanceF = 0.01
res@vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector up
base = gsn_csm_vector_map(wks,iuq, ivq, res)
draw(base)
frame(wks)
end
================分割线===========================================
老师们好
最近在画整层水汽通量图，出现了一些问题想请大家帮忙看看
数据来自era5的月均数据，引入变量的时候提前计算了夏季的年平均，
我用era5下载的比湿数据存在add_offset和scale_factor，直接用的话会报错不足三层，使用short2flt 转成实际数据（单位是kg/kg），
问题：绘图问题，选择用自己计算出的比湿计算整层水汽通量，出的图是这样的图1，而别人的图是图2，需要的趋势全都看不出来，
别人同样是气候态的图，差别太大了，贴一下我的脚本，麻烦大家帮我看看是哪里出错了
begin
fu = addfile("ujja1979-2018.nc","r");m/s
fv = addfile("vjja1979-2018.nc","r");m/s
fq = addfile("shum.nc","r");kg/kg
fps = addfile("pres.nc", "r");Pa=(kg/(m*s2))
ft =addfile("airjja.nc", "r");K
frh =addfile("rhumjja.nc", "r");%
ys=1979
ye=2010
Time = fu->time
year = cd_calendar(Time,-1)/100 ;获得年变量
iyear = ind(year.ge.ys .and. year .le.ye) ;找到年代为1979到2010的数据位置
u =fu->uwnd(iyear,:,:,:);m/s
v =fv->vwnd(iyear,:,:,:);m/s
ps =short2flt(fps->sp(iyear,:,:));Pa=(kg/(m*s2))
t=ft->air(iyear,:,:,:)
rh=frh->rhum(iyear,:,:,:)
p=fv->level
p@units="hPa"
q = mixhum_ptrh(conform(t,p,1), t ,rh, 2 );kg/kg 2 = specific humidity (kg/kg).
copy_VarMeta(v,q)
q@units="kg/kg"
printMinMax(q,False)
pbot=1000
ptop=300
g= 9.80665 ; m/s2
uq = u*q ; (:,:,:,:)
copy_VarCoords(u,uq) ; (time,level,lat,lon)
;printMinMax(uq,False)
vq = v*q ; (:,:,:,:)
copy_VarCoords(v,vq) ; (time,level,lat,lon)
vint_qu=100*vibeta(p,uq(time|:,lat|:,lon|:,level|:),1,ps,pbot,ptop)/g
copy_VarMeta(ps,vint_qu)
vint_qv=100*vibeta(p,vq(time|:,lat|:,lon|:,level|:),1,ps,pbot,ptop)/g
copy_VarMeta(ps,vint_qv)
mvq=dim_avg_n_Wrap(vint_qv,0)
muq=dim_avg_n_Wrap(vint_qu,0)
printMinMax(mvq, False)
printMinMax(muq, False)
wks = gsn_open_wks("X11","X1") ; send graphics to PNG file
res = True ; plot mods desired
res@gsnLeftString = "Wind"
res@gsnDraw = False ; don't draw yet
res@gsnFrame =False ; don't advance frame yet
res@gsnAddCyclic = False ; plotted dataa are not cyclic
res@gsnLeftString=" "
res@gsnRightString=" "
res@pmTickMarkDisplayMode="Always";在坐标标签上添加°符号
res@mpFillOn = True ; turn off map fill
res@mpMinLatF = 0. ; zoom in on map
res@mpMaxLatF =60.
res@mpMinLonF =50.
res@mpMaxLonF = 145.
res@mpDataSetName = "Earth..4" ; This new database contains; divisions for other countries.
res@mpDataBaseVersion = "MediumRes" ; Medium resolution database
res@mpOutlineOn = True ; Turn on map outlines
res@mpAreaMaskingOn = True
res@mpMaskAreaSpecifiers = (/"land"/)
res@mpOutlineSpecifiers = (/"China"/)
res@mpLandFillColor = "white"
res@mpInlandWaterFillColor = "white"
res@mpOceanFillColor = "white"
res@vcRefMagnitudeF =100. ; define vector ref mag
res@vcRefLengthF = 0.04 ; define length of vec ref
res@vcGlyphStyle = "CurlyVector" ; turn on curly vectors
res@vcMinDistanceF = 0.02
res@vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector up
base = gsn_csm_vector_map(wks,muq,mvq,res)
draw(base)
frame(wks)
end
|
-
图1
-
图2
|
窥视卡
雷达卡
发表于 2021-11-12 12:30:10
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
{:eb511:}{:eb511:}{:eb511:}{:eb511:}{:eb511:}
楼主