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

begin

; for use of high resolution map database, execute
; 'setenv NCARG_RANGS /network/rit/lab/atm419lab/zhou/rangs' (on snow)
; 'setenv NCARG_RANGS /atm419lab/zhou/rangs' (on ash, reed)
; before running the script

c="./cmap_shear.txt"
cmap=asciiread(c,(/103,3/),"float")
cmap=cmap/255

center=asciiread("../../DATA_FAILURE_END/latlon_failure_end.txt",(/82,2/),"float")
;print(center(:,1))
tinfo=asciiread("../../DATA_FAILURE_END/time_failure_end.txt",(/82,2/),"integer")
yyyymmddhh=tinfo(:,1)
;print(yyyymmddhh)
name=asciiread("../../DATA_FAILURE_END/name_failure.txt",(/82/),"string")
;print(name)

dcase=dimsizes(yyyymmddhh)
U_comp=new((/dcase,161,161/),"float")
V_comp=new((/dcase,161,161/),"float")
Z_comp=new((/dcase,161,161/),"float")
avg_shear=new((/dcase/),"string")

do case=0,dcase-1,1

print("working on "+yyyymmddhh(case)+" -- "+name(case))

data=addfile("../../DATA_FAILURE_END/"+yyyymmddhh(case)+"_UV.grib","r")
;names=getfilevarnames(data)
;print(names)

U:=data->U_GDS0_ISBL
U!0="lev"
U!1="lat"
U!2="lon"

V:=data->V_GDS0_ISBL
V!0="lev"
V!1="lat"
V!2="lon"

Z:=data->Z_GDS0_ISBL
;printVarSummary(Z)
Z!0="lev"
Z!1="lat"
Z!2="lon"

;print(U&lev) ; levs are 200, 500, 850
;print(U&lat)

lat=U&lat
lon=U&lon
dlat=dimsizes(lat)
dlon=dimsizes(lon)
centlat_ind=ind_nearest_coord(center(case,0),lat,0)
centlon_ind=ind_nearest_coord(center(case,1),lon,0)
centlat=lat(centlat_ind)
centlon=lon(centlon_ind)

levels=(/0,2/) ;200hPa, 850hPa, respectively, sequence matters for calculation
dlev=dimsizes(levels)
u_buff=new((/2,dlat,dlon/),"float") ;storing vortex removed array for calculating shear
v_buff=new((/2,dlat,dlon/),"float")
z_buff=new((/2,dlat,dlon/),"float")

do z=0,dlev-1,1

UVAR=U(levels(z),:,:)
VVAR=V(levels(z),:,:)
MAG=sqrt(UVAR^2+VVAR^2)
ZVAR=Z(levels(z),:,:)/9.8/10.0

copy_VarMeta(U(levels(z),:,:),UVAR)
copy_VarMeta(U(levels(z),:,:),VVAR)
copy_VarMeta(U(levels(z),:,:),MAG)
copy_VarMeta(Z(levels(z),:,:),ZVAR)

pi=4.0*atan(1.0)
mlist=(/2,3,4,2,5,6,7,2,8,9,2/)
dm=dimsizes(mlist)

UVAR_1deg=UVAR(::4,::4)
VVAR_1deg=VVAR(::4,::4)
ZVAR_1deg=ZVAR(::4,::4)
;printVarSummary(UVAR_1deg)
;print(UVAR_1deg&lon)
lat_1deg=UVAR_1deg&lat
lon_1deg=UVAR_1deg&lon
dlat_1deg=dimsizes(lat_1deg)
dlon_1deg=dimsizes(lon_1deg)

UVAR_1deg_smthx=UVAR_1deg
VVAR_1deg_smthx=VVAR_1deg
ZVAR_1deg_smthx=ZVAR_1deg

