;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 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"


load "./geolocation_circle.ncl"

begin
 dir = "./"
 filename1 = "uwnd.2008.nc"
 f  = addfile(dir+filename1,"r" )
 u  = f->uwnd(:,{200},:,:)        ; u at 200 hpa level
 u1 = f->uwnd(:,{850},:,:)        ; u at 850 hpa level

 filename1 = "vwnd.2008.nc"
 f  = addfile(dir+filename1,"r" )
 v  = f->vwnd(:,{200},:,:)        ; v at 200 hpa level
 v1 = f->vwnd(:,{850},:,:)        ; v at 850 hpa level

; Calculating Shear

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

  vShear = (v1 - v)*1.96
  copy_VarMeta( v1, vShear )
  vShear@long_name = "Meridional Wind Shear"
  vShear@units     = "knots"

  windspeed = sqrt( uShear^2 + vShear^2 )
  copy_VarCoords( uShear, windspeed )
  windspeed@long_name = "Wind Shear Magnitude"
  windspeed@units     = uShear@units   
  printVarSummary(windspeed)
  printMinMax(windspeed,0)

; Draw cyclone centre and extend to 800km radii

;---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      = "Knots"         ; 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

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

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

; add circles
    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
         plotA@$circ_id$ = gsn_add_polyline(wks, plotA, circle_lon(nl,nr,:), circle_lat(nl,nr,:), plres)
      end do       ; nr
    end do         ; nl


  draw(plotA)
  frame(wks)
 end