; 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  = "snow2day.dat"
sitedatafile = asciiread(sitefile,-1,"string")
delim = ","
sitelats  = tofloat(str_get_field(sitedatafile,1,delim))
sitelons  = tofloat(str_get_field(sitedatafile,2,delim))
sitedata  = tofloat(str_get_field(sitedatafile,3,delim))

;;; this part adds fake stations around the periphery of NY so that the map is smoother near the border
size = dimsizes(sitelats)-1
do i = 1,size
      if sitelats(i).gt.44.87
          sitelats := array_append_record(sitelats,sitelats(i)+1,0)
          sitelons := array_append_record(sitelons,sitelons(i),0)
          sitedata := array_append_record(sitedata,sitedata(i),0)
      end if
      if (sitelats(i).lt.42.1) .and. (sitelons(i).lt.(-75.9))
          sitelats := array_append_record(sitelats,sitelats(i)-1,0)
          sitelons := array_append_record(sitelons,sitelons(i),0)
          sitedata := array_append_record(sitedata,sitedata(i),0)
      end if
      if (sitelats(i).gt.41.5) .and. (sitelons(i).gt.(-73.6))
          sitelats := array_append_record(sitelats,sitelats(i),0)
          sitelons := array_append_record(sitelons,sitelons(i)+1,0)
          sitedata := array_append_record(sitedata,sitedata(i),0)
      end if
      if (sitelons(i).lt.(-78.74))
          sitelats := array_append_record(sitelats,sitelats(i),0)
          sitelons := array_append_record(sitelons,sitelons(i)-1,0)
          sitedata := array_append_record(sitedata,sitedata(i),0)
      end if
end do


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,400)
newlats = fspan(latmin,latmax,400)

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

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

;newdata = cssgrid(sitelats,sitelons,sitedata,newlats,newlons)
;newdata2 = dsgrid2(sitelats,sitelons,sitedata,newlats,newlons)
newdata = natgrid(sitelats,sitelons,sitedata,newlats,newlons)
;newdata = obj_anal_ic(sitelats,sitelons,sitedata,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,"plot")

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 = 9
res@mpUSStateLineThicknessF = 9
res@mpNationalLineThicknessF = 9
res@mpCountyLineThicknessF = 2.5

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

gsn_define_colormap(wks,"wh-bl-gr-ye-re")

maxvalue = max(sitedata)
minvalue = 0
;minvalue = min(floor(sitedata))

;minvalue = minvalue - mod(minvalue,5)
;maxvalue = ceil((5-mod(maxvalue,5))+maxvalue)

if ceil(maxvalue).le.4
colors := round(fspan(2,199,toint(maxvalue*15)+5),3)
res@cnMinLevelValF = .5 ;-- minimum contour level
res@cnLevelSpacingF = 0.1
res@lbLabelStride = 5.0
end if

if (ceil(maxvalue).lt.10).and.(ceil(maxvalue).ge.4)
colors := round(fspan(2,199,toint(maxvalue*4)+10),3)
res@cnMinLevelValF = 1 ;-- minimum contour level
res@cnLevelSpacingF = 0.2
res@lbLabelStride = 5.0
end if

if (ceil(maxvalue).ge.10).and.(ceil(maxvalue).lt.25)
colors := round(fspan(2,199,toint(4*maxvalue)+2),3)
res@cnMinLevelValF = 1 ;-- minimum contour level
res@cnLevelSpacingF = 0.25
res@lbLabelStride = 16.0
end if

if ceil(maxvalue).ge.25
colors := round(fspan(2,199,toint(2*maxvalue)+2),3)
res@cnMinLevelValF = 2 ;-- minimum contour level
res@cnLevelSpacingF = 0.5
res@lbLabelStride = 8.0
end if

;colors := array_append_record(0,colors,0)

res@cnMaxLevelValF = (ceil(maxvalue)) ;-- maximum contour level
res@lbOrientation = "Vertical"
res@cnFillColors = colors
res@lbBoxLinesOn = False

res@lbTitleString         = " "
res@lbTitleFontHeightF    = 0.012
res@lbLabelFontHeightF    = 0.01
res@pmLabelBarOrthogonalPosF = -0.20 
res@pmLabelBarParallelPosF   = 0.65
res@pmLabelBarHeightF        = .3

cres  = True
cres@cnLevelSelectionMode = "ExplicitLevels"

if (ceil(maxvalue)).lt.4
    cres@cnLevels = .5*ispan(1,10,1) 
end if
if ((ceil(maxvalue)).ge.4) .and. ((ceil(maxvalue)).le.10)
    cres@cnLevels := ispan(1,11,1)
end if
if (ceil(maxvalue).ge.10).and.(ceil(maxvalue).lt.25)
    cres@cnLevels := ispan(3,27,3)
end if
if (ceil(maxvalue)).ge.25 .and.(ceil(maxvalue).lt.50)
    cres@cnLevels := ispan(5,800,5)
end if
if (ceil(maxvalue)).ge.50
    cres@cnLevels := ispan(10,800,10)
end if

cres@cnLineThicknessF               = 12.5
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        = "CenterCenter"

label = "Total Snowfall (in)"
gsn_text_ndc(wks,label,.32,.35 ,txres)

txres@txFontHeightF = .022 
label = "Days Of: " + systemfunc("date -u -d -2days +%F") + "  And  " + systemfunc("date -u -d -1days +%F")
gsn_text_ndc(wks,label,.32,.29 ,txres)

frame(wks)

end