do m=0,dm-1,1
  ;print("max u ="+max(UVAR_1deg_smthx)+"min u ="+min(UVAR_1deg_smthx))
  K=0.5*(1.0-cos(2.0*pi/mlist(m)))^(-1.0)

  UVAR_1deg_smthx_inter=UVAR_1deg_smthx
  VVAR_1deg_smthx_inter=VVAR_1deg_smthx
  ZVAR_1deg_smthx_inter=ZVAR_1deg_smthx

    do i=1,dlon_1deg-2,1

      UVAR_1deg_smthx(:,i)=UVAR_1deg_smthx_inter(:,i)+K*(UVAR_1deg_smthx_inter(:,i-1)+UVAR_1deg_smthx_inter(:,i+1)-2.0*UVAR_1deg_smthx_inter(:,i))
      VVAR_1deg_smthx(:,i)=VVAR_1deg_smthx_inter(:,i)+K*(VVAR_1deg_smthx_inter(:,i-1)+VVAR_1deg_smthx_inter(:,i+1)-2.0*VVAR_1deg_smthx_inter(:,i))
      ZVAR_1deg_smthx(:,i)=ZVAR_1deg_smthx_inter(:,i)+K*(ZVAR_1deg_smthx_inter(:,i-1)+ZVAR_1deg_smthx_inter(:,i+1)-2.0*ZVAR_1deg_smthx_inter(:,i))

    end do

end do

UVAR_1deg_smthy=UVAR_1deg_smthx
VVAR_1deg_smthy=VVAR_1deg_smthx
ZVAR_1deg_smthy=ZVAR_1deg_smthx

do m=0,dm-1,1
  ;print("max u ="+max(UVAR_1deg_smthy)+"min u ="+min(UVAR_1deg_smthy))
  K=0.5*(1.0-cos(2.0*pi/mlist(m)))^(-1)

  UVAR_1deg_smthy_inter=UVAR_1deg_smthy
  VVAR_1deg_smthy_inter=VVAR_1deg_smthy
  ZVAR_1deg_smthy_inter=ZVAR_1deg_smthy

    do j=1,dlat_1deg-2,1

      UVAR_1deg_smthy(j,:)=UVAR_1deg_smthy_inter(j,:)+K*(UVAR_1deg_smthy_inter(j-1,:)+UVAR_1deg_smthy_inter(j+1,:)-2.0*UVAR_1deg_smthy_inter(j,:))
      VVAR_1deg_smthy(j,:)=VVAR_1deg_smthy_inter(j,:)+K*(VVAR_1deg_smthy_inter(j-1,:)+VVAR_1deg_smthy_inter(j+1,:)-2.0*VVAR_1deg_smthy_inter(j,:))
      ZVAR_1deg_smthy(j,:)=ZVAR_1deg_smthy_inter(j,:)+K*(ZVAR_1deg_smthy_inter(j-1,:)+ZVAR_1deg_smthy_inter(j+1,:)-2.0*ZVAR_1deg_smthy_inter(j,:))

    end do

end do

UVAR_1deg_smth=UVAR_1deg_smthy
VVAR_1deg_smth=VVAR_1deg_smthy
ZVAR_1deg_smth=ZVAR_1deg_smthy

;print("max u ="+max(UVAR_1deg_smth)+"min u ="+min(UVAR_1deg_smth))
;print("max v ="+max(VVAR_1deg_smth)+"min v ="+min(VVAR_1deg_smth))
;print("max z ="+max(ZVAR_1deg_smth)+"min z ="+min(ZVAR_1deg_smth))

UVAR_smth=linint2(lon_1deg,lat_1deg(::-1),UVAR_1deg_smth,True,lon,lat(::-1),0) ;lat should be increasing
VVAR_smth=linint2(lon_1deg,lat_1deg(::-1),VVAR_1deg_smth,True,lon,lat(::-1),0) ;as long as two sets of lat are aligned,
ZVAR_smth=linint2(lon_1deg,lat_1deg(::-1),ZVAR_1deg_smth,True,lon,lat(::-1),0) ;original data sequence doesn't matter

;print(num(ismissing(UVAR_smth)))
;print(num(ismissing(VVAR_smth)))
;print(num(ismissing(ZVAR_smth)))
;printVarSummary(UVAR_smth)

