; load libraries, twice for good measure

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/wrf/WRFUserARW.ncl"
;load "/home/torn/ncl/wrf_map_info.ncl"
;load "/home/torn/ncl/array_stat.ncl"
;load "/home/torn/ncl/image_conv.ncl"

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/wrf/WRFUserARW.ncl"
;load "/home/torn/ncl/wrf_map_info.ncl"
;load "/home/torn/ncl/array_stat.ncl"
;load "/home/torn/ncl/image_conv.ncl"
;load "wrf_dart_util.ncl"

begin

; ------ Load Data

rows = numAsciiRow("usp01.dat")
columns = numAsciiCol("usp01.dat")

data01=asciiread("usp01.dat",(/rows,columns/),"float")
templat = data01(:,2)
templon = data01(:,3)
tempdata= data01(:,4)

p = onedtond(tempdata,(/7,9/))
lat = onedtond(templat,(/7,9/))
lon = onedtond(templon,(/7,9/))
p(0,0)=p(0,0)+.01

; ------ Assign Coordinates
uniquelat = lat(:,1)
uniquelon = lon(1,:)

uniquelat@units="degrees_north"
uniquelon@units="degrees_east"

latmin=min(lat)
latmax=max(lat)
lonmin=min(lon)
lonmax=max(lon)
fillvarlevel = p
fillvarlevel!0="lat"
fillvarlevel!1="lon"
fillvarlevel&lon = uniquelon
fillvarlevel&lat = uniquelat

;---Read in site info

sitefile  = "threemonthly.dat"
sitedata = asciiread(sitefile,-1,"string")
delim = ","
sitelats  = tofloat(str_get_field(sitedata,1,delim))
sitelons  = tofloat(str_get_field(sitedata,2,delim))
siteanom  = tofloat(str_get_field(sitedata,3,delim))

wks_type = "png"
wks_type@wkWidth = 2000
wks_type@wkHeight = 1600 
wks = gsn_open_wks(wks_type,"random")

res = True
res@gsnFrame      = False
res@mpProjection = "CylindricalEquidistant"
;res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
res@gsnMaximize = True
res@lbBoxMinorExtentF = 0.25
res@cnFillOn = False   ; this controls whether or not a label bar exists
res@cnLinesOn = False   ;no contours on top of fill
res@cnFillMode = "RasterFill"
res@cnRasterSmoothingOn = True
res@mpMinLonF = lonmin ;-- min longitude
res@mpMaxLonF = lonmax ;-- max longitude
res@mpMinLatF = latmin ;-- min latitude
res@mpMaxLatF = latmax ;-- max latitude
res@tfDoNDCOverlay = True ;this is a resource for overlaying multiple plots
res@cnInfoLabelOn   = False  ; turn off contour label

res@mpDataBaseVersion           = "Ncarg4_1"
res@mpDataSetName               = "Earth..4"   ; For counties
res@mpFillAreaSpecifiers        = (/"Land", "United States:States","Water"/) ; for states only
res@mpFillAreaSpecifiers        = (/"Land", "Land:North America:United States:Conterminous US:New York","Water"/)
res@mpSpecifiedFillColors       = (/"transparent","transparent",       "transparent"/)
res@mpInlandWaterFillColor      = 238
res@mpLandFillColor             = -1     ; transparent
res@mpOutlineOn                 = True
res@mpMaskOutlineSpecifiers =   (/"Land"/)
res@mpOutlineSpecifiers =   (/ "Land:North America:United States:Conterminous US:New York"/)

res@mpGeophysicalLineThicknessF = 4
res@mpUSStateLineThicknessF = 4
res@mpNationalLineThicknessF = 4
res@mpGridLonSpacingF = 2
res@mpDataBaseVersion = "MediumRes" ;-- better map resolution
res@gsnAddCyclic = False

res@cnGridBoundFillColor = "Transparent"


  res@tmXBBorderOn      = False
  res@tmXBMinorPerMajor = 1
  res@tmXBOn            = False
  res@tmXTBorderOn      = False
  res@tmXTOn            = False
  res@tmYLBorderOn      = False
  res@tmYLMinorPerMajor = 1
  res@tmYLOn            = False
  res@tmYRBorderOn      = False
  res@tmYROn            = False

res@lbLabelBarOn        = False  
res@gsnDraw             = False

size = dimsizes(sitelons)
res@mpOutlineSpecifiers =   (/ "Land:North America:United States:Conterminous US:New York"/)
res@mpGeophysicalLineThicknessF = 6
res@mpUSStateLineThicknessF = 6
res@mpNationalLineThicknessF = 6
res@mpOutlineDrawOrder   = "PreDraw"       ;      = "PostDraw" 
plot = gsn_csm_contour_map(wks,fillvarlevel,res)


;;;;;;;;;;;; this part tells the marker what color to be based on the anomaly
markercolors = new(size(0),"integer")
do n=0,size(0)-1,1
  markercolors(n) = round((siteanom(n)*15),3)+325
end do
markercolors=where((markercolors.lt.17),17,markercolors)
markercolors=where((markercolors.gt.255),237,markercolors)
stationmarker = new(size(0),graphic)
do n=0,size(0)-1,1
mres                = True         ; marker resources for best track
mres@gsMarkerIndex  = 16           ; marker style (filled circle)
mres@gsMarkerSizeF  = 23.0          ; marker size
mres@gsMarkerColor  = markercolors(n)      ; marker color
mres@tfPolyDrawOrder= "PostDraw"

stationmarker(n) = gsn_add_polymarker(wks,plot,sitelons(n),sitelats(n),mres)

tres                      = True                ; text mods desired
tres@tfPolyDrawOrder  = "PostDraw"
tres@txFont = "helvetica-bold"
tres@txFontHeightF        = 0.014
    if any(siteanom.lt.-9.5)
      tres@txFontHeightF        = 0.011
    end if
stationtext  =     gsn_add_text(wks,plot,round(siteanom(n),3),sitelons(n),sitelats(n),tres)
end do

tres@txFontHeightF        = 0.020
dum = gsn_add_text(wks,plot,"Mean Temperature Anomaly March 10~S~th~N~-12~S~th",-77.5,41.3,tres)
tres@txFontHeightF        = 0.015 
dum = gsn_add_text(wks,plot,"(data for all 3 days required)",-77.5,41.0,tres)


draw(plot)
frame(wks)


end