; 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  = "maxwc2day.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 (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,"BlAqGrYeOrReVi200")

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

minvalue = minvalue - where(minvalue.ge.0,mod(minvalue,5),mod(5-abs(minvalue),5))-where(minvalue.le.-5,5,0)
maxvalue = ceil((5-mod(maxvalue,5))+maxvalue)

diff = toint((maxvalue-minvalue))+2

if diff.gt.99
colors = round(fspan(2,199,diff),3)
res@cnLevelSpacingF = 1.0 ;-- contour level spacing
res@lbLabelStride = 10.0
end if

if (diff.lt.99).and.(diff.ge.49)
colors = round(fspan(2,199,2*diff),3)
res@cnLevelSpacingF = 0.5 ;-- contour level spacing
res@lbLabelStride = 10.0
end if

if diff.lt.49
colors = round(fspan(2,199,4*diff),3)
res@cnLevelSpacingF = 0.25 ;-- contour level spacing
res@lbLabelStride = 20.0
end if


res@cnMinLevelValF = minvalue ;-- minimum contour level
res@cnMaxLevelValF = 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 diff.ge.25
    cres@cnLevels = ispan(-70,120,10)
end if
if diff.lt.25
    cres@cnLevels = ispan(-70,120,5)
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 = "Minimum Wind Chill (~S~o~N~F)"
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