copy_VarMeta(UVAR,UVAR_smth) ;otherwise, lat will be flipped
copy_VarMeta(VVAR,VVAR_smth)
copy_VarMeta(ZVAR,ZVAR_smth)

UVAR_dev=UVAR-UVAR_smth
VVAR_dev=VVAR-VVAR_smth
ZVAR_dev=ZVAR-ZVAR_smth

copy_VarMeta(UVAR,UVAR_dev) 
copy_VarMeta(VVAR,VVAR_dev)
copy_VarMeta(ZVAR,ZVAR_dev)

;print("max Udev = "+max(UVAR_dev)+" min Udev = "+min(UVAR_dev))

;rn=31 ;radius of interpolation domain, set bigger than needed (7.5 deg instead of 5 deg)
rn=61 ;15 deg
tn=72 
dr=0.25 ; deg lat/lon
dtheta=5 ;deg
theta=ispan(0,360-1,dtheta)

lon_cyl2d=new((/rn,tn/),float)
lat_cyl2d=new((/rn,tn/),float)

do r=0,rn-1,1

  lon_cyl2d(r,:)=lon(centlon_ind)+r*dr*cos(theta*pi/180)
  lat_cyl2d(r,:)=lat(centlat_ind)+r*dr*sin(theta*pi/180)

end do

lon_cyl1d=ndtooned(lon_cyl2d)
lat_cyl1d=ndtooned(lat_cyl2d)

;print(lon_cyl1d)

UVAR_dev_cyl1d=linint2_points(lon,lat(::-1),UVAR_dev(::-1,:),True,lon_cyl1d,lat_cyl1d,0) ;need to reverse original data 
VVAR_dev_cyl1d=linint2_points(lon,lat(::-1),VVAR_dev(::-1,:),True,lon_cyl1d,lat_cyl1d,0) ;corresponding to lat array
ZVAR_dev_cyl1d=linint2_points(lon,lat(::-1),ZVAR_dev(::-1,:),True,lon_cyl1d,lat_cyl1d,0)

UVAR_dev_cyl2d=reshape(UVAR_dev_cyl1d,(/rn,tn/))
VVAR_dev_cyl2d=reshape(VVAR_dev_cyl1d,(/rn,tn/))
ZVAR_dev_cyl2d=reshape(ZVAR_dev_cyl1d,(/rn,tn/))

;print("max Udev = "+max(UVAR_dev)+" min Udev = "+min(UVAR_dev))
;print("max Udev.interp = "+max(UVAR_dev_cyl1d)+" min Udev.interp = "+min(UVAR_dev_cyl1d))

;r0=21 ;radius of filtering domain (5 deg)
r0=41
l=r0/5.0 ;parameter controling filtering shape

UVAR_vortex_cyl2d=UVAR_dev_cyl2d
VVAR_vortex_cyl2d=VVAR_dev_cyl2d
ZVAR_vortex_cyl2d=ZVAR_dev_cyl2d

;print("max Uvortex = "+max(UVAR_vortex_cyl2d)+" min Uvortex = "+min(UVAR_vortex_cyl2d))

do r=0,rn-1,1
  
  if (r.le.r0) then

    E=(exp(-((r0-r)^2/l^2))-exp(-((r0)^2/l^2)))/(1.0-exp(-((r0)^2/l^2)))
    ;print("E = "+E)
    UVAR_vortex_cyl2d(r,:)=UVAR_dev_cyl2d(r,:)-UVAR_dev_cyl2d(r0,:)*E-avg(UVAR_dev_cyl2d(r0,:))*(1-E)
    VVAR_vortex_cyl2d(r,:)=VVAR_dev_cyl2d(r,:)-VVAR_dev_cyl2d(r0,:)*E-avg(VVAR_dev_cyl2d(r0,:))*(1-E)
    ZVAR_vortex_cyl2d(r,:)=ZVAR_dev_cyl2d(r,:)-ZVAR_dev_cyl2d(r0,:)*E-avg(ZVAR_dev_cyl2d(r0,:))*(1-E)

  else

    UVAR_vortex_cyl2d(r,:)=0
    VVAR_vortex_cyl2d(r,:)=0
    ZVAR_vortex_cyl2d(r,:)=0

  end if

