; 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

;---Read in site info

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

lonmin = -81 ;round(min(sitelons),3)-2
lonmax = -71 ;round(max(sitelons),3)+2
latmin = 38  ;round(min(sitelats),3)-2
latmax = 48  ;round(max(sitelats),3)+2

newlons = fspan(lonmin,lonmax,200)
newlats = fspan(latmin,latmax,200)

latsize=(dimsizes(newlats))
lonsize=(dimsizes(newlons))

newlons@units="degrees_east"
newlats@units="degrees_north"

;newdata = cssgrid(sitelats,sitelons,siteprecip,newlats,newlons)
;newdata = dsgrid2(sitelats,sitelons,siteprecip,newlats,newlons)
newdata = natgrid(sitelats,sitelons,siteprecip,newlats,newlons)
;newdata = obj_anal_ic(sitelats,sitelons,siteprecip,newlats,newlons,(/2,1,.5/),False)

;newdata = (newdata+newdata2)/2
newdata!0="lat"
newdata!1="lon"
newdata&lon = newlons
newdata&lat = newlats

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

res = True
res@gsnAddCyclic = False
res@mpProjection = "Mercator"
res@mpProjection = "CylindricalEquidistant"
res@mpLimitMode       = "Corners"
;res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
res@mpLeftCornerLatF  = 40.45
res@mpLeftCornerLonF  = -80
res@mpRightCornerLatF = 45.2
res@mpRightCornerLonF = -71.8
res@cnFillOn = True
res@cnFillMode = "RasterFill"
res@cnRasterSmoothingOn = True
res@cnLinesOn = False   ;no contours on top of fill
res@gsnMaximize = True
res@gsnFrame           = False 

  res@tmXBBorderOn      = True
  res@tmBorderLineColor = "white"
  res@tmXBMinorPerMajor = 0
  res@tmXBOn            = False
  res@tmXTBorderOn      = True
  res@tmXBMajorLineColor= "white"
  res@tmXTOn            = False
  res@tmYLBorderOn      = True
  res@tmYLMajorLineColor= "white"
  res@tmYLMinorPerMajor = 0
  res@tmYLOn            = False
  res@tmYRBorderOn      = True
  res@tmYRMajorLineColor= "white"
  res@tmYROn            = False

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

res@mpGeophysicalLineThicknessF = 8
res@mpUSStateLineThicknessF = 8
res@mpNationalLineThicknessF = 8
res@cnFillDrawOrder             = "PreDraw"
res@gsnDraw               = False  
res@cnLevelSelectionMode = "ManualLevels" ;-- set contour levels manually

gsn_define_colormap(wks,"BlAqGrYeOrReVi200")

colors = round(fspan(2,199,102),3)

res@cnMinLevelValF = 75 ; minvalue ;-- minimum contour level
res@cnMaxLevelValF = 95 ; maxvalue ;-- maximum contour level
res@cnLevelSpacingF = .2 ;-- contour level spacing
res@lbLabelStride = 25.0
res@lbOrientation = "Vertical"
res@cnFillColors = colors
res@lbBoxLinesOn = False

res@lbTitleString         = " "
res@lbTitleFontHeightF    = 0.015
res@lbLabelFontHeightF    = 0.01
res@pmLabelBarOrthogonalPosF = -0.18 
res@pmLabelBarParallelPosF   = 0.65
res@pmLabelBarHeightF        = .4

cres  = True
cres@cnLevelSelectionMode = "ExplicitLevels"
cres@cnLevels = (/75,80,85,90,95/)   
cres@cnLineThicknessF               = 10
cres@cnInfoLabelOn     = False
cres@cnLineColors      = "red"
cres@cnLineLabelFontColors  = "red"
cres@cnMonoLineLabelFontColor = True
cres@cnMonoLineThickness  = True
cres@cnMonoLineColor  = True
cres@cnMonoLevelFlag  = False
cres@cnLineLabelsOn   = True
cres@cnLineLabelBackgroundColor     = -1    ; transparent
cres@cnLineLabelPlacementMode   = "Constant"
cres@cnLineLabelInterval   =  1
cres@gsnAddCyclic = False
cres@cnLineDashSegLenF = .1
cres@cnLineDrawOrder             = "PreDraw"
plot = gsn_csm_contour_map_overlay(wks,newdata,newdata,res,cres)


draw(plot)

txres               = True                            
txres@txFontHeightF = .03          ; Set the font height
txres@txJust        = "CenterLeft"

label = "High Temperature (F)"
gsn_text_ndc(wks,label,.15,.30 ,txres)
frame(wks)

end