- 积分
- 12311
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-11-23
- 最后登录
- 1970-1-1
|
发表于 2020-5-10 11:23:22
|
显示全部楼层
码一下我改的楼主的脚本，0：7层不是我们需要的1000：300，要单独从里面挑出来。
; http://bbs.06climate.com/forum.php?mod=viewthread&tid=28890
; http://bbs.06climate.com/forum.php?mod=viewthread&tid=35437
; 资料是FNl的grib2，调了一晚上都没成功，跪求大神帮忙，红色是报错
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
;************************************************
; open file and read in variable
;***********************************************
;levels 1000 925 850 700 600 500 400 300
; 层数 0 1 2 3 4 5 6 7
;***********************************************
nlat = 181
nlon = 360
nyear = 1
nlev = 8
grb_file = addfile("/home/sma/data/work/ping/fnl_20190620_00_00.grib2","r")
q = grb_file->SPFH_P0_2L108_GLL0(::-1,:)
u = grb_file->UGRD_P0_L100_GLL0(:,::-1,:)
v = grb_file->VGRD_P0_L100_GLL0(:,::-1,:)
slp = grb_file->PRES_P0_L1_GLL0(::-1,:)
;printVarSummary(q)
u8 = u({30000},:,:)
u7 = u({40000},:,:)
u6 = u({50000},:,:)
u5 = u({60000},:,:)
u4 = u({70000},:,:)
u3 = u({85000},:,:)
u2 = u({92500},:,:)
u1 = u({100000},:,:)
uu = new((/nlev,nlat,nlon/),float)
uu = (/u1,u2,u3,u4,u5,u6,u7,u8/)
v8 = v({30000},:,:)
v7 = v({40000},:,:)
v6 = v({50000},:,:)
v5 = v({60000},:,:)
v4 = v({70000},:,:)
v3 = v({85000},:,:)
v2 = v({92500},:,:)
v1 = v({100000},:,:)
vv = new((/nlev,nlat,nlon/),float)
vv = (/v1,v2,v3,v4,v5,v6,v7,v8/)
lev = (/1000,925,850,700,600,500,400,300/)
lat = fspan(90,-90,181)
lon = fspan(0,360,360)
lon!0 = "lon"
lon@units = "degrees_east"
lat!0 = "lat"
lat@units = "degrees_north"
lev!0 = "lev"
lev@units = "hPa"
uu!0 = "lev"
uu!1 = "lat"
uu!2 = "lon"
uu&lon = lon
uu&lat = lat
uu&lev = lev
vv!0 = "lev"
vv!1 = "lat"
vv!2 = "lon"
vv&lon = lon
vv&lat = lat
vv&lev = lev
;printVarSummary(uu) ;(uu(:,{38:43},{92:100}))
qu1 = new((/nlev,nlat,nlon/),float)
qv1 = new((/nlev,nlat,nlon/),float)
do iz=0,7
qu1(iz,:,:) = uu(iz,:,:)*q(:,:)
qv1(iz,:,:) = vv(iz,:,:)*q(:,:)
end do
copy_VarMeta(uu, qu1)
copy_VarMeta(vv, qv1)
qu1!0 = "lev"
qu1!1 = "lat"
qu1!2 = "lon"
qu1&lon = lon
qu1&lat = lat
qu1&lev = lev
qv1!0 = "lev"
qv1!1 = "lat"
qv1!2 = "lon"
qv1&lon = lon
qv1&lat = lat
qv1&lev = lev
; print(qu1(:,{38:43},{92:100}))
; print(qv1(:,{38:43},{92:100}))
; exit
;qu1 = latFlip
;asciiwrite("qu1.txt",qu1(0:7,{38:43},{92:100}))
;print(qu1(0:7,{38:43},{92:100}))
;exit
;**************************************************************************************
qu11 = qu1(lat|:, lon|:, lev|:)
delete(qu1)
qv11 = qv1(lat|:, lon|:, lev|:)
delete(qv1)
p = (/ 1000.,925.,850.,700.,600.,500.,400.,300. /)
linlog = 2
pbot = 1000.
ptop = 300.
; 1为线性积分 2为对数积分
data11 = (vibeta(p,qu11,linlog,slp,1000,300))/9.8 ; 返回的数值是2维的(lat,lon)
data12 = (vibeta(p,qv11,linlog,slp,1000,300))/9.8
qu = new((/nlat,nlon/),float)
qv = new((/nlat,nlon/),float)
mm = new((/nlat,nlon/),float)
mm1 = sqrt(data11^2+data12^2)
mm = mm1
qu = data11
qv = data12
; copy_VarMeta(q, mm)
; copy_VarMeta(uu(0,:,:), qu)
; copy_VarMeta(vv(0,:,:), qv)
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"
; printVarSummary(qu)
; printVarSummary(qv)
; printVarSummary(mm)
print(qu({38:43},{92}))
print(qv({38},{92:100}))
print(qu({38:43},{100}))
print(qv({43},{92:100}))
; exit
wks = gsn_open_wks("pdf","vapour_flux")
gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
;***************************************************************************
vcres = True ; plot mods desired
vcres@gsnAddCyclic = False
vcres@gsnDraw = False ; don't draw yet
vcres@gsnFrame = False ; don't advance frame yet
vcres@vcMinFracLengthF = 0.33 ;设置了矢量对于参考矢量的长度的最小数量级(default=0)0.33
vcres@vcRefAnnoOrthogonalPosF = -1.0
vcres@vcRefAnnoArrowLineColor = "black"
vcres@vcRefMagnitudeF = 20 ;20 风速矢量长度参考
vcres@vcRefLengthF = 0.045 ;0.045 0.02
vcres@vcGlyphStyle = "CurlyVector" ; turn on curly vectors
vcres@vcLineArrowThicknessF = 1.5
;vcres@vcRefAnnoOn = True
vcres@vcRefAnnoArrowUseVecColor = False ; don't use vec color for ref
vcres@vcLineArrowColor = "black" ; change vector color
vcres@vcVectorDrawOrder = "PostDraw"
vcres@vcRefAnnoString1 = "20"
vcres@vcRefAnnoSide = "Top"
vcres@vcRefAnnoString2On = False
vcres@vcRefAnnoPerimOn = False
vcres@vcRefAnnoOrthogonalPosF = -0.105
vcres@vcRefAnnoParallelPosF = 0.999
vector = gsn_csm_vector(wks,qu,qv,vcres) ; create plot
;***************************************************************************
;
;***************************************************************************
res = True ; plot mods desired
res@gsnMaximize = False
res@gsnAddCyclic = False
res@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
res@mpDataSetName = "Earth..4" ; This new database contains
res@mpDataBaseVersion = "MediumRes" ; Medium resolution database
res@mpAreaMaskingOn = True
res@mpMaskAreaSpecifiers = (/"China"/)
res@mpOutlineSpecifiers = (/"China","China:Provinces"/)
res@mpLandFillColor = "white"
res@mpOceanFillColor = "white"
res@mpFillBoundarySets = "NoBoundaries"
res@mpOutlineBoundarySets = "NoBoundaries"
res@mpNationalLineColor = "black"
res@mpProvincialLineColor = "black"
res@mpGeophysicalLineColor = "black"
res@mpGeophysicalLineThicknessF = 2.
res@mpNationalLineThicknessF = 2.
res@mpOutlineOn = True ; Turn on map outlines
res@mpLimitMode = "LatLon"
res@mpMinLonF = 75 ;92
res@mpMaxLonF = 135 ;100
res@mpMinLatF = 25 ;38
res@mpMaxLatF = 50 ;43
res@mpShapeMode = "FreeAspect"
res@vpHeightF = 0.55
res@vpWidthF = 0.8
res@pmTickMarkDisplayMode = "Always"
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 = (/20,24,28,32,36,40,44/) ;(/14,18,22,26,30,34,38,42,46,50,52/)
res@gsnLeftString =""
res@gsnRightString =""
res@tiMainString =""
res@lbLabelBarOn = True
res@lbLabelOffsetF = 0.06
res@pmLabelBarWidthF = 0.1
res@lbOrientation = "Vertical" ; vertical label bars
res@lbAutoManage = True
res@lbLabelFontHeightF = 0.009 ; make labels smaller
mm = smth9(mm,0.5,0.25,False)
;print(mm(128:133,92:100))
;mm = mm*100
;print(mm(128:133,92:100))
contour = gsn_csm_contour_map(wks,mm,res)
overlay(contour,vector)
x=(/92,100,100,92,92/)
y=(/43,43,38,38,43/)
lnres = True
lnres@gsLineColor="red"
lnres@gsLineThicknessF= 4
lnres@gsLineDashPattern= 0
dum=gsn_add_polyline(wks,contour,x,y,lnres)
draw(contour)
frame(wks)
end |
|
窥视卡
雷达卡
显身卡
楼主