end do

;print("max Uvortex = "+max(UVAR_vortex_cyl2d)+" min Uvortex = "+min(UVAR_vortex_cyl2d))
;small discrepancy should be a result of TC location difference

;lon_sub=lon(centlon_ind-20:centlon_ind+20)
;lat_sub=lat(centlat_ind+20:centlat_ind-20)

lon_sub=lon(centlon_ind-40:centlon_ind+40)
lat_sub=lat(centlat_ind+40:centlat_ind-40)

UVAR_vortex_sub:=transpose(dsgrid2(lon_cyl1d,lat_cyl1d,ndtooned(UVAR_vortex_cyl2d),lon_sub,lat_sub)) ;dsgrid returns dx*dy array
VVAR_vortex_sub:=transpose(dsgrid2(lon_cyl1d,lat_cyl1d,ndtooned(VVAR_vortex_cyl2d),lon_sub,lat_sub)) ;instead of dy*dx !!!!!!!!!
ZVAR_vortex_sub:=transpose(dsgrid2(lon_cyl1d,lat_cyl1d,ndtooned(ZVAR_vortex_cyl2d),lon_sub,lat_sub))

;print("max Uvortex_sub = "+max(UVAR_vortex_sub)+" min Uvortex_sub = "+min(UVAR_vortex_sub))
;printVarSummary(UVAR_vortex_sub)

;copy_VarMeta(UVAR(centlat_ind+20:centlat_ind-20,centlon_ind-20:centlon_ind+20),UVAR_vortex_sub)
;copy_VarMeta(VVAR(centlat_ind+20:centlat_ind-20,centlon_ind-20:centlon_ind+20),VVAR_vortex_sub)
;copy_VarMeta(ZVAR(centlat_ind+20:centlat_ind-20,centlon_ind-20:centlon_ind+20),ZVAR_vortex_sub)

copy_VarMeta(UVAR(centlat_ind+40:centlat_ind-40,centlon_ind-40:centlon_ind+40),UVAR_vortex_sub)
copy_VarMeta(VVAR(centlat_ind+40:centlat_ind-40,centlon_ind-40:centlon_ind+40),VVAR_vortex_sub)
copy_VarMeta(ZVAR(centlat_ind+40:centlat_ind-40,centlon_ind-40:centlon_ind+40),ZVAR_vortex_sub)

;printVarSummary(UVAR_vortex_sub)

UVAR_vortex_rm=UVAR
VVAR_vortex_rm=VVAR
ZVAR_vortex_rm=ZVAR

;print("max Uvortex_rm = "+max(UVAR_vortex_rm)+" min Uvortex_rm = "+min(UVAR_vortex_rm))
;printVarSummary(UVAR_vortex_rm)
;printVarSummary(UVAR_vortex_rm(centlat_ind+20:centlat_ind-20,centlon_ind-20:centlon_ind+20))
;printVarSummary(UVAR_vortex_sub)

;UVAR_vortex_rm(centlat_ind+20:centlat_ind-20,centlon_ind-20:centlon_ind+20)= \
;  UVAR_vortex_rm(centlat_ind+20:centlat_ind-20,centlon_ind-20:centlon_ind+20)-UVAR_vortex_sub
;VVAR_vortex_rm(centlat_ind+20:centlat_ind-20,centlon_ind-20:centlon_ind+20)= \
;  VVAR_vortex_rm(centlat_ind+20:centlat_ind-20,centlon_ind-20:centlon_ind+20)-VVAR_vortex_sub
;ZVAR_vortex_rm(centlat_ind+20:centlat_ind-20,centlon_ind-20:centlon_ind+20)= \
;  ZVAR_vortex_rm(centlat_ind+20:centlat_ind-20,centlon_ind-20:centlon_ind+20)-ZVAR_vortex_sub

