- 积分
- 3025
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2019-1-19
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2019-8-31 22:13:12
|
显示全部楼层
脚本如下：
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
a = addfile("d:/cygwin/home/0228_87.5/wrfout_d03_2018-02-27_12_00_00.nc","r")
a1 = addfile("d:/cygwin/home/0401_87.5/wrfout_d03_2018-04-01_00_00_00.nc","r")
a2 = addfile("d:/cygwin/home/0513_87.5/wrfout_d03_2018-05-13_00_00_00.nc","r")
type = "pdf"
wks = gsn_open_wks(type,"wrf_panel_dv_vector1")
gsn_define_colormap(wks, "ncl_default")
res = True
res@MainTitle = ""
res@InitTime = False
res@ValidTime = False
res@Footer = False
pltres = True
pltres@NoTitles = True
pltres@PanelPlot = True
;set lat/lon
FirstTime = True
FirstTimeMap = True
times = wrf_user_getvar(a,"times",-1) ; 读取所有时次
ntimes = dimsizes(times) ; number of times in the file
plots = new (3, graphic)
mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
dNames = getfilevardims(a,"P") ; NCL 6.4.0 and earlier
nd = dimsizes(mdims) ; 查看mdims为几围数组
print(nd)
print(mdims)
mdims1 = getfilevardimsizes(a1,"P") ; get some dimension sizes for the file
dNames1 = getfilevardims(a1,"P") ; NCL 6.4.0 and earlier
nd1 = dimsizes(mdims1) ; 查看mdims为几围数组
print(nd1)
print(mdims1)
mdims2 = getfilevardimsizes(a2,"P") ; get some dimension sizes for the file
dNames2 = getfilevardims(a2,"P") ; NCL 6.4.0 and earlier
nd2 = dimsizes(mdims2) ; 查看mdims为几围数组
print(nd2)
print(mdims2)
;20180228
xlat = wrf_user_getvar(a, "XLAT",0)
xlon = wrf_user_getvar(a, "XLONG",0)
ter = wrf_user_getvar(a, "HGT",0)
u10 = wrf_user_getvar(a,"U10",0) ; u at 10 m, mass point
v10 = wrf_user_getvar(a,"V10",0) ; v at 10 m, mass point
u10@units = "m/s"
v10@units = "m/s"
lat = xlat(:,0)
lon = xlon(0,:)
;20180401
xlat1 = wrf_user_getvar(a1, "XLAT",0)
xlon1 = wrf_user_getvar(a1, "XLONG",0)
ter1 = wrf_user_getvar(a1, "HGT",0)
u101 = wrf_user_getvar(a1,"U10",0) ; u at 10 m, mass point
v101 = wrf_user_getvar(a1,"V10",0) ; v at 10 m, mass point
u101@units = "m/s"
v101@units = "m/s"
lat1 = xlat1(:,0)
lon1 = xlon1(0,:)
;20180513
xlat2 = wrf_user_getvar(a2, "XLAT",0)
xlon2 = wrf_user_getvar(a2, "XLONG",0)
ter2 = wrf_user_getvar(a2, "HGT",0)
u102 = wrf_user_getvar(a2,"U10",0) ; u at 10 m, mass point
v102 = wrf_user_getvar(a2,"V10",0) ; v at 10 m, mass point
u102@units = "m/s"
v102@units = "m/s"
lat2 = xlat2(:,0)
lon2 = xlon2(0,:)
it = 25
print("Working on time: " + times(it) )
;20180228
z = wrf_user_getvar(a, "z",28) ; grid point height
u = wrf_user_getvar(a,"ua",28) ; 3D U at mass points
v = wrf_user_getvar(a,"va",28) ; 3D V at mass points
w = wrf_user_getvar(a,"wa",28)
;20180401
z1 = wrf_user_getvar(a1, "z",it) ; grid point height
u1 = wrf_user_getvar(a1,"ua",it) ; 3D U at mass points
v1 = wrf_user_getvar(a1,"va",it) ; 3D V at mass points
w1 = wrf_user_getvar(a1,"wa",it)
;20180513
z2 = wrf_user_getvar(a2, "z",it) ; grid point height
u2 = wrf_user_getvar(a2,"ua",it) ; 3D U at mass points
v2 = wrf_user_getvar(a2,"va",it) ; 3D V at mass points
w2 = wrf_user_getvar(a2,"wa",it)
if ( FirstTime ) then ; get height info for labels
zmin = 0.
zmax = 15. ; We are only interested in the first 6km
nz = floattoint(zmax + 1)
end if
;plot AB垂直剖面的位置
plane = new(2,float)
plane = (/ mdims(nd-1)/2, mdims(nd-2)/2 /) ; pivot point is center of domain (x,y)
plane1 = new(2,float)
plane1 = (/ mdims1(nd1-1)/2, mdims1(nd1-2)/2 /) ; pivot point is center of domain (x,y)
plane2 = new(2,float)
plane2 = (/ mdims2(nd2-1)/2, mdims2(nd2-2)/2 /) ; pivot point is center of domain (x,y)
opts = False
angle = 0.
PI=3.1415926
alfa = (90-angle)*PI/180
X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
X_plane1 = wrf_user_intrp2d(xlat1,plane1,angle,opts)
X_plane2 = wrf_user_intrp2d(xlat2,plane2,angle,opts)
uv = u*sin(alfa)+v*cos(alfa)
uv1 = u1*sin(alfa)+v1*cos(alfa)
uv2 = u2*sin(alfa)+v2*cos(alfa)
;X_desc = "latitude"
;计算涡度和散度
u@_FillValue = 9.96921e+36
scale=1.e05
dv = new ( dimsizes(u), typeof(u), u@_FillValue )
dv@long_name = "divengence"
dv@units = "scaled"
dv = uv2dv_cfd (u,v,lat,lon, 2) * scale ;散度
dv1 = uv2dv_cfd (u1,v1,lat1,lon1, 2) * scale ;散度
dv2 = uv2dv_cfd (u2,v2,lat2,lon2, 2) * scale ;散度
copy_VarCoords(u, dv)
copy_VarCoords(u1, dv1)
copy_VarCoords(u2, dv2)
;把数据插值到各个高度上
;20180228
dv_plane = wrf_user_intrp3d(dv,z,"v",plane,angle,opts)
ter_plane = wrf_user_intrp2d(ter,plane,angle,opts)
u_plane = wrf_user_intrp3d(uv,z,"v",plane,angle,opts)
w_plane = wrf_user_intrp3d(w,z,"v",plane,angle,opts)
;20180401
dv_plane1 = wrf_user_intrp3d(dv1,z1,"v",plane1,angle,opts)
ter_plane1 = wrf_user_intrp2d(ter1,plane1,angle,opts)
u_plane1 = wrf_user_intrp3d(uv1,z1,"v",plane1,angle,opts)
w_plane1 = wrf_user_intrp3d(w1,z1,"v",plane1,angle,opts)
;20180513
dv_plane2 = wrf_user_intrp3d(dv2,z2,"v",plane2,angle,opts)
ter_plane2 = wrf_user_intrp2d(ter2,plane2,angle,opts)
u_plane2 = wrf_user_intrp3d(uv2,z2,"v",plane2,angle,opts)
w_plane2 = wrf_user_intrp3d(w2,z2,"v",plane2,angle,opts)
;20180228
if ( FirstTime ) then
zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts)
b = ind(zz(:,0) .gt. zmax*1000. )
zmax_pos = b(0) - 1
if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then
zspan = b(0) - 1
else
zspan = b(0)
end if
delete(zz)
delete(b)
FirstTime = False
end if
print(zmax_pos)
print(zspan)
;20180401
FirstTime = True
if ( FirstTime ) then
zz = wrf_user_intrp3d(z1,z1,"v",plane1,angle,opts)
b = ind(zz(:,0) .gt. zmax*1000. )
zmax_pos1 = b(0) - 1
if ( abs(zz(zmax_pos1,0)-zmax*1000.) .lt. abs(zz(zmax_pos1+1,0)-zmax*1000.) ) then
zspan1 = b(0) - 1
else
zspan1 = b(0)
end if
delete(zz)
delete(b)
FirstTime = False
end if
print(zmax_pos1)
print(zspan1)
;201805013
FirstTime = True
if ( FirstTime ) then
zz = wrf_user_intrp3d(z2,z2,"v",plane2,angle,opts)
b = ind(zz(:,0) .gt. zmax*1000. )
zmax_pos2 = b(0) - 1
if ( abs(zz(zmax_pos2,0)-zmax*1000.) .lt. abs(zz(zmax_pos2+1,0)-zmax*1000.) ) then
zspan2 = b(0) - 1
else
zspan2 = b(0)
end if
delete(zz)
delete(b)
FirstTime = False
end if
print(zmax_pos2)
print(zspan2)
; X-axis lables
dimsX = dimsizes(X_plane)
xmin = X_plane(0)
xmax = X_plane(dimsX(0)-1)
xspan = dimsX(0)-1
nx = floattoint( (xmax-xmin)/2 + 1)
nx = 10
;print("X_plane"+X_plane)
print("dimsX" + dimsX)
print("xmin and xmax in plot is " + xmin+" and "+xmax)
print("xspan and nx in plot is " + xspan+" and "+nx)
dimsX1 = dimsizes(X_plane1)
xmin1 = X_plane1(0)
xmax1 = X_plane1(dimsX1(0)-1)
xspan1 = dimsX1(0)-1
nx1 = floattoint( (xmax1-xmin1)/2 + 1)
nx1 = 10
;print("X_plane"+X_plane)
print("dimsX1" + dimsX)
print("xmin1 and xmax1 in plot is " + xmin1+" and "+xmax1)
print("xspan1 and nx1 in plot is " + xspan1+" and "+nx1)
dimsX2 = dimsizes(X_plane2)
xmin2 = X_plane2(0)
xmax2 = X_plane2(dimsX2(0)-1)
xspan2 = dimsX2(0)-1
nx2 = floattoint( (xmax2-xmin2)/2 + 1)
nx2 = 10
;print("X_plane"+X_plane)
print("dimsX2" + dimsX)
print("xmin2 and xmax2 in plot is " + xmin2+" and "+xmax2)
print("xspan2 and nx2 in plot is " + xspan2+" and "+nx2)
;---------------------------------------------------------------
;20180228
; Options for XY Plots
opts_xy = res
;opts_xy@tiXAxisString = X_desc
opts_xy@tiYAxisString = "Height (km)"
opts_xy@tmXTOn = False
opts_xy@tmYROn = False
;;-------------------想要绘制的数据的经纬度范围-------------------------------------------------
;;-----------------------------------------------------------------------------------------------
opts_xy@tmXBOn = True
opts_xy@tmYLOn = True
opts_xy@tmXBMode = "Explicit"
opts_xy@tmXBValues = fspan(0,xspan,nx) ; Create tick marks
;opts_xy@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nx)) ; Create labels
opts_xy@tmXBLabels = (/"41.9~S~o~N~N","42.2~S~o~N~N","42.7~S~o~N~N","43.2~S~o~N~N","43.7~S~o~N~N","44.2~S~o~N~N","44.7~S~o~N~N","45.2~S~o~N~N","45.7~S~o~N~N","46.2~S~o~N~N"/)
opts_xy@tmXBLabelFontHeightF = 0.010
opts_xy@tmYLMode = "Explicit"
opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks
opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels
opts_xy@tiXAxisFontHeightF = 0.020
opts_xy@tiYAxisFontHeightF = 0.020
opts_xy@tmXBMajorLengthF = 0.02
opts_xy@tmYLMajorLengthF = 0.02
opts_xy@tmYLLabelFontHeightF = 0.015
; Plotting options for terrain
opts_ter = opts_xy
opts_ter@gsnYRefLine = 0.0
opts_ter@gsnAboveYRefLineColor = "black"
opts_ter@gsnDraw = False
opts_ter@gsnFrame = False
opts_ter@trYMaxF = zmax*1000 ;地形数据为m，因此这里要放大倍数为km
;;plot for dv
opts_dv = opts_xy
opts_dv@gsnDraw = False
opts_dv@gsnFrame = False
opts_dv@cnFillOn = True
opts_dv@cnLinesOn = False
opts_dv@cnInfoLabelOn = False
opts_dv@cnLineLabelsOn = False
opts_dv@cnSmoothingOn = True
opts_dv@lbLabelBarOn = False
opts_dv@lbTitleOn = False
opts_dv@cnSmoothingDistanceF=0.015
opts_dv@cnSmoothingTensionF=0.001
opts_dv@cnLineLabelInterval=1
opts_dv@cnLineLabelPlacementMode="Computed"
opts_dv@cnLineLabelAngleF=0.5 ;大于0则平行于X轴，小于0为曲线切线方向
opts_dv@cnLineLabelDensityF=1.0
opts_dv@cnLineThicknessF = 4.0
opts_dv@cnInfoLabelOrthogonalPosF=0.00
opts_dv@ContourParameters=(/-10, 10, 2 /)
opts_dv@gsnContourPosLineDashPattern=0
opts_dv@cnLineDashPattern=0
vcres = opts_xy
vcres@vcGlyphStyle = "CurlyVector";"LineArrow";
vcres@vcPositionMode = "ArrowTail"
vcres@vcRefAnnoOn = True
vcres@vcMinDistanceF = 0.03
vcres@vcRefMagnitudeF = 8. ;单位长度
vcres@vcLineArrowThicknessF = 3.0
vcres@vcMapDirection = False
vcres@vcFillArrowsOn = True
vcres@vcVectorDrawOrder = "PreDraw"
vcres@vcRefLengthF = 0.015 ;矢量长度
vcres@vcFillArrowHeadXF = 0.2
vcres@vcFillArrowHeadYF = 0.2
vcres@vcFillArrowHeadInteriorXF = 0.5
;---------------------------------------------------------------
;20180401
; Options for XY Plots
opts_xy1 = res
;opts_xy1@tiXAxisString = X_desc
;opts_xy1@tiYAxisString = "Height (km)"
opts_xy1@tmXTOn = False
opts_xy1@tmYROn = False
;;-------------------想要绘制的数据的经纬度范围-------------------------------------------------
;;-----------------------------------------------------------------------------------------------
opts_xy1@tmXBOn = True
opts_xy1@tmYLOn = False
opts_xy1@tmXBMode = "Explicit"
opts_xy1@tmXBValues = fspan(0,xspan1,nx) ; Create tick marks
;opts_xy1@tmXBLabels = sprintf("%.1f",fspan(xmin1,xmax1,nx)) ; Create labels
opts_xy1@tmXBLabels = (/"41.9~S~o~N~N","42.2~S~o~N~N","42.7~S~o~N~N","43.2~S~o~N~N","43.7~S~o~N~N","44.2~S~o~N~N","44.7~S~o~N~N","45.2~S~o~N~N","45.7~S~o~N~N","46.2~S~o~N~N"/)
opts_xy1@tmXBLabelFontHeightF = 0.010
opts_xy1@tmYLMode = "Explicit"
opts_xy1@tmYLValues = fspan(0,zspan1,nz) ; Create tick marks
opts_xy1@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels
opts_xy1@tiXAxisFontHeightF = 0.020
opts_xy1@tiYAxisFontHeightF = 0.020
opts_xy1@tmXBMajorLengthF = 0.02
opts_xy1@tmYLMajorLengthF = 0.02
opts_xy1@tmYLLabelFontHeightF = 0.015
; Plotting options for terrain
opts_ter1 = opts_xy1
opts_ter1@gsnYRefLine = 0.0
opts_ter1@gsnAboveYRefLineColor = "black"
opts_ter1@gsnDraw = False
opts_ter1@gsnFrame = False
opts_ter1@trYMaxF = zmax*1000 ;地形数据为m，因此这里要放大倍数为km
;;plot for dv
opts_dv1 = opts_xy1
opts_dv1@gsnDraw = False
opts_dv1@gsnFrame = False
opts_dv1@cnFillOn = True
opts_dv1@cnLinesOn = False
opts_dv1@cnInfoLabelOn = False
opts_dv1@cnLineLabelsOn = False
opts_dv1@cnSmoothingOn = True
opts_dv1@lbLabelBarOn = False
opts_dv1@lbTitleOn = False
opts_dv1@cnSmoothingDistanceF=0.015
opts_dv1@cnSmoothingTensionF=0.001
opts_dv1@cnLineLabelInterval=1
opts_dv1@cnLineLabelPlacementMode="Computed"
opts_dv1@cnLineLabelAngleF=0.5 ;大于0则平行于X轴，小于0为曲线切线方向
opts_dv1@cnLineLabelDensityF=1.0
opts_dv1@cnLineThicknessF = 4.0
opts_dv1@cnInfoLabelOrthogonalPosF=0.00
opts_dv1@ContourParameters=(/-10, 10, 2 /)
opts_dv1@gsnContourPosLineDashPattern=0
opts_dv1@cnLineDashPattern=0
vcres1 = opts_xy1
vcres1@vcGlyphStyle = "CurlyVector";"LineArrow";
vcres1@vcPositionMode = "ArrowTail"
vcres1@vcRefAnnoOn = True
vcres1@vcMinDistanceF = 0.03
vcres1@vcRefMagnitudeF = 8. ;单位长度
vcres1@vcLineArrowThicknessF = 3.0
vcres1@vcMapDirection = False
vcres1@vcFillArrowsOn = True
vcres1@vcVectorDrawOrder = "PreDraw"
vcres1@vcRefLengthF = 0.015 ;矢量长度
vcres1@vcFillArrowHeadXF = 0.2
vcres1@vcFillArrowHeadYF = 0.2
vcres1@vcFillArrowHeadInteriorXF = 0.5
;---------------------------------------------------------------
;20180513
; Options for XY Plots
opts_xy2 = res
;opts_xy2@tiXAxisString = X_desc
;opts_xy2@tiYAxisString = "Height (km)"
opts_xy2@tmXTOn = False
opts_xy2@tmYROn = False
;;-------------------想要绘制的数据的经纬度范围-------------------------------------------------
;;-----------------------------------------------------------------------------------------------
opts_xy2@tmXBOn = True
opts_xy2@tmYLOn = False
opts_xy2@tmXBMode = "Explicit"
opts_xy2@tmXBValues = fspan(0,xspan2,nx) ; Create tick marks
;opts_xy2@tmXBLabels = sprintf("%.1f",fspan(xmin2,xmax2,nx)) ; Create labels
opts_xy2@tmXBLabels = (/"41.9~S~o~N~N","42.2~S~o~N~N","42.7~S~o~N~N","43.2~S~o~N~N","43.7~S~o~N~N","44.2~S~o~N~N","44.7~S~o~N~N","45.2~S~o~N~N","45.7~S~o~N~N","46.2~S~o~N~N"/)
opts_xy2@tmXBLabelFontHeightF = 0.010
opts_xy2@tmYLMode = "Explicit"
opts_xy2@tmYLValues = fspan(0,zspan2,nz) ; Create tick marks
opts_xy2@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels
opts_xy2@tiXAxisFontHeightF = 0.020
opts_xy2@tiYAxisFontHeightF = 0.020
opts_xy2@tmXBMajorLengthF = 0.02
opts_xy2@tmYLMajorLengthF = 0.02
opts_xy2@tmYLLabelFontHeightF = 0.015
; Plotting options for terrain
opts_ter2 = opts_xy2
opts_ter2@gsnYRefLine = 0.0
opts_ter2@gsnAboveYRefLineColor = "black"
opts_ter2@gsnDraw = False
opts_ter2@gsnFrame = False
opts_ter2@trYMaxF = zmax*1000 ;地形数据为m，因此这里要放大倍数为km
;;plot for dv
opts_dv2 = opts_xy2
opts_dv2@gsnDraw = False
opts_dv2@gsnFrame = False
opts_dv2@cnFillOn = True
opts_dv2@cnLinesOn = False
opts_dv2@cnInfoLabelOn = False
opts_dv2@cnLineLabelsOn = False
opts_dv2@cnSmoothingOn = True
opts_dv2@lbLabelBarOn = False
opts_dv2@lbTitleOn = False
opts_dv2@cnSmoothingDistanceF=0.015
opts_dv2@cnSmoothingTensionF=0.001
opts_dv2@cnLineLabelInterval=1
opts_dv2@cnLineLabelPlacementMode="Computed"
opts_dv2@cnLineLabelAngleF=0.5 ;大于0则平行于X轴，小于0为曲线切线方向
opts_dv2@cnLineLabelDensityF=1.0
opts_dv2@cnLineThicknessF = 4.0
opts_dv2@cnInfoLabelOrthogonalPosF=0.00
opts_dv2@ContourParameters=(/-10, 10, 2 /)
opts_dv2@gsnContourPosLineDashPattern=0
opts_dv2@cnLineDashPattern=0
vcres2 = opts_xy2
vcres2@vcGlyphStyle = "CurlyVector";"LineArrow";
vcres2@vcPositionMode = "ArrowTail"
vcres2@vcRefAnnoOn = True
vcres2@vcMinDistanceF = 0.03
vcres2@vcRefMagnitudeF = 8. ;单位长度
vcres2@vcLineArrowThicknessF = 3.0
vcres2@vcMapDirection = False
vcres2@vcFillArrowsOn = True
vcres2@vcVectorDrawOrder = "PreDraw"
vcres2@vcRefLengthF = 0.015 ;矢量长度
vcres2@vcFillArrowHeadXF = 0.2
vcres2@vcFillArrowHeadYF = 0.2
vcres2@vcFillArrowHeadInteriorXF = 0.5
countour_dv = wrf_contour(a, wks, dv_plane(0:zmax_pos,:), opts_dv)
countour_dv1 = wrf_contour(a1, wks, dv_plane1(0:zmax_pos1,:), opts_dv1)
countour_dv2 = wrf_contour(a2, wks, dv_plane2(0:zmax_pos2,:), opts_dv2)
contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter)
contour_ter1 = gsn_csm_xy(wks,X_plane1,ter_plane1,opts_ter1)
contour_ter2 = gsn_csm_xy(wks,X_plane2,ter_plane2,opts_ter2)
vector = gsn_vector(wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:)*10.,vcres)
vector1 = gsn_vector(wks,u_plane1(0:zmax_pos1,:),w_plane1(0:zmax_pos1,:)*10.,vcres1)
vector2 = gsn_vector(wks,u_plane2(0:zmax_pos2,:),w_plane2(0:zmax_pos2,:)*10.,vcres2)
plots(0) = wrf_overlays(a,wks,(/countour_dv,vector,contour_ter/),True) ; plot x-section
plots(1) = wrf_overlays(a1,wks,(/countour_dv1,vector1,contour_ter/),True) ; plot x-section
plots(2) = wrf_overlays(a2,wks,(/countour_dv2,vector2,contour_ter2/),True) ; plot x-section
; Delete options and fields, so we don't have carry over
delete(opts_xy)
delete(ter_plane)
delete(opts_xy1)
delete(ter_plane1)
delete(opts_xy2)
delete(ter_plane2)
;set panel
pnlres = res
pnlres@gsnPanelCenter = True
pnlres@gsnPanelLabelBar = True
pnlres@lbOrientation = "Vertical"
pnlres@lbLabelPosition = "Right"
pnlres@pmLabelBarParallelPosF = 0.0
pnlres@pmLabelBarOrthogonalPosF = 0.0
pnlres@pmLabelBarWidthF = 0.05
pnlres@pmLabelBarHeightF = 0.15
pnlres@gsnPanelRowSpec = True
pnlres@gsnPanelYWhiteSpacePercent = 5
pnlres@gsnPanelScalePlotIndex = 1
pnlres@amJust = "TopLeft"
gsn_panel(wks, (/plots/), (/3/), pnlres)
end
有点长 非常感谢！！ |
|
窥视卡
雷达卡
发表于 2019-8-31 20:03:03
显身卡