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/csm/shea_util.ncl"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;	Developed by Matthew Sienkiewicz
;	Stony Brook University
;	11 June 2015
;
;	This NCL script reads in an ASCII file containing four columns of aircraft data:
;		latitude
;		longitude
;		u
;		v
;
;	The aircraft flight track is plotted in blue along with 30-second averaged wind barbs in
;	knots (plotted every 60 seconds).
;
;	This script requires a custom colormap as well as the RANGS dataset, so be sure to set:
;		export NCARG_RANGS=/path/to/RANGS/
;		export NCARG_COLORMAPS=/path/to/colormaps/
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

begin

	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;IMPORT FLIGHT FILE;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	flightStuff				= asciiread("nybj_01_01.txt",-1,"float")
	m					= dimsizes(flightStuff)
	dimLength				= m / 4;
	delete(m)
	flightPath				= onedtond(flightStuff, (/dimLength,4/) )
	delete(dimLength)
	delete(flightStuff)
	lat					= flightPath(:,0)
	lon					= flightPath(:,1)
	u					= flightPath(:,2)
	v					= flightPath(:,3)
	delete(flightPath)
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;SET-UP OUTPUT AND COLORMAP;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	wks_type				= "png"
	wks_type@wkWidth			= 2500
	wks_type@wkHeight			= 2500
	wks					= gsn_open_wks(wks_type,"20140723_02_flightpath")
	gsn_define_colormap(wks,"flight")
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;SET-UP MAP RESOURCES;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	mpres					= True
	mpres@gsnFrame				= False
	mpres@gsnMaximize			= True
	mpres@gsnPaperOrientation		= "landscape"
	mpres@mpGeophysicalLineColor		= "black"
	mpres@mpPerimOn				= True
	;mpres@mpGridAndLimbOn			= True
	;mpres@mpGridLineDashPattern		= 2
	;mpres@mpGridLatSpacingF		= 0.5
	;mpres@mpGridLonSpacingF		= 0.5
	;mpres@mpGridLineColor			= "black"
	mpres@mpOutlineOn			= True
	mpres@mpUSStateLineColor		= "black"
	mpres@mpDataBaseVersion			= "HighRes"
	mpres@mpFillOn				= True
	mpres@mpFillDrawOrder			= "PreDraw"
	mpres@pmTickMarkDisplayMode		= "NoCreate"
        mpres@mpOceanFillColor			= "paleturquoise"
	mpres@mpInlandWaterFillColor		= "paleturquoise"
	mpres@mpLandFillColor			= "palegreen"
	mpres@mpPerimLineThicknessF		= 20.0
	mpres@mpGeophysicalLineThicknessF 	= 6.0
	mpres@mpUSStateLineThicknessF		= 6.0
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;DEFINE CORNERS OF DOMAIN;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	mpres@mpLimitMode			= "Corners"
	mpres@mpLeftCornerLatF			= 39.69613
	mpres@mpLeftCornerLonF			= -74.48114
	mpres@mpRightCornerLatF			= 40.90421
	mpres@mpRightCornerLonF			= -72.48978
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;DEFINE LCC PROJECTION;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	mpres@mpProjection			= "LambertConformal"
	mpres@mpLambertParallel1F		= 40.79045
	mpres@mpLambertParallel2F		= 40.79045
	mpres@mpLambertMeridianF		= -72.29144
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;PLOT MAP;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	map					= gsn_csm_map(wks,mpres)
	draw(map)
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;ADD STATE BOUNDARIES TO MAP;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	mpres@mpDataBaseVersion			= "MediumRes"
	mpres@mpOutlineBoundarySets		= "USStates"
        mpres@mpOceanFillColor			= "transparent"
	mpres@mpInlandWaterFillColor		= "transparent"
	mpres@mpLandFillColor			= "transparent"
	mpres@mpGeophysicalLineThicknessF	= 0.0
	mpres@mpOutlineDrawOrder		= "PostDraw"
	map2					= gsn_csm_map(wks,mpres)
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;DRAW FLIGHT PATH ON MAP;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	pres					= True
	pres@gsLineThicknessF			= 6.0
	pres@gsLineColor			= 4
	pres@gsLineOpacityF 			= 0.8
	gsn_polyline(wks,map,lon,lat,pres)
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;DRAW WIND BARBS;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	wmsetp("wbs",0.03)
	wmsetp("wbd",0.15)
	wmsetp("wbt",0.5)
	wmsetp("blw",8.0)
	length = dimsizes(lat)
	do i=14,length-15,60
		wsp = round((sqrt(((avg(u(i-14:i+14))*1.9438444924574)^2)+((avg(v(i-14:i+14))*1.9438444924574)^2))/5),3)*5
		if (wsp.ge.60) then
			wmsetp("col",14)
		else if (wsp.eq.55) then
			wmsetp("col",13)
		else if (wsp.eq.50) then
			wmsetp("col",12)
		else if (wsp.eq.45) then
			wmsetp("col",11)
		else if (wsp.eq.40) then
			wmsetp("col",10)
		else if (wsp.eq.35) then
			wmsetp("col",9)
		else if (wsp.eq.30) then
			wmsetp("col",8)
		else if (wsp.eq.25) then
			wmsetp("col",7)
		else if (wsp.eq.20) then
			wmsetp("col",6)
		else if (wsp.eq.15) then
			wmsetp("col",5)
		else if (wsp.eq.10) then
			wmsetp("col",4)
		else if (wsp.eq.5) then
			wmsetp("col",3)
		else if (wsp.eq.0) then
			wmsetp("col",2)
		end if
		end if
		end if
		end if
		end if
		end if
		end if
		end if
		end if
		end if
		end if
		end if
		end if
		wmbarbmap(wks,lat(i),lon(i),avg(u(i-14:i+14))*1.9438444924574,avg(v(i-14:i+14))*1.9438444924574)
	end do
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;ADD TEXT;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	tres					= True
	tres@txJust				= "TopLeft"
	tres@txFontHeightF			= 0.020
	tres@txBackgroundFillColor		= "transparent"
	tres@txFontColor			= 1
	tres@txFont				= 25
	gsn_text_ndc(wks,"ACFT: 23-Jul-2014",0.030,0.86,tres)
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;ADD COLORBAR;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	cres					= True
	cres@lbPerimOn				= False
	cres@lbOrientation			= "Horizontal"
	cres@vpWidthF				= 0.4
	cres@vpHeightF				= 0.09
	cres@lbLabelFontHeightF			= 0.009
	cres@lbLabelAlignment			= "BoxCenters"
	cres@lbMonoFillPattern			= True
	cres@lbFillColors			= (/2,3,4,5,6,7,8,9,10,11,12,13,14/)
	cres@lbTitleString			= "Wind Speed [KTS]"
	cres@lbTitleFontHeightF			= 0.012
	cres@lbTitleFont			= 25
	gsn_labelbar_ndc(wks,13,(/"0","5","10","15","20","25","30","35","40","45","50","55","60+"/),0.51,0.39,cres)
	;drawNDCGrid(wks)
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;ADVANCE FRAME;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	frame(wks)
	delete(map)
	delete(wks)
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

end