- 积分
- 2139
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-4-15
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2018-3-23 13:08:23
|
显示全部楼层
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/wrf/WRFUserARW.ncl"
begin
a = addfile("~/urban_slope/wrfout_d04_2015-12-29_00_00_00.nc","r")
wks = gsn_open_wks("png","uw_theta")
; Set some basic resources
res = True
res@Footer = False
pltres = True
ter_res = True
opts_ter = ter_res
opts_ter@gsnYRefLine = 0.0
opts_ter@gsnAboveYRefLineColor = "black"
opts_ter@gsnDraw = False
opts_ter@gsnFrame = False
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FirstTime = True
FirstTimeMap = True
times = wrf_user_getvar(a,"times",-1) ; get times in the file
ntimes = dimsizes(times) ; number of times in the file
mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
nd = dimsizes(mdims)
xlat = wrf_user_getvar(a, "XLAT",0)
xlon = wrf_user_getvar(a, "XLONG",0)
ter = wrf_user_getvar(a, "HGT",0)
;---------------------------------------------------------------
do it = 16,39,1 ; TIME LOOP
; it = 31
print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots
u = wrf_user_getvar(a,"ua",it)
; v = wrf_user_getvar(a,"va",it)
w = wrf_user_getvar(a,"wa",it)
theta = wrf_user_getvar(a,"theta",it)
z = wrf_user_getvar(a, "z",it)
if ( FirstTime ) then ; get height info for labels
zmin = 0.
zmax = 6. ; We are only interested in the first 6km
nz = floattoint(zmax + 1)
end if
opts_ter@trYMaxF = zmax*1000
;---------------------------------------------------------------
do ip = 1, 3 ; we are doing 3 plots
plane = new(2,float)
plane = (/106,104/) ; pivot point is center of domain (x,y)
opts = False
if(ip .eq. 1) then
angle = 90.
X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
X_desc = "longitude"
end if
if(ip .eq. 2) then
angle = 0.
X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
X_desc = "latitude"
end if
if(ip .eq. 3) then
angle = 45.
X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
X_desc = "longitude"
end if
theta_plane = wrf_user_intrp3d(theta,z,"v",plane,angle,opts)
u_plane = wrf_user_intrp3d(u,z,"v",plane,angle,opts)
; v_plane = wrf_user_intrp3d(v,z,"v",plane,angle,opts)
w_plane = wrf_user_intrp3d(w,z,"v",plane,angle,opts)
ter_plane = wrf_user_intrp2d(ter,plane,angle,opts)
print("Max terrain height in plot " + max(ter_plane))
theta_plane2 = theta_plane
w_plane2 = w_plane
cross_dims = dimsizes(theta_plane2)
rank = dimsizes(cross_dims)
iz_do = 25
do iz = 0,24
iz_do = iz_do-1
do ix = 0,cross_dims(rank-1)-1
if ( ismissing(theta_plane2(iz_do,ix)) ) then
theta_plane2(iz_do,ix) = theta_plane2(iz_do+1,ix)
end if
if ( ismissing(w_plane2(iz_do,ix)) ) then
w_plane2(iz_do,ix) = w_plane2(iz_do+1,ix)
end if
end do
end do
; Find the index where 6km is - only need to do this once
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
; 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)
;---------------------------------------------------------------
; Options for XY Plots
opts_xy = res
opts_xy@tiXAxisString = X_desc
opts_xy@tiYAxisString = "Height (km)"
opts_xy@cnMissingValPerimOn = True
opts_xy@cnMissingValFillColor = 0
opts_xy@cnMissingValFillPattern = 11
opts_xy@tmXTOn = False
opts_xy@tmYROn = False
opts_xy@tmXBMode = "Explicit"
opts_xy@tmXBValues = fspan(0,xspan,6) ; Create tick marks
opts_xy@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,6))
opts_xy@tmXBLabelFontHeightF = 0.015
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
opts_xy@PlotOrientation = w_plane@Orientation
; Plotting options for RH
opts_theta = opts_xy
opts_theta@pmLabelBarOrthogonalPosF = -0.18
opts_theta@cnFillOn = True
opts_theta@cnLevelSelectionMode = "ExplicitLevels"
opts_theta@cnLevels = (/286,287,288,289,290,292,294,296,300,304,308,312/)
; Plotting options for Temperature
;opts_vec = opts_xy
vcres = opts_xy
vcres@vcGlyphStyle = "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@vcRefLengthF = 0.02
vcres@vcFillArrowHeadXF = 0.2
vcres@vcFillArrowHeadYF = 0.2
vcres@vcFillArrowHeadInteriorXF = 0.15
;opts_vec@cnInfoLabelZone = 1
;opts_vec@cnInfoLabelSide = "Top"
;opts_vec@cnInfoLabelPerimOn = True
;opts_vec@cnInfoLabelOrthogonalPosF = -0.00005
;opts_vec@ContourParameters = (/ 5. /)
;Contour terrain cross section
contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter)
; Get the contour info for the rh and temp
vector = wrf_vector(a,wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:)*20.,vcres)
contour_theta = wrf_contour(a,wks,theta_plane(0:zmax_pos,:),opts_theta)
vector2 = wrf_vector(a,wks,u_plane(0:zmax_pos,:),w_plane2(0:zmax_pos,:)*20.,vcres)
contour_theta2 = wrf_contour(a,wks,theta_plane2(0:zmax_pos,:),opts_theta)
;---------------------------------------------------------------
; MAKE PLOTS
if (FirstTimeMap) then
lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
mpres = True
pltres = True
pltres@FramePlot = False
optsM = res
optsM@NoHeaderFooter = True
optsM@cnFillOn = True
optsM@lbTitleOn = False
optsM@cnLevelSelectionMode = "ExplicitLevels"
optsM@cnLevels = (/0,2,100,300,500,700,1000,1500,2000,2500,3000,3500,4000,4500,5000,5600/)
contour = wrf_contour(a,wks,ter,optsM)
plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
lnres = True
lnres@gsLineThicknessF = 3.0
lnres@gsLineColor = "Red"
gsn_polyline(wks,plot,(/lon_plane(0),lon_plane(217)/),(/30.67,30.67/),lnres)
; gsn_polyline(wks,plot,(/104.41,104.41/),(/lat_plane(0),lat_plane(211)/),lnres)
frame(wks)
delete(lon_plane)
delete(lat_plane)
pltres@FramePlot = True
end if
;plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres)
;plot = wrf_overlays(a,wks,(/contour_rh,contour_tc,contour_ter/),pltres)
plot = wrf_overlays(a,wks,(/contour_theta2,contour_ter,vector/),pltres)
; Delete options and fields, so we don't have carry over
delete(opts_xy)
delete(vcres)
delete(opts_theta)
delete(w_plane)
delete(theta_plane)
delete(w_plane2)
delete(theta_plane2)
delete(X_plane)
delete(ter_plane)
; end do ; make next cross section
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FirstTimeMap = False
end do ; END OF TIME LOOP
end
完全是按照之前的网址改的，具体还要自己调整，希望能帮助到你
|
|
窥视卡
雷达卡
发表于 2018-3-23 11:08:33
显身卡