;****************************************************
; xy_28.ncl
;
; Concepts illustrated:
;   - Drawing a hodograph plot
;   - Reading data from an ASCII file
;   - Drawing XY plot curves with both lines and markers
;   - Putting labels on an XY plot
;
;****************************************************
;
; This script was contributed by Haoming Chen from
; LASG/IAP in Beijing, China.
;
; A hodograph plot is a suitable approach in
; describing diurnal variations in the wind vector
; at a given point.
;

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/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin
; Read data of an ASCII file
  fin  = "KENX_1832.dat"
  nlvl = numAsciiRow(fin) 
  ncol = 3
  dat  = asciiread(fin,(/nlvl,ncol/),"float")
  badv = -999.0

  alt  = dat(:,0)
   u    = dat(:,1)*1.94
   v    = dat(:,2)*1.94
 ;u    = dat(:,1)
 ;v    = dat(:,2)
 ul = 0.0

 print (alt)
 print (u)
 print (v)

 shear13u = u(ind(alt(:).eq.100))-u(ind(alt(:).eq.30))

 shear13v = v(ind(alt(:).eq.100))-v(ind(alt(:).eq.30))

 shear16u = u(ind(alt(:).eq.200))-u(ind(alt(:).eq.30))

 shear16v = v(ind(alt(:).eq.200))-v(ind(alt(:).eq.30))

 speed13 = (shear13u^2+shear13v^2)^.5
 dir13	  = wind_direction(shear13u,shear13v,0)
 speed16 = (shear16u^2+shear16v^2)^.5
 dir16	  = wind_direction(shear16u,shear16v,0)

 print(speed16)
 print(dir16)
  
  do i=0,nlvl-2
  	dis = (abs((u(i)-u(i+1)))^2+abs((v(i)-v(i+1))^2))^.5
  	ul = ul + dis
  end do

  print(ul)

;====================================================

  wks  = gsn_open_wks("png",fin+"_hodograph_KENX")

  res                      = True 

  res@gsnMaximize          = True    ; Maximize plot in window.

  res@gsnDraw              = False   ; Don't draw plot
  res@gsnFrame             = False   ; Don't advance frame

;  res@gsnPaperOrientation  = "portrait"

  res@gsnLeftString        = ""
  res@gsnRightString       = ""

  res@xyMarkLineMode       = "MarkLines"  ; markers and lines
  res@xyMarkers            = 16           ; filled dot
  res@xyMarkerSizeF        = 0.0035
  res@xyMarkerThicknessF   = 0.6

  res@tiMainString         = fin+"  1-3km shear: "+round(speed13,0)+"; 1-6km shear: "+round(speed16,0)
  res@tiYAxisString        = "V-Wind"
  res@tiXAxisString        = "U-Wind"
  res@trXMinF		   = -50
  res@trXMaxF		   =  50 
  res@trYMinF		   = -50
  res@trYMaxF		   =  50

  plot = gsn_csm_xy (wks,u,v,res)    ; Create plot, but don't draw it 


; Attach some text strings next to the markers.
  txres               = True
  txres@txFontHeightF = 0.0075
  txres@txJust        = "TopLeft"

; labels            = " " + ispan(2000,30000,2000)
  labels            = " " + alt(:) 
  labels(0:nlvl-1:2) 	    = ""   ; Don't do the last one.

  txid = gsn_add_text(wks,plot,labels,u,v,txres)

  draw(plot)    ; Drawing the plot will also cause text strings to 

  frame(wks)    ; Advance the frame.
end