I have modified one of the cross section example scripts to display a
WRF cross section of theta along a line between two lat/lon points. The
vertical coordinate is height in km. Currently the script displays the
entire vertical depth of all WRF levels, but I'd like to display only
the lowest 4 km or so. What changes do I need to make to do that?
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Script to plot WRF forecast theta on a cross section
; This script will plot data from a a given point A to point B
; Vertical coordinate is height
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
a =
addfile("/media/disk-1/taku_00011612/wrfout_d03_2000-01-16_21:00:00.nc","r")
; We generate plots, but what kind do we prefer?
type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
wks = gsn_open_wks(type,"CrossSection2")
gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map
; Set some Basic Plot options
ARWres = True
ARWres@MainTitle = "REAL-TIME WRF"
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; What times and how many time steps are in the data set?
FirstTime = True
times = wrf_user_list_times(a) ; get times in the file
ntimes = dimsizes(times) ; number of times in the file
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
do it = 0,ntimes-1,2 ; TIME LOOP
print("Working on time: " + times(it) )
ARWres@TimeLabel = times(it) ; Set Valid time to use on plots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need
w = wrf_user_getvar(a,"wa",it) ; T in C
th = wrf_user_getvar(a,"th",it) ; relative humidity
z = wrf_user_getvar(a, "z",it) ; grid point height
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if (FirstTime) then
; THIS IS NEEDED FOR Y LABLES - ALWAYS DO
zz = z/100.
dims = dimsizes(zz)
zmax = -999.
do imax = 0,dims(0)-1
if ( .not.ismissing(zz(imax,1,1)) .and. zz(imax,1,1) .gt. zmax
) then
zmax = zz(imax,1,1)
end if
end do
imax = floattoint(zmax)
zmax = imax/10.
print("Y labels set to: 0 - " + zmax + " km")
FirstTime = False
; END OF ALWAYS DO
end if
do ip = 1, 1
opts = True
plane = new(4,float)
if(ip .eq. 1) then ; Specify the start/end points of the
cross section in lat/lon
lats = (/58.20, 58.38/) ;
lons = (/-134.46, -134.18/) ;
loc = wrf_user_ll_to_ij(a, lons, lats, True)
; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y)
; subtract one since we want to use it as an index in NCL
x_start = loc(0,0) - 1
x_end = loc(0,1) - 1
y_start = loc(1,0) - 1
y_end = loc(1,1) - 1
plane = (/ x_start,y_start , x_end,y_end /) ; approx. start
x;y and end x;y point
end if
if(ip .eq. 2) then
plane = (/ 130,1 , 130,162 /) ; approx. start x;y and end
x;y point
end if
if(ip .eq. 3) then
plane = (/ 49,1 , 210,162 /) ; approx. start x;y and end
x;y point
end if
; Interpolate data vertically (in z)
th_plane = wrf_user_intrp3d(th,z,"v",plane,0.,opts)
w_plane = wrf_user_intrp3d(w,z,"v",plane,0.,opts)
; Options for XY Plots
opts_xy = ARWres
opts_xy@tiXAxisString = "Number of Grid Points"
opts_xy@tiYAxisString = "Height (km)"
opts_xy@AspectRatio = .75
opts_xy@cnMissingValPerimOn = True
opts_xy@cnMissingValFillColor = 0
opts_xy@cnMissingValFillPattern = 11
opts_xy@tmYLMode = "Explicit"
opts_xy@tmYLValues = ispan(0,100,10) ; over 0-100
span of Y label
opts_xy@tmYLLabels = fspan(0,zmax,11);
Corresponding 11 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 theta
opts_th = opts_xy
opts_th@cnLevelSelectionMode = "ManualLevels" ;set manual
contour lvls
opts_th@cnMinLevelValF = 260. ; set min
contour level
opts_th@cnMaxLevelValF = 290. ; set max
contour level
opts_th@cnLevelSpacingF = 1 ; set contour spacing
opts_th@cnFillOn = True
opts_th@gsnSpreadColors = True
; Plotting options for vertical vel
opts_w = opts_xy
opts_w@ContourParameters = (/ 1. /)
; Get the contour info for the th and temp
contour_th = wrf_contour(a,wks,th_plane,opts_th)
contour_w = wrf_contour(a,wks,w_plane,opts_w)
; MAKE PLOTS
wrf_overlay(wks,(/contour_th,contour_w/),True)
; Delete options and fields, so we don't have carry over
delete(opts_w)
delete(opts_th)
delete(w_plane)
delete(th_plane)
end do ; make next cross section
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end do ; END OF TIME LOOP
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Apr 25 18:22:41 2011
This archive was generated by hypermail 2.1.8 : Tue May 03 2011 - 14:47:35 MDT