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

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

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

;newdata = cssgrid(sitelats,sitelons,siteprecip,newlats,newlons)
newdata2 = 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,"weeklyheatingmap")

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"/)
res@mpDataBaseVersion = "MediumRes" ;-- better map resolution

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

;gsn_define_colormap(wks,"MPL_YlGn")

colors = ispan(5,255,1)

res@cnMinLevelValF =  0 ; minvalue ;-- minimum contour level
res@cnMaxLevelValF =  500 ; maxvalue ;-- maximum contour level
res@cnLevelSpacingF = 2 ;-- contour level spacing
res@lbLabelStride = 25
res@lbOrientation = "Vertical"
res@cnFillColors = colors
res@lbBoxLinesOn = False

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

hghtfile  = "topo15.dat"
hghtdatafile  = asciiread(hghtfile,-1,"string")
hghtlats  = tofloat(str_get_field(hghtdatafile,3,delim))
hghtlons  = tofloat(str_get_field(hghtdatafile,4,delim))
hght      = tofloat(str_get_field(hghtdatafile,5,delim))*3.28084

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

hghtdata = (hghtdata+hghtdata2)/2
hghtdata!0="lat"
hghtdata!1="lon"
hghtdata&lon = newlons
hghtdata&lat = newlats

cres  = True
cres@cnLevelSelectionMode = "ExplicitLevels"
cres@cnLevels = (/500,750,1000,1250,1500,1750,2000,2250,2500/)   
;cres@cnLineThicknesses               = (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21/)
cres@cnLineThicknesses               = (/5,6,7,8,9,10,11,12,13,14,15,20,22/)
cres@cnInfoLabelOn     = False
cres@cnLineColors      = (/"gray2","gray10","gray18","gray26","gray34","gray42","gray50","gray58","gray64","gray72"/)
cres@cnLineLabelFontColors  = "red"
cres@cnMonoLineLabelFontColor = True
cres@cnMonoLineThickness  = False
cres@cnMonoLineColor  = False
cres@cnMonoLevelFlag  = False
cres@cnLineLabelsOn   = True
cres@cnLineLabelBackgroundColor     = -1    ; transparent
cres@cnLineLabelPlacementMode   = "Constant"
cres@cnLineLabelInterval   =  0
cres@gsnAddCyclic = False
cres@cnLineDashSegLenF = .3
cres@cnLineDrawOrder             = "PreDraw"
plot = gsn_csm_contour_map_overlay(wks,newdata,hghtdata,res,cres)

;mres                = True         ; marker resources for best track
;mres@gsMarkerIndex  = 16           ; marker style (filled circle)
;mres@gsMarkerSizeF  = 5.0          ; marker size
;mres@gsMarkerColor  = "gray"       ; marker color
;markers = gsn_add_polymarker(wks,plot,sitelons,sitelats,mres)

draw(plot)

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

label = "Weekly Heating Degree Days (HDD, fill) and Elevation"
gsn_text_ndc(wks,label,.03,.34 ,txres)

label = "Time Period: "+systemfunc("date -u -d -7days +%a")+"., "+systemfunc("date -u -d -7days +%F")+" through "+ systemfunc("date -u -d -1days +%a")+"., "+systemfunc("date -u -d -1days +%F")
gsn_text_ndc(wks,label,.03,.30 ,txres)

label = (dimsizes(sitelats))+" Stations Were Used To Make This Plot"
gsn_text_ndc(wks,label,.03,.26 ,txres)

frame(wks)

end