NCL cross section with wind barbs

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

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
Received on Wed Apr 21 08:39:57 2010