UVAR_vortex_rm(centlat_ind+40:centlat_ind-40,centlon_ind-40:centlon_ind+40)= \
  UVAR_vortex_rm(centlat_ind+40:centlat_ind-40,centlon_ind-40:centlon_ind+40)-UVAR_vortex_sub
VVAR_vortex_rm(centlat_ind+40:centlat_ind-40,centlon_ind-40:centlon_ind+40)= \
  VVAR_vortex_rm(centlat_ind+40:centlat_ind-40,centlon_ind-40:centlon_ind+40)-VVAR_vortex_sub
ZVAR_vortex_rm(centlat_ind+40:centlat_ind-40,centlon_ind-40:centlon_ind+40)= \
  ZVAR_vortex_rm(centlat_ind+40:centlat_ind-40,centlon_ind-40:centlon_ind+40)-ZVAR_vortex_sub

MAG_vortex_rm=sqrt(UVAR_vortex_rm^2+VVAR_vortex_rm^2)
copy_VarMeta(UVAR_vortex_rm,MAG_vortex_rm)

;printVarSummary(UVAR_vortex_rm)

;print(" "+centlat+" "+lat(centlat_ind+28)+" "+lat(centlat_ind-28))
;print(" "+centlon+" "+lon(centlon_ind+28)+" "+lon(centlon_ind-28))

u_buff(z,:,:)=UVAR_vortex_rm
v_buff(z,:,:)=VVAR_vortex_rm
z_buff(z,:,:)=ZVAR_vortex_rm

end do

;printVarSummary(u_buff)

shear_u=u_buff(0,:,:)-u_buff(1,:,:)
shear_v=v_buff(0,:,:)-v_buff(1,:,:)
shear_mag=sqrt(shear_u^2+shear_v^2)

copy_VarMeta(UVAR,shear_u)
copy_VarMeta(UVAR,shear_v)
copy_VarMeta(UVAR,shear_mag)
;printVarSummary(shear_u)

shear_u_cyl1d=linint2_points(lon,lat(::-1),shear_u(::-1,:),True,lon_cyl1d,lat_cyl1d,0) ;for calculating shear within 5 deg of TC center
shear_v_cyl1d=linint2_points(lon,lat(::-1),shear_v(::-1,:),True,lon_cyl1d,lat_cyl1d,0)
shear_mag_cyl1d=linint2_points(lon,lat(::-1),shear_mag(::-1,:),True,lon_cyl1d,lat_cyl1d,0)

shear_u_cyl2d=reshape(shear_u_cyl1d,(/rn,tn/))
shear_v_cyl2d=reshape(shear_v_cyl1d,(/rn,tn/))
shear_mag_cyl2d=reshape(shear_mag_cyl1d,(/rn,tn/))

;printVarSummary(shear_u_cyl2d)

avg_shear_u=avg(shear_u_cyl2d(0:20,:)) ;5 deg
avg_shear_v=avg(shear_v_cyl2d(0:20,:))
avg_shear_mag=avg(shear_mag_cyl2d(0:20,:))

avg_shear(case)=sprintf("%9.2f",avg_shear_u)+sprintf("%9.2f",avg_shear_v)+sprintf("%9.2f",avg_shear_mag)

U_comp(case,:,:)=shear_u(centlat_ind+80:centlat_ind-80,centlon_ind-80:centlon_ind+80)
V_comp(case,:,:)=shear_v(centlat_ind+80:centlat_ind-80,centlon_ind-80:centlon_ind+80)

wks=gsn_open_wks("png","failure_end_case_"+case)
gsn_define_colormap(wks,cmap)

res=True

