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"

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_diff.txt"
cmap=asciiread(c,(/83,3/),"float")
cmap=cmap/255

data1=addfile("./UVZ_composite_sub.nc","r")
data2=addfile("../../../NCL_FAILURE_END/UV200/analysis/UVZ_composite_sub.nc","r")

U1=data1->U
U2=data2->U
U1!0="case"
U1!1="lat"
U1!2="lon"

V1=data1->V
V2=data2->V
Z1=data1->Z
Z2=data2->Z

;U1&case=ispan(0,91,1)
U1&lat=fspan(-20,20,161)
U1&lon=fspan(-20,20,161)

copy_VarMeta(U1,U2)
copy_VarMeta(U1,V1)
copy_VarMeta(U1,V2)
copy_VarMeta(U1,Z1)
copy_VarMeta(U1,Z2)

dsizes=dimsizes(U1)
;print(dsizes)
dcase=dsizes(0)
dlat=dsizes(1)
dlon=dsizes(2)

U1_avg=new((/dlat,dlon/),"float")
V1_avg=new((/dlat,dlon/),"float")
Z1_avg=new((/dlat,dlon/),"float")

U2_avg=new((/dlat,dlon/),"float")
V2_avg=new((/dlat,dlon/),"float")
Z2_avg=new((/dlat,dlon/),"float")

U1_avg=dim_avg_n_Wrap(U1,0)
V1_avg=dim_avg_n_Wrap(V1,0)
Z1_avg=dim_avg_n_Wrap(Z1,0)
U2_avg=dim_avg_n_Wrap(U2,0)
V2_avg=dim_avg_n_Wrap(V2,0)
Z2_avg=dim_avg_n_Wrap(Z2,0)

MAG1_avg=sqrt(U1_avg^2+V1_avg^2)
MAG2_avg=sqrt(U2_avg^2+V2_avg^2)

copy_VarMeta(U1(0,:,:),MAG1_avg)
copy_VarMeta(U2(0,:,:),MAG2_avg)

;U_avg_diff=U1_avg-U2_avg
;V_avg_diff=V1_avg-V2_avg
Z_avg_diff=Z2_avg-Z1_avg ;height fall during ERC
;MAG_avg_diff=sqrt(U_avg_diff^2+V_avg_diff^2)

;copy_VarMeta(U1_avg,U_avg_diff)
;copy_VarMeta(V1_avg,V_avg_diff)
copy_VarMeta(Z1_avg,Z_avg_diff)
;copy_VarMeta(U1_avg,MAG_avg_diff)

wks=gsn_open_wks("png","composite_height_fall_sub")
;gsn_define_colormap(wks,"MPL_RdBu")
;gsn_reverse_colormap(wks)
gsn_define_colormap(wks,cmap)

res=True

res@gsnAddCyclic=False
res@gsnFrame=False
res@gsnDraw=False
res@gsnMaximize=False
res@gsnStringFont=21
res@gsnLeftString="200hPa Wind (kt), Geopotential Height (gpdm,contour), Height Fall (gpdm,shading)"
res@gsnLeftStringOrthogonalPosF=0.008
res@gsnLeftStringFontHeightF=0.013
res@gsnRightString=" "
;res@gsnRightString="composite of 53 cases"
;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@trYMinF=-20
res@trYMaxF=20
res@trXMinF=-20
res@trXMaxF=20

;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(-15,15,5)
res@tmXBLabels            = res@tmXBValues
;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(-15,15,5)
res@tmYLLabels            = res@tmYLValues
;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(-10,10,41)
res@cnFillOpacityF=1
res@cnLinesOn=False

res@lbBoxSeparatorLinesOn=False
res@lbLabelStride=10
res@lbLabelFontHeightF=0.015

res@pmLabelBarHeightF=0.1

plot1=gsn_csm_contour(wks,Z_avg_diff,res)

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

zmin=min(Z2_avg)
zmax=max(Z2_avg)

res@cnLevels:=ispan(floattoint(zmin),floattoint(zmax),5)
res@cnLineLabelPlacementMode="constant"
res@cnLineDashSegLenF=0.25
res@cnLineLabelInterval=1
res@cnLineThicknessF=4.0
res@cnLineLabelFontHeightF=0.010
res@cnLineLabelFontThicknessF=2.0

;res@gsnContourNegLineDashPattern=2

plot2=gsn_csm_contour(wks,smth9_Wrap(Z2_avg,0.5,0.25,False),res)

U2_avg_kts=U2_avg/0.51444
V2_avg_kts=V2_avg/0.51444

copy_VarMeta(U2_avg,U2_avg_kts)
copy_VarMeta(V2_avg,V2_avg_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@vcRefAnnoArrowUseVecColor=False
res@vcRefAnnoArrowLineColor=1

res@vcRefAnnoString2On=False

res@vcRefAnnoOn=False

plot3=gsn_csm_vector(wks,U2_avg_kts(::10,::10),V2_avg_kts(::10,::10),res)

resc=True
resc@gsLineColor="red"
resc@gsLineThicknessF=10
resc@gsLineDashPattern=2
resc@gsLineDashSegLenF=0.3
circle_ll(wks,plot2,0,0,5,resc)

draw(plot1)
draw(plot2)
draw(plot3)

frame(wks)

end