NCL画水汽通量图出来了，但是有问题，不知道哪里出错了，求指点！！！！

xxc

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"

begin
;读取文件
f1=addfile("/mnt/f/shum.mon.mean.nc","r") ;比湿数据
;print(f1)
f2=addfile("/mnt/f/uwnd.mon.mean.nc","r")
f3=addfile("/mnt/f/vwnd.mon.mean.nc","r")
f4=addfile("/mnt/f/pres.sfc.mon.mean .nc","r")
;print(f4)
;print(f1)
shum  = f1->shum(1664:1666,0:14,:,:)
;printVarSummary(shum)
uwnd  =f2->uwnd(1664:1666,0:14,:,:)
;printVarSummary(uwnd)
vwnd  = f3->vwnd(1664:1666,0:14,:,:)
pres  = f4->pres(1664:1666,:,:)
;print(pres)
lev=f1->level(0:14)
;print(lev)
lat=f1->lat(:)
lon=f1->lon(:)
;print(lat)
;print(lon)

qu=uwnd*shum
  qv=vwnd*shum
  copy_VarCoords(vwnd,qv)
  copy_VarCoords(uwnd,qu)

  v1 =  new((/15,91,180/), "float")
  u1 =  new((/15,91,180/), "float")

u1=dim_avg_n_Wrap(uwnd,0)
v1=dim_avg_n_Wrap(vwnd,0)
u2=dim_avg_n_Wrap(u1,0)
v2=dim_avg_n_Wrap(v1,0)

  ;print(qu2)
;d_qu= dimsizes(qu) ; 67 8 73 144

Qu1=qu;
Qv1=qv;
lev1=lev

;print(dimsizes(Qu1))
Qu1!0="time"
Qu1!1="level"
Qu1!2="lat"
Qu1!3="lon"

Qv1!0="time"
Qv1!1="level"
Qv1!2="lat"
Qv1!3="lon"  

Qv1&lon=(/lon/)
Qv1&lat=(/lat/)
Qv1&level=(/lev/)

Qu1&lon=(/lon/)
Qu1&lat=(/lat/)
Qu1&level=(/lev1/)

;print(mean_anoQu1(:,1,1))
;Qu=new((/d_qu(1),d_qu(2),d_qu(0)/),float)
;Qv=new((/d_qu(1),d_qu(2),d_qu(0)/),float)

Qu= Qu1(time|:,lat|:,lon|:,level|:)  
Qv= Qv1(time|:,lat|:,lon|:,level|:)  
psfc=pres(time|:,lat|:,lon|:)
qu_new=new((/3,91,180,15+4/),double)
qv_new=new((/3,91,180,15+4/),double)
do i=0,14
qu_new(:,:,:,i)= Qu(:,:,:,i)
qv_new(:,:,:,i)= Qv(:,:,:,i)
end do

do j=15,18
qu_new(:,:,:,j)= 0.0
  qv_new(:,:,:,j)= 0.0
end do
;print(Qu(1,1,:))
pa=(/1000,950,900,850,800,750,700,650,600,550,500,450,400,350,300,250,200,150,100/)
;qqu=new((/73,144,67/),double)
;qqv=new((/73,144,67/),double)

;计算整层水汽通量
ptop=300
pbot=1000
qqu=(vibeta(pa,qu_new,1,psfc,pbot,ptop))/9.8*100
qqv=(vibeta(pa,qv_new,1,psfc,pbot,ptop))/9.8*100
print(dimsizes(qqu))  
qu2=dim_avg_n_Wrap(qqu, 0)
qv2=dim_avg_n_Wrap(qqv, 0)
copy_VarCoords(uwnd(0,:,:,:),qu2)
copy_VarCoords(vwnd(0,:,:,:),qv2)
print(qu2)



  lat1              =  fspan(-10,40,91)    ;设定全球图的网格点经纬度信息
  lat1!0            =  "lat"
  lat1@units        =  "degrees_north"
  lat1@long_name    =  "longitude"
  lon1              =  fspan(40,130,180)
  lon1!0            =  "lon"
  lon1@units        =  "degrees_east"
  lon1@long_name    =  "longitude"

  qu2!0   =  "lat"
  qu2!1   =  "lon"
  qu2&lat =  (/lat1/)
  qu2&lon =  (/lon1/)
  qv2!0   =  "lat"
  qv2!1   =  "lon"
  qv2&lat =  (/lat1/)
  qv2&lon =  (/lon1/)



wks=gsn_open_wks("x11","uv200")
res1 = True
res1@gsnAddCyclic = False
res1@vcGlyphStyle = "CurlyVector" ;设置箭头参考样式
res1@vcRefMagnitudeF =200.0 ;矢量箭头参考量
res1@vcRefLengthF = 0.045 ;矢量箭头参考长度
res1@vcMinDistanceF = 0.03 ;间隔距离
res1@gsnMajorLatSpacing = 5;设置坐标纬度间隔
res1@gsnMajorLonSpacing = 10;设置坐标经度间隔
res1@mpOutlineOn             =  False
res1@mpMaxLatF = 40
res1@mpMinLatF = -10
res1@mpMaxLonF = 130
res1@mpMinLonF = 40
res1@mpCenterLonF = 180
vector = gsn_csm_vector_map(wks,qu2,qv2,res1)
end

凉快了吧

请问您图出了吗，图是啥样啊

眼里有星辰

总得把图放出来看看吧

guoguohh

想问一下 你遇到过报错warning:vibeta: there must be at least three levels with data above the surface的情况吗，我用你的脚本一直报错

李乐乐

guoguohh 发表于 2021-11-9 21:22
想问一下 你遇到过报错warning:vibeta: there must be at least three levels with data above the surface ...

我也一直这样报错，请问你解决了吗？

李乐乐

guoguohh 发表于 2021-11-9 21:22
想问一下 你遇到过报错warning:vibeta: there must be at least three levels with data above the surface ...

感觉是地形的原因

guoguohh

李乐乐 发表于 2021-11-15 16:03
感觉是地形的原因

找到原因了，我用era5下载的比湿数据存在add_offset和scale_factor，直接用的话会报错不足三层，用short2flt 转成实际数据就不报错了