res@gsnAddCyclic=False
res@gsnFrame=False
res@gsnDraw=False
res@gsnMaximize=False
res@gsnStringFont=21
res@gsnLeftString="200-850hPa shear (kt)"
res@gsnLeftStringOrthogonalPosF=0.01
res@gsnLeftStringFontHeightF=0.015
res@gsnRightString=name(case)+" @ "+sprinti("%10.0i",yyyymmddhh(case))
;res@gsnRightString="+"+sprinti("%0.2i",counter)+"h"
;res@gsnRightString="+"+sprintf("%4.1f",counter*0.5)+"h"
res@gsnRightStringOrthogonalPosF=0.01
res@gsnRightStringFontHeightF=0.015
res@gsnScalarContour=True

res@vpXF=0.15
res@vpYF=0.95
res@vpHeightF=0.7
res@vpWidthF=0.7

;res@stLineColor=0
;res@stArrowLengthF=0.004
;res@stMinArrowSpacingF=0.004
;res@stMinDistanceF=0.1
;res@stLineThicknessF=2.0
;res@stStepSizeF=0.01

res@mpDataResolution="FinestResolution"
res@mpDataBaseVersion="HighRes"
res@mpLimitMode="LatLon"
;res@mpCenterLatF=centlat
;res@mpCenterLonF=centlon
res@mpMinLatF=centlat-20
res@mpMaxLatF=centlat+20
res@mpMinLonF=centlon-20
res@mpMaxLonF=centlon+20
res@mpMonoFillColor=False
res@mpFillColors=(/-1,-1,-1/)
res@mpOutlineOn=True

res@tmXTOn                = False
res@tmYROn                = False

res@tmXBMode              = "Explicit"
res@tmXBValues            = ispan(100,180,10)
res@tmXBLabels            = (/"100~S~o~N~E","110~S~o~N~","120~S~o~N~E",\
                              "130~S~o~N~","140~S~o~N~E","150~S~o~N~",\
                              "160~S~o~N~E","170~S~o~N~","180~S~o~N~E"/)
res@tmXBMajorLengthF      = 0.01
res@tmXBMajorOutwardLengthF=0.0
res@tmXBMajorThicknessF   = 3.0
;res@tmXBMinorOn           = False
res@tmXBMinorPerMajor     = 4
res@tmXBMinorLengthF      = 0.004
res@tmXBMinorOutwardLengthF=0.0
res@tmXBMinorThicknessF   = 2.0
res@tmXBLabelFont         = 21
res@tmXBLabelFontHeightF  = 0.015
res@tmXBLabelDeltaF       = -0.6

res@tmYLMode              = "Explicit"
res@tmYLValues            = ispan(0,50,10)
res@tmYLLabels            = (/"0~S~o~N~","10~S~o~N~N","20~S~o~N~",\
                              "30~S~o~N~N","40~S~o~N~","50~S~o~N~N"/)
res@tmYLMajorLengthF      = 0.01
res@tmYLMajorOutwardLengthF=0.0
res@tmYLMajorThicknessF   = 3.0
;res@tmYLMinorOn           = False
res@tmYLMinorPerMajor     = 4
res@tmYLMinorLengthF      = 0.004
res@tmYLMinorOutwardLengthF=0.0
res@tmYLMinorThicknessF   = 2.0
res@tmYLLabelFont         = 21
res@tmYLLabelFontHeightF  = 0.015
res@tmYLLabelDeltaF       = -0.3

res@cnFillOn=True
res@cnLevelSelectionMode="ExplicitLevels"
;res@cnLevels:=fspan(0,80,81)
res@cnLevels:=fspan(0,100,101)
res@cnFillOpacityF=1
res@cnLinesOn=False

res@lbBoxSeparatorLinesOn=False
res@lbLabelStride=10

res@pmLabelBarHeightF=0.1

shear_mag_kts=shear_mag/0.51444

copy_VarMeta(UVAR,shear_mag_kts)

plot1=gsn_csm_contour_map(wks,shear_mag_kts,res)

;res@cnFillOn=False
;res@cnLinesOn=True
;res@cnInfoLabelOn=False

