Hi,
I am trying to plot vertical cross section (height x latitude) from a
WRF output file. The problem is that NCL not plotting anything in the
first 200 m while I have at least 11 vertical levels in the first 200 m
asl. I am using terrain following coordinates system and most of the
cross section area is either sea or terrain less than 50 m high asl. I
will greatly appreciate if somebody could help me fix this problem. The
ncl code is copied below and the output of the ncl code is attached.
Regards,
Basit
Centre for Atmospheric Research
Department of Geography,
University of Canterbury, Private Bag 4800,
Christchurch 8140, New Zealand.
Tel: +643 364 2987 ext. 3088
email: basit.khan@canterbury.ac.nz
------------------------------------------------------------------------
-------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
;
a =
addfile("/hpc/bluefern/WRF/basit/basit-archive/b-haze/bh30/wrfout_d03_20
06-05-30_12: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,"plt_CrossSection4")
; Set some basic resources
res = True
res@MainTitle = "REAL-TIME WRF"
res@Footer = False
pltres = True
;;;;;;;;;;;;;;;; HIGHRES RANGS/GSHHS COASTLINE PLOT RESOURCES
;;;;;;;;;;;;;;;;;;;;;;;;
mpres =True
mpres@gsnMaximize = True
mpres@mpOutlineOn = True
mpres@mpFillOn = False
mpres@mpProjection = "Mercator"
mpres@mpLimitMode = "LatLon"
mpres@mpMinLonF = 174.4061
mpres@mpMaxLonF = 175.2502
mpres@mpMinLatF = -37.27655
mpres@mpMaxLatF = -36.60195
mpres@mpOutlineOn = True
mpres@mpGridLineColor = -1
mpres@gsnDraw = False ; don't draw the plots
mpres@gsnFrame = False ; don't advance the
frame
mpres@tmYROn = False
mpres@tmXTOn = False
mpres@mpGeophysicalLineColor= "Red"
mpres@mpLimbLineColor = "Blue"
mpres@mpGeophysicalLineThicknessF = 1.5
mpres@mpDataBaseVersion = "HighRes" ; Use the high-res database
mpres@pmTickMarkDisplayMode = "Always" ; Turn on map tickmarks.
; mpres@tiMainString = "Using the RANGS-GSHHS coastal
database"
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FirstTime = True
FirstTimeMap = True
times = wrf_user_list_times(a) ; 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; print(mdims)=97x50x150x150
nd = dimsizes(mdims) ; incase of 'P' the print(nd)=4
cause P has four dims
xlat = wrf_user_getvar(a, "XLAT",0)
xlon = wrf_user_getvar(a, "XLONG",0)
ter = wrf_user_getvar(a, "HGT",0)
;---------------------------------------------------------------
do it = 58,74,1 ; TIME LOOP
print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots
tc = wrf_user_getvar(a,"tc",it) ; T in C
rh = wrf_user_getvar(a,"rh",it) ; relative humidity
z = wrf_user_getvar(a, "height",it) ; grid point height
if ( FirstTime ) then ; get height info for labels
zmin = 0
zmax = 1. ; We are only interested in the
first 1000 m
nz = floattoint(zmax + 10)
end if
;---------------------------------------------------------------
plane = new(2,float)
plane = (/ mdims(nd-1)/2, mdims(nd-2)/2 /) ; pivot point is
center of domain (x,y)
opts = False
angle = 0.
X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
X_desc = "latitude"
rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts) ;
"v" is for vertical and "h" is for horinzontal interporlation
tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts) ;
Array of vertical coordninates to interpolate
; 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)
; printVarSummary(zz) = dimsizesV x Dimsize H = 100 x 149
b = ind(zz(:,0) .gt. zmax*1000. )
; (zmax *1000 --> conver km into metres
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 + 10)
; nx = (xmax-xmin)/2
;---------------------------------------------------------------
; print("the xmin value is "+xmin)
; print("the xmax value is "+xmax)
; 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,nx)
; Create tick marks
; opts_xy@tmXBValues = fspan(0,xspan,nx)
; Create tick marks
opts_xy@tmXBLabels =
sprintf("%6.2f",fspan(xmin,xmax,nx)) ; Create labels
opts_xy@tmXBLabelFontHeightF = 0.015
opts_xy@tmYLMode = "Explicit"
opts_xy@tmYLValues = fspan(0,zspan,nz)
; Create tick marks
opts_xy@tmYLLabels =
sprintf("%3.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 = tc_plane@Orientation
; Plotting options for RH
opts_rh = opts_xy
opts_rh@ContourParameters = (/ 10., 90., 10. /)
opts_rh@pmLabelBarOrthogonalPosF = -0.1
opts_rh@cnFillOn = True
opts_rh@cnFillColors = (/"White","White","White", \
"White","Chartreuse","Green", \
"Green3","Green4", \
"ForestGreen","PaleGreen4"/)
; Plotting options for Temperature
opts_tc = opts_xy
opts_tc@cnInfoLabelZone = 1
opts_tc@cnInfoLabelSide = "Top"
opts_tc@cnInfoLabelPerimOn = True
opts_tc@cnInfoLabelOrthogonalPosF = -0.00005
opts_tc@ContourParameters = (/ 1. /)
; Get the contour info for the rh and temp
contour_tc = wrf_contour(a,wks,tc_plane(0:zmax_pos,:),opts_tc)
contour_rh = wrf_contour(a,wks,rh_plane(0:zmax_pos,:),opts_rh)
;---------------------------------------------------------------
; 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
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"
do ii = 0,dimsX(0)-2
gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),
lat_plane(ii+1)/),lnres)
end do
frame(wks)
delete(lon_plane)
delete(lat_plane)
pltres@FramePlot = True
end if
plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres) ;
plot x-section
; Delete options and fields, so we don't have carry over
delete(opts_xy)
delete(opts_tc)
delete(opts_rh)
delete(tc_plane)
delete(rh_plane)
delete(X_plane)
end do ; make next cross section
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FirstTimeMap = False
; end do ; END OF TIME LOOP
end
This email may be confidential and subject to legal privilege, it may
not reflect the views of the University of Canterbury, and it is not
guaranteed to be virus free. If you are not an intended recipient,
please notify the sender immediately and erase all copies of the message
and any attachments.
Please refer to http://www.canterbury.ac.nz/emaildisclaimer for more
information.
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
This archive was generated by hypermail 2.1.8 : Tue Feb 15 2011 - 09:43:19 MST