Re: NCL not plotting first 200 m in the vertical cross section

From: Mary Haley <haley_at_nyahnyahspammersnyahnyah>
Date: Mon Feb 14 2011 - 08:17:41 MST

This might be due to a bug we fixed awhile back. What version of NCL (ncl -V) are you running?
If you can provide me with the file, then I can try your script and see if the bug is fixed.

--Mary

On Feb 13, 2011, at 3:15 PM, Basit Khan wrote:

> 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_2006-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.
>
> <vcs_ncl.eps>_______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Feb 14 08:17:51 2011

This archive was generated by hypermail 2.1.8 : Tue Feb 15 2011 - 09:43:19 MST