Re: NCL cross section with wind barbs

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Wed Apr 21 2010 - 13:52:17 MDT

Hello,

This is being forwarded to ' wrfhelp@ucar.edu '

Good luck

On 04/21/2010 08:20 AM, Luis C. Cana Cascallar wrote:
> Hi all,
> I'm trying to make a vertical cross-section plot with wind barbs using
> NCL. When the plot is done, the wind barbs all have a southerly
> orientation. Can anyone telling me what I'm doing wrong? Thanks in advance,
>
> Luis Cana
>
> Dr. Luis C. Cana-Cascallar
> Departamento de Física
> ULPGC
> 35017-Las Palmas de Gran Canaria
>
> lcana@dfis.ulpgc.es <mailto:lcana@dfis.ulpgc.es>
>
> Here is the script I am using:
> ========================================
>
> ; Plot data on a cross section (Modified from plt_CrossSection4.ncl)
> ; This script will plot data at a set angle through a specified point
> ; This script adds lon/lat info along X-axis
>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.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("./wrfout_d02_NoahLSM_51.nc","r")
>
>
> ; We generate plots, but what kind do we prefer?
> type = "x11"
> ; type = "pdf"
> ; type = "ps"
> ; type = "ncgm"
> wks = gsn_open_wks(type,"Seccion_zoom")
>
>
> ; Set some basic resources
> res = True
> res@MainTitle = "REAL-TIME WRF"
> res@Footer = False
> pltres = True
>
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> 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
> nd = dimsizes(mdims)
>
> xlat = wrf_user_getvar(a, "XLAT",0)
> xlon = wrf_user_getvar(a, "XLONG",0)
> ter = wrf_user_getvar(a, "HGT",0)
>
> ;---------------------------------------------------------------
>
> do it = 44,64,1 ; My 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",0) ; 3D componente U en mass points
> v = wrf_user_getvar(a,"va",0) ; 3D componente V en mass points
> tc = wrf_user_getvar(a,"tc",it) ; T in C
> ; Smooth T
> wrf_smooth_2d(tc,3)
> rh = wrf_user_getvar(a,"rh",it) ; relative humidity
> ; Smooth RH
> wrf_smooth_2d(rh,3)
> z = wrf_user_getvar(a, "z",it) ; grid point height
>
> if ( FirstTime ) then ; get height info for labels
> zmin = 0.
> zmax = 3. ; First 3km
> nz = floattoint(zmax + 1)
> end if
>
>
> ;************************************************************************************************
> lat = (/ 28.7, 28.7 /) ; user specified coordinate pairs W-E
> Fuerteventura Island; cambiar si no
> lon = (/ -14.6, -13.2 /) ; user specified coordinate pairs
> ;************************************************************************************************
> nm = getind_latlon2d (xlat,xlon, lat, lon) ; return 2d subscripts
>
> angle=0
> plane = new(4,float)
> plane = (/ nm(0,1), nm(0,0), nm(1,1), nm(1,0) /) ; approx. start x;y and
> end x;y point
> opts = True
>
> X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
> X_desc = "Longitude"
> ; Find the index where 3km is - only need to do this once
> 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-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 + 1)
>
> ;*************************************************************************************************
> ; 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@tmXBLabels = sprintf("%.1f",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("%.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
> ;*************************************************************************************************
> rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts)
> tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts)
> u_plane = wrf_user_intrp3d(u,z,"v",plane,0,opts_xy)
> v_plane = wrf_user_intrp3d(v,z,"v",plane,0,opts_xy)
>
>
> ; 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","White","White", \
> "White","Steelblue1", \
> "Steelblue2","Steelblue3"/)
> ;******************************************************************************************************
> ; Plotting options for Wind Vectors
> opts_wv = opts_xy
> opts_wv@FieldTitle = "Wind" ; overwrite Field Title
> opts_wv@NumVectors = 50 ; density of wind barbs
> opts_wv@vcGlyphStyle ="WindBarb"
> vector =
> wrf_vector(a,wks,u_plane(0:zmax_pos,:),v_plane(0:zmax_pos,:),opts_wv)
> delete(opts_wv)
> ;******************************************************************************************************
> ; 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,vector/),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)
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> FirstTimeMap = False
> end do ; END OF TIME LOOP
>
> end
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-- 
======================================================
Dennis J. Shea                  tel: 303-497-1361    |
P.O. Box 3000                   fax: 303-497-1333    |
Climate Analysis Section                             |
Climate & Global Dynamics Div.                       |
National Center for Atmospheric Research             |
Boulder, CO  80307                                   |
USA                        email: shea 'at' ucar.edu |
======================================================
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Apr 21 13:54:00 2010

This archive was generated by hypermail 2.1.8 : Fri Apr 23 2010 - 14:40:07 MDT