;res@cnLevels:=ispan(floattoint(min(ZVAR_vortex_rm)),floattoint(max(ZVAR_vortex_rm)),5)
;res@cnLineLabelPlacementMode="constant"
;res@cnLineDashSegLenF=0.4
;res@cnLineColor="black"
;res@cnLineLabelFontColor="black"
;res@cnLineLabelFontHeightF=0.010
;res@cnLineLabelFontThicknessF=2.0
;res@cnLineLabelInterval=1
;res@cnLineThicknessF=5.0

;plot2=gsn_csm_contour_map(wks,smth9_Wrap(ZVAR_vortex_rm,0.5,0.25,False),res)

shear_u_kts=shear_u/0.51444
shear_v_kts=shear_v/0.51444

copy_VarMeta(UVAR,shear_u_kts)
copy_VarMeta(VVAR,shear_v_kts)

res@vcRefMagnitudeF=65
res@vcRefLengthF=0.025
;res@vcGlyphStyle="CurlyVector"
res@vcGlyphStyle="WindBarb"
;res@vcLineArrowThicknessF=3.0
res@vcWindBarbLineThicknessF=2.0
res@vcWindBarbTickLengthF=0.4
res@vcRefAnnoOrthogonalPosF=-0.079
res@vcRefAnnoParallelPosF=1.0
;res@vcMinDistanceF=0.03
res@vcMonoLineArrowColor=True
;res@vcLineArrowColor=1
res@vcWindBarbColor="white"
res@vcRefAnnoArrowUseVecColor=False
res@vcRefAnnoArrowLineColor="black"

res@vcWindBarbTickSpacingF=0.125
res@vcRefAnnoOn=True
res@vcMonoWindBarbColor=True
res@vcGlyphOpacityF=1.0

plot3=gsn_csm_vector_map(wks,shear_u_kts(::10,::10),shear_v_kts(::10,::10),res)

shear_vec_u=UVAR
shear_vec_v=VVAR

shear_vec_u=0
shear_vec_v=0

shear_vec_u(centlat_ind,centlon_ind)=avg_shear_u/0.51444
shear_vec_v(centlat_ind,centlon_ind)=avg_shear_v/0.51444

shear_vec_u@_FillValue=0
shear_vec_v@_FillValue=0

res@vcRefLengthF=0.090
res@vcWindBarbLineThicknessF=15.0
res@vcWindBarbColor="cyan"
res@vcWindBarbTickSpacingF=0.2

res@vcRefAnnoOn=False
res@vcScalarMissingValColor=-1
res@vcGlyphOpacityF=1.0

plot4=gsn_csm_vector_map(wks,shear_vec_u,shear_vec_v,res)

resc=True
resc@gsLineColor="cyan" ;RGB for "cyan" is (0 255 255)
resc@gsLineThicknessF=10
resc@gsLineDashPattern=2
resc@gsLineDashSegLenF=0.3
circle_ll(wks,plot1,centlon,centlat,5,resc)

draw(plot1)
;draw(plot2)
draw(plot4)
draw(plot3)

;angle=asin(avg_shear_v/avg_shear_u)*180/pi+180
;wmsetp("ARD",angle)
;size=0.05*avg_shear_mag/5.0
;wmsetp("ARS",size)
;length=0.5*avg_shear_mag/5.0
;wmsetp("ARL",length)
;new_color=NhlNewColor(wks,0.0,1.0,1.0)
;;wmsetp("ASC",103)
;wmsetp("AWC",103)

;;wmlabs(wks,centlon,centlat,"Arrow")

;wmsetp("VCH",0.02)
;wmsetp("VCC",103)
;wmsetp("VRS",0.3)

;wmvect(wks,centlon,centlat,avg_shear_u,avg_shear_v)

frame(wks)

end do

rm=systemfunc("rm -rf UVZ_composite_all.nc")

fout=addfile("./shear_composite_all.nc","c")
fout->U=U_comp
fout->V=V_comp
;fout->Z=Z_comp

rm=systemfunc("rm -rf avg_shear.txt")
asciiwrite("./avg_shear.txt",avg_shear)

end