;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Description: Calculate vertical wind shear and plot
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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/contrib/time_axis_labels.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"


begin
 dir = "./"
 filename1 = "200.nc"
 f = addfile(dir+filename1,"r" )
 u = f->U_GRD_L100_Avg        ; u at 200 hpa level

 filename2 = "200.nc"
 f = addfile(dir+filename2,"r" )
 v = f->V_GRD_L100_Avg         ;v at 200 hpa level

 filename3 = "850.nc"
 f = addfile(dir+filename3,"r" )
 u1 = f->U_GRD_L100_Avg        ; u at 850 hpa level

 filename4 = "850.nc"
 f = addfile(dir+filename4,"r" )
 v1 = f->V_GRD_L100_Avg        ; v at 850 hpa level

; Calculating Shear

  uShear = u1 - u
  copy_VarMeta( u1, uShear )
  uShear@long_name = "Zonal Wind Shear"
  uShear@units = "m/s"

  vShear = v1 - v
  copy_VarMeta( v1, vShear )
  vShear@long_name = "Meridional Wind Shear"
  vShear@units = "m/s"

  windspeed = sqrt( uShear^2 + vShear^2 )
 copy_VarMeta( uShear, windspeed )
 windspeed@long_name = "Wind Shear Magnitude"
 windspeed@units = "m/s"
;************************************************
; create plots
;************************************************
  wks = gsn_open_wks("png","Vertical_Wind_shear")               ; send graphics to PNG file
                         
;************************************************
; create plots
;************************************************
  res                     = True
  res@gsnDraw             = False           ; don't draw
  res@gsnFrame            = False           ; don't advance frame
  res@cnFillOn            = True            ; turn on color
  res@cnFillPalette       = "gui_default"   ; set color map
  res@cnLinesOn           = False           ; no contour lines
  res@gsnLeftString       = "Wind Speed"         ; change left string
  res@gsnRightString      = u@units         ; assign right string
  res@mpFillOn            = False           ; no map fill
  ;res@vcMinDistanceF = 0.02   ; thin the vectors
  ;res@vcLineArrowThicknessF = 3.0     ; 3x as thick (1.0 is the default)

  res@mpMinLatF                  = -34           ; zoom map for North Indian Ocean
  res@mpMaxLatF                  =  34
  res@mpMinLonF                  = 0
  res@mpMaxLonF                  = 130

  vecres                  = True            ; vector only resources
  vecres@gsnDraw          = False           ; don't draw
  vecres@gsnFrame         = False           ; don't advance frame
  vecres@vcGlyphStyle     = "LineArrow"     ; Line vectors
  vecres@vcRefLengthF     = 0.045           ; define length of vec ref
  vecres@gsnRightString   = " "             ; turn off right string
  vecres@gsnLeftString    = " "             ; turn off left string
  vecres@tiXAxisString    = " "             ; turn off axis label 
  vecres@vcRefAnnoOrthogonalPosF = -.535    ; move ref vector into plot


  plotA   = gsn_csm_contour_map_ce(wks,windspeed(0,:,:),res)
  plotB   = gsn_csm_vector(wks,uShear(0,:,:),vShear(0,:,:),vecres)
  overlay(plotA,plotB)                      ; result will be plotA
  draw(plotA)
  frame(wks)
 end