;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 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 = "cdas1.20130514.pgrbl.grb2.nc"
 f = addfile(dir+filename1,"r" )
 u = f->U_GRD_L100(:,1,:,:)        ; u at 200 hpa level

 filename2 = "cdas1.20130514.pgrbl.grb2.nc"
 f = addfile(dir+filename2,"r" )
 v = f->V_GRD_L100(:,1,:,:)          ;v at 200 hpa level

 filename3 = "cdas1.20130514.pgrbl.grb2.nc"
 f = addfile(dir+filename3,"r" )
 u1 = f->U_GRD_L100(:,0,:,:)         ; u at 850 hpa level


 filename4 = "cdas1.20130514.pgrbl.grb2.nc"
 f = addfile(dir+filename4,"r" )
 v1 = f->V_GRD_L100(:,0,:,:)        ; v at 850 hpa level


; Calculating Shear

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

  vShear = (v1 - v)*1.96
  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 = "Knot"


; Draw cyclone centre and extend to 800km radii

load "./geolocation_circle.ncl"

;---One centre with multiple radii

    slat      =   15                     ; lat:storm
    slon      =  87                      ; lon:  storm
    srad      = (/200, 400, 800/)  ; station radii (km)
    srad_unit = 1                            ; km
    N         = 180                          ; # of points; more points nicer 'circle'
    opt       = False
;---End input

    nLoc = dimsizes(slat)                   ; # locations
    nRad = dimsizes(srad)                   ; # radii

    circle    = geolocation_circle(slat, slon, srad, srad_unit, N, opt) ; circle is type list
    printVarSummary(circle)

;---Extract variables from 'list'
    circle_lat = circle[0]                   ; clarity: explicitly extract list variables
    circle_lon = circle[1]

;---Print and examine returned variables
    printVarSummary(circle_lat)
    print("------------------")
    printVarSummary(circle_lon)
    print("------------------")

;************************************************
; create plots
;************************************************
  wks = gsn_open_wks("png","VWS_800_radii2")               ; 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       = "VWS within 800 radii (Mahasen 14MAY2013)"         ; change left string
  res@gsnRightString      = "Knot"         ; 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                  = 5           ; zoom map for North Indian Ocean
  res@mpMaxLatF                  =  30
  res@mpMinLonF                  = 65
  res@mpMaxLonF                  = 100
  res@mpCenterLonF = (res@mpMinLonF  + res@mpMaxLonF )*0.5 

  vecres                  = True            ; vector only resources
  vecres@gsnDraw          = False           ; don't draw
  vecres@gsnFrame         = False           ; don't advance frame
  vcPositionMode          = True
  vecres@vcGlyphStyle     = "WindBarb"     ; windbar vectors
 ;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

; add polygon
    plres  = True
    plres@gsMarkerIndex = 2           ; Plus sign ( + )
    plres@gsMarkerSizeF = 10.       
    plres@gsMarkerColor = "black"  
    ;plot@$"marker"$     = gsn_add_polymarker(wks, plotA, slon, slat, plres) ; "marker" is arbitrary

    colors = (/ "magenta", "orange", "red"/)
    thick  = (/    1     ,  1.5    ,  2.0 /) 

    do nl=0,nLoc-1
      do nr=0,nRad-1
         plres@gsLineColor      = colors(nr)
         plres@gsLineThicknessF = thick(nr)
         circ_id = "radii_"+nl+"_"+nr           ; any unique name
         ;plot@$circ_id$ = gsn_add_polyline(wks, plotA, circle_lon(nl,nr,:), circle_lat(nl,nr,:), plres)
      end do       ; nr
    end do         ; nl



  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
  d = gsn_add_polyline(wks, plotA, circle_lon(nl,nr,:), circle_lat(nl,nr,:), plres)
  d1= gsn_add_polymarker(wks, plotA, slon, slat, plres)
  draw(plotA)
  frame(wks)
 end