<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
</head>
<body>
Dear WRF users and NCL guys,<br>
<br>
when I use wrf_user_intrp3d in NCL to make a vertical interpolation (in meter), I found the first interpolating level is always missing value. For example my model top is 30 km, then the first interpolating level is 300 meter, and it is always a missing value.
As you can from the figure-Z_Vertical_CrossSection.000004.png, the white part in the surface represents the missing value which contaminates the true terrain height.<br>
<br>
I will be very appreciated if anyone give me some suggestion.<br>
<br>
I copy the script as follows.<br>
<br>
Best regards,<br>
Guiting<br>
<br>
<br>
<br>
<br>
<br>
<br>
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"<br>
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"<br>
;load "./WRFUserARW.ncl"<br>
<br>
begin<br>
;<br>
; The WRF ARW input file. <br>
; This needs to have a ".nc" appended, so just do it.<br>
a = addfile("./<a href="http://wrfout.nc">wrfout.nc</a>","r")<br>
<br>
<br>
; We generate plots, but what kind do we prefer?<br>
type = "x11"<br>
; type = "pdf"<br>
; type = "ps"<br>
; type = "ncgm"<br>
type = "png"<br>
wks = gsn_open_wks(type,"Z_Vertical_CrossSection")<br>
<br>
<br>
; Set some basic resources<br>
res = True<br>
res@MainTitle = "REAL-TIME WRF"<br>
res@Footer = False<br>
<br>
pltres = True<br>
<br>
<br>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<br>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<br>
<br>
; What times and how many time steps are in the data set?<br>
FirstTime = True<br>
FirstTimeMap = True<br>
times = wrf_user_list_times(a) ; get times in the file<br>
ntimes = dimsizes(times) ; number of times in the file<br>
<br>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<br>
<br>
xlat = wrf_user_getvar(a, "XLAT",0)<br>
xlon = wrf_user_getvar(a, "XLONG",0)<br>
ter = wrf_user_getvar(a, "HGT",0)<br>
<br>
do it = 0,ntimes-1,2 ; TIME LOOP<br>
<br>
print("Working on time: " + times(it) )<br>
res@TimeLabel = times(it) ; Set Valid time to use on plots<br>
<br>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<br>
; First get the variables we will need <br>
<br>
tc = wrf_user_getvar(a,"tc",it) ; T in C<br>
rh = wrf_user_getvar(a,"rh",it) ; relative humidity<br>
z = wrf_user_getvar(a, "z",it) ; grid point height<br>
<br>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<br>
<br>
do ip = 1, 3 ; we are doing 3 plots<br>
; all with the pivot point (plane) in the center of the domain<br>
; at angles 0, 45 and 90<br>
; |<br>
; |<br>
; angle=0 is |, angle=90 is ------<br>
; |<br>
; |<br>
<br>
; Build plane (pivot point) through which the cross section will go <br>
; OR set to zero, if start and end points are specified<br>
; IF using plane and angle, set opts in wrf_user_intrp3d to False<br>
dimsrh = dimsizes(rh)<br>
plane = new(2,float)<br>
plane = (/ dimsrh(2)/2, dimsrh(1)/2 /) ; pivot point is center of domain<br>
opts = False<br>
<br>
if(ip .eq. 1) then<br>
angle = 90.<br>
X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)<br>
X_desc = "longitude"<br>
end if<br>
if(ip .eq. 2) then<br>
angle = 0.<br>
X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)<br>
X_desc = "latitude"<br>
end if<br>
if(ip .eq. 3) then<br>
angle = 45.<br>
X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)<br>
X_desc = "longitude"<br>
end if<br>
<br>
; X-axis lables<br>
dimsX = dimsizes(X_plane)<br>
xmin = X_plane(0)<br>
xmax = X_plane(dimsX(0)-1)<br>
xspan = dimsX(0)-1<br>
nxlabs = floattoint( (xmax-xmin)/2 + 1)<br>
<br>
if (FirstTimeMap) then<br>
lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)<br>
lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)<br>
mpres = True<br>
mpres@mpGeophysicalLineColor = "Black"<br>
mpres@mpNationalLineColor = "Black"<br>
mpres@mpUSStateLineColor = "Black"<br>
mpres@mpGridLineColor = "Black"<br>
mpres@mpLimbLineColor = "Black"<br>
mpres@mpPerimLineColor = "Black"<br>
pltres = True<br>
pltres@FramePlot = False<br>
optsM = res<br>
optsM@NoHeaderFooter = True<br>
optsM@cnFillOn = True<br>
optsM@lbTitleOn = False<br>
contour = wrf_contour(a,wks,ter,optsM)<br>
plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)<br>
lnres = True<br>
lnres@gsLineThicknessF = 3.0<br>
lnres@gsLineColor = "black";"Red"<br>
do ii = 0,dimsX(0)-2<br>
gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres)<br>
end do<br>
frame(wks)<br>
delete(lon_plane)<br>
delete(lat_plane)<br>
pltres@FramePlot = True<br>
end if<br>
<br>
<br>
if (FirstTime) then<br>
; THIS IS NEEDED FOR LABLES - ALWAYS DO (Z axis only needed once. X everytime plane changes)<br>
<br>
; Y-axis labels<br>
zmax = 6000. ; We only want to see the first 6 km<br>
zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts)<br>
dims = dimsizes(zz)<br>
do imax = 0,dims(0)-1<br>
if ( .not.ismissing(zz(imax,0)) .and. zz(imax,0) .lt. zmax ) then<br>
zmax_pos = imax<br>
end if<br>
end do<br>
zspan = zmax_pos<br>
zmin = z(0,0,0)<br>
zmax = zz(zmax_pos,0)<br>
print(zmax)<br>
zmin=zmin/1000.<br>
zmax=zmax/1000. <br>
nzlabs = floattoint(zmax + 1)<br>
<br>
FirstTime = False<br>
; END OF ALWAYS DO<br>
end if<br>
<br>
; Interpolate data vertically (in z)<br>
<br>
rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts)<br>
tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts)<br>
<br>
; Options for XY Plots<br>
opts_xy = res<br>
opts_xy@tiXAxisString = X_desc<br>
opts_xy@tiYAxisString = "Height (km)"<br>
opts_xy@cnMissingValPerimOn = True<br>
opts_xy@cnMissingValFillColor = 0<br>
opts_xy@cnMissingValFillPattern = 11<br>
opts_xy@tmXTOn = False<br>
opts_xy@tmYROn = False<br>
opts_xy@tmXBMode = "Explicit"<br>
opts_xy@tmXBValues = fspan(0,xspan,nxlabs) ; Create the correct tick marks<br>
opts_xy@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nxlabs)); Create labels<br>
opts_xy@tmXBLabelFontHeightF = 0.015<br>
opts_xy@tmYLMode = "Explicit"<br>
opts_xy@tmYLValues = fspan(0,zspan,nzlabs) ; Create the correct tick marks<br>
opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nzlabs)); Create labels<br>
opts_xy@tiXAxisFontHeightF = 0.020<br>
opts_xy@tiYAxisFontHeightF = 0.020<br>
opts_xy@tmXBMajorLengthF = 0.02<br>
opts_xy@tmYLMajorLengthF = 0.02<br>
opts_xy@tmYLLabelFontHeightF = 0.015<br>
opts_xy@PlotOrientation = tc_plane@Orientation<br>
<br>
<br>
; Plotting options for RH<br>
opts_rh = opts_xy<br>
opts_rh@ContourParameters = (/ 10., 90., 10. /)<br>
opts_rh@pmLabelBarOrthogonalPosF = -0.1 <br>
opts_rh@cnFillOn = True<br>
opts_rh@cnFillColors = (/"White","White","White", \<br>
"White","Chartreuse","Green", \<br>
"Green3","Green4", \<br>
"ForestGreen","PaleGreen4"/)<br>
<br>
; Plotting options for Temperature<br>
opts_tc = opts_xy<br>
opts_tc@cnInfoLabelZone = 1<br>
opts_tc@cnInfoLabelSide = "Top"<br>
opts_tc@cnInfoLabelPerimOn = True<br>
opts_tc@cnInfoLabelOrthogonalPosF = -0.00005<br>
opts_tc@ContourParameters = (/ 5. /)<br>
<br>
<br>
; Get the contour info for the rh and temp<br>
contour_tc = wrf_contour(a,wks,tc_plane(0:zmax_pos,:),opts_tc)<br>
contour_rh = wrf_contour(a,wks,rh_plane(0:zmax_pos,:),opts_rh)<br>
<br>
<br>
; MAKE PLOTS <br>
plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres)<br>
<br>
; Delete options and fields, so we don't have carry over<br>
delete(opts_xy)<br>
delete(opts_tc)<br>
delete(opts_rh)<br>
delete(tc_plane)<br>
delete(rh_plane)<br>
delete(X_plane)<br>
<br>
end do ; make next cross section<br>
<br>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<br>
<br>
FirstTimeMap = False<br>
end do ; END OF TIME LOOP<br>
<br>
end<br>
<br>
<hr>
<font face="Arial" color="Gray" size="2">CONFIDENTIALITY: This email is intended solely for the person(s) named and may be confidential and/or privileged. If you are not the intended recipient, please delete it, notify us and do not copy, use, or disclose its
content. Thank you.<br>
<br>
Towards A Sustainable Earth: Print Only When Necessary<br>
</font>
</body>
</html>