- 积分
- 3751
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-10-15
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2016-11-17 21:18:57
|
显示全部楼层
脚本如下
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
a = addfile("/wrfout_24_00:00:00.nc","r")
; We generate plots, but what kind do we prefer?
type = "png"
; type = "pdf"
; type = "ps"
; type = "ncgm"
wks = gsn_open_wks(type,"dbz115.95-12km")
gsn_define_colormap(wks,"radar_1")
; Cross section
minlat= 30.5
maxlat= 32.0
minlon = 115.95
maxlon = 115.95
loc = wrf_user_ll_to_ij(a,(/minlon,maxlon/),(/minlat,maxlat/),True)
FirstTime = True
opts = True
angle= 0.
plane = new(4,float)
plane = (/ loc(0,0)-1, loc(0,1)-1, loc(1,0)-1,loc(1,1)-1 /)
; Set some basic resources
res = True
res@MainTitle = ""
pltres = True
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FirstTime = 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)
;---------------------------------------------------------------
xlon = wrf_user_getvar(a, "XLONG",0)
xlat = wrf_user_getvar(a,"XLAT",0)
do it = 31,31 ; TIME LOOP
print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots
; u = wrf_user_getvar(a,"ua",it) ; theta
; w = wrf_user_getvar(a,"wa",it) ; relative humidity
z = wrf_user_getvar(a, "z",it) ; grid point height
p = wrf_user_getvar(a, "pressure",it) ; grid point height
dbz= wrf_user_getvar(a,"dbz",it)
ter = wrf_user_getvar(a, "HGT",0)
; u_plane = wrf_user_intrp3d(u,z,"v",plane,0.,opts)
; w_plane = wrf_user_intrp3d(w,z,"v",plane,0.,opts)
dbz_plane= wrf_user_intrp3d(dbz,z,"v",plane,0.,opts)
; Let's create nice labels - only have to do this once
if ( FirstTime ) then
zmin = 0.
zmax = 12.
nz = floattoint(zmax+1)
end if
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_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
; X_desc = "longitude"
X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
X_desc = "latitude"
dimsX = dimsizes(X_plane)
printVarSummary(X_plane)
xmin = X_plane(0)
xmax = X_plane(dimsX-1)
xspan = dimsX(0)-1
nx = floattoint( (xmax-xmin)/0.5 + 1)
; Options for XY Plots
opts_xy = res
opts_xy@tiYAxisString = "Height (km)"
opts_xy@AspectRatio = 0.75
opts_xy@cnMissingValPerimOn = True
opts_xy@cnMissingValFillColor = 0
opts_xy@cnMissingValFillPattern = 11
opts_xy@tmXBMode = "Explicit"
opts_xy@tmXBValues = fspan(0,xspan,nx) ; Create tick marks
opts_xy@tmXBLabels = sprintf("%.02f",fspan(xmin,xmax,nx)) ; Create labels
; opts_xy@trXReverse = True
opts_xy@tmYLMode = "Explicit"
opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks
opts_xy@tmYLLabels = sprintf("%.02f",fspan(zmin,zmax,nz)) ; Create labels
opts_xy@tiXAxisFontHeightF = 0.050
opts_xy@tiYAxisFontHeightF = 0.050
opts_xy@tmXBMajorLengthF = 0.05
opts_xy@tmYLMajorLengthF = 0.05
opts_xy@tmYLLabelFontHeightF = 0.05
opts_xy@PlotOrientation = dbz_plane@Orientation
opts_dbz = opts_xy
opts_dbz@pmLabelBarOrthogonalPosF = -0.06
opts_dbz@cnFillOn = True
opts_dbz@cnLevelSelectionMode = "ExplicitLevels"
opts_dbz@cnLevels = (/10,15,20,25,30,35,40,45,50,55,60,65,70/)
opts_dbz@cnFillColors = (/ 0,11,10,13,16,17,18,19,20,21,22,23,4,24/)
; opts_dbz@lbLabelFontHeightF =0.035
; opts_dbz@lbTitleFontHeightF =0.035
contour_dbz = wrf_contour(a,wks,dbz_plane(0:zmax_pos,:),opts_dbz)
; MAKE PLOTS
plot = wrf_overlays(a,wks,(/contour_dbz/),pltres)
; Delete options and fields, so we don't have carry over
delete(opts_dbz)
delete(dbz_plane)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end do ; END OF TIME LOOP
end
|
|