; 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

    fasc = "../*txt"   ; a unique identifier for file
    
    DASC = "./"         ; output dir
    FASC = "BIG.txt"    ; output file name

    system ("/bin/rm -f "+DASC+FASC)   ; rm any pre-existing file
    ncol = stringtoint(systemfunc("ls -1 ../*txt | wc -l"))

; Use UNIX "cat" to concatenate the files into one file.
    system ("cat "+fasc+" > "+DASC+FASC)

; You can now read the file via "asciiread".

    data  = asciiread(DASC+FASC,(/ncol,365/),"float")

    system ("/bin/rm -f "+DASC+FASC)   ; rm any pre-existing file

    fils = systemfunc ("ls *txt")  

;;;;;;;;;Here's the order of data

;sdtmax(0-3),sdtmin(4-7),tavg(8-11),tmax(12-15),tmin(16-19),ytdp(20-23)

; meta files with info

meta = asciiread("station.ascii", -1, "string")
metadata = stringtofloat(meta(2:4))

waterfile  = "water.dat"
waterdata = asciiread(waterfile,-1,"string")
delim = ","

waterlats    =  tofloat(str_get_field(waterdata,3,delim))
waterlons    =  tofloat(str_get_field(waterdata,4,delim))
wateryes     =  tofloat(str_get_field(waterdata,5,delim))

distance =  new((/dimsizes(waterlats)/),float)
distance = 10000
do i = 0,dimsizes(waterlats)-1
distance(i) = gc_latlon(metadata(0),metadata(1),waterlats(i),waterlons(i),2,4)
distance(i) = distance(i)*wateryes(i)
end do
distance@_FillValue = 0
print(min(distance))

;days      =ispan(1,365,1)

temps=new((/8,365/), "float")
temps(0,:)=data(15,:)+metadata(2)*data(16,:)+metadata(0)*data(17,:)+metadata(1)*data(18,:)+min(distance)*data(19,:)      ; tmax
temps(1,:)=temps(0,:)+(data(0,:)+metadata(2)*data(1,:)+metadata(0)*data(2,:)+metadata(1)*data(3,:))+min(distance)*data(4,:)   ; tmax+
temps(2,:)=temps(0,:)-(data(0,:)+metadata(2)*data(1,:)+metadata(0)*data(2,:)+metadata(1)*data(3,:))+min(distance)*data(4,:)   ; tmax-
temps(3,:)=data(20,:)+metadata(2)*data(21,:)+metadata(0)*data(22,:)+metadata(1)*data(23,:)+min(distance)*data(24,:)      ; tmin
temps(4,:)=temps(3,:)+(data(5,:)+metadata(2)*data(6,:)+metadata(0)*data(7,:)+metadata(1)*data(8,:))+min(distance)*data(9,:)   ; tmin+
temps(5,:)=temps(3,:)-(data(5,:)+metadata(2)*data(6,:)+metadata(0)*data(7,:)+metadata(1)*data(8,:))+min(distance)*data(9,:)   ; tmin-
temps(6,:)=data(10,:)+metadata(2)*data(11,:)+metadata(0)*data(12,:)+metadata(1)*data(13,:)+min(distance)*data(14,:)        ; tavg
temps(7,:)=.1*(data(25,:)+metadata(2)*data(26,:)+metadata(0)*data(27,:)+metadata(1)*data(28,:)) ;  +min(distance)*data(29,:) ; precip - use last part if you have fixed problem

; figure out what days to plot

daystring=new((/86/), "string")
daynumber=new((/86/), "integer")
official=new((/3,86/), "float")
do i=-75,10,1
   daystring(i+75) = tostring(systemfunc("date +%j -d "+i+"days"))
   daynumber(i+75) = tointeger(systemfunc("date +%j -d "+i+"days"))
end do
print(tointeger(systemfunc("date +%j -d "+0+"days")))
;reading in station data
official(:,:) = -996
do i=-73,0,1
   official(:,i+74) = asciiread(meta(1)+daystring(i+74)+".txt", -1, "float")
end do

; getting right climo data
days = daynumber
window = temps(:,(daynumber(0)-1):(daynumber(85)-1))

official@_FillValue = -996

wks_type = "png"
wks_type@wkWidth = 2048
wks_type@wkHeight = 2048 

wks   = gsn_open_wks(wks_type,meta(1))

res                    = True
res@gsnFrame      = False
res@tiMainFontHeightF  = 0.018
res@tiXAxisFontHeightF = 0.015
res@tiYAxisFontHeightF = 0.015
res@tiXAxisString     = "Month"
res@tmXBTickSpacingF  = 1.
res@tmXBTickStartF    = 0.
res@tmXMajorGrid      = True
res@tmYLMinorPerMajor = 1
res@tmYMajorGrid      = True
res@trXMinF           = days(0)
res@trXMaxF           = days(85)
res@tmXBTickSpacingF  = 5
res@tmXBMode          = "Explicit"  			; explicitly set X-axis labels
res@tmXBValues        = (/"1","32","60","91","152","182","213","244","274","305","335"/)
res@tmXBLabels        = (/"Jan","Feb","Mar","Apr","Jun","Jul","Aug","Sep","Oct","Nov","Dec","Jan"/)
res@tmXBLabelFontHeightF = .015
res@tmYROn               = False 


res@tiMainString      = "Recent Climate Values For "  + meta(0)       ; add title
res@tiYAxisString     = "Average Daily High/Low/Mean (Red/Blue/Green) +/- 1 SD"
res@xyLineColors      = (/"red","red","red","blue","blue","blue","green","green"/)          ; change line color
res@xyLineThicknesses = (/30.0,8.0,8.0,24.0,4.0,4.0,24.0,24.0/) 
res@tiYAxisFontHeightF = 0.015
res@trYMinF           =  min(window)-30
res@trYMaxF           =  max(window)+10
res@tmYLPrecision     = 2


;res@xyMonoMarkLineMode= False
;res@xyMarkLineModes    = (/"Lines","Lines","Lines","Markers","Markers","Markers","Markers"/)
res@xyMarkLineModes    = "Lines"
;res@xyMonoDashPattern = False
;res@xyDashPattern     = (/"SolidLine","SolidLine","SolidLine",4,4,4,4/)
res@xyDashPattern     = (/0,0,0,4,4,4,4,4/)
;res@xyMarkers = 16 ; some dot 
;res@xyMarkerColors = "red"
;res@xyMarkerSizes = 0.055


res@gsnXYFillColors =   180 ; "Light Orange"
res@xyLineColor     = -1 
plot  = gsn_csm_xy(wks,days,window(0:2,:),res); create plot

res@xyLineColors      = (/"blue","blue","blue","blue","blue","blue","green","green"/)          ; change line color

res@gsnXYFillColors = 50 ; "pink1"
res@xyLineColor     = -1 
plot  = gsn_csm_xy(wks,days,window(3:5,:),res); create plot

overlap=new((/2,86/), "float")
overlap(0,:)=window(2,:)
overlap(1,:)=window(4,:)

overlap(0,:)=where(overlap(0,:).gt.overlap(1,:),overlap(1,:),overlap(0,:))
overlap(1,:)=where(overlap(0,:).gt.overlap(1,:),overlap(0,:),overlap(1,:))

res@xyLineThicknesses = (/0.001,0.001,1,24,4,4,24,24/) 
res@xyLineColors      = (/"blue","red","blue","blue","blue","blue","green","green"/)          ; change line color
res@gsnXYFillColors = 255 ; "pink1"
res@xyLineColor     = -1 
plot  = gsn_csm_xy(wks,days,overlap,res); create plot

precip = new((/2,365/), "float")
precip(1,:) = temps(7,:)
do i=1,364,1
   precip(0,i) = 0
   precip(1,i) = temps(7,i)-temps(7,i-1)
end do

res@xyLineThicknesses = (/40.0,8.0,8.0,24.0,4.0,4.0,24.0,24.0/) 
res@xyLineColors      = (/"green","blue","blue","blue","blue","blue","green","green"/)  
res@gsnXYFillColors = 255 ; "pink1"
res@xyLineColor     = -1 
res@tiMainString      = ""          ; add title

res2 = res
res2@tmYUseLeft = False
res2@xyLineThicknesses = (/0.01,0.01,8.0,24.0,4.0,4.0,24.0,24.0/) 
res2@xyLineColors      = (/"black","black","blue","blue","blue","blue","green","green"/)  
res2@gsnXYFillColors = 85 ; "pink1"
res2@tiYAxisString     = "Expected Daily Precipitation (in)"
res2@tmYMajorGrid      = False
res2@xyMarkLineModes   = "Lines"                ; choose which have markers
res2@trYMinF           = 0
res2@trYMaxF           = max(precip)*5
res2@tmYMinorPerMajor  = 0
res2@tmYRLabelsOn      = True                ; turn on the YR tick mark labels
res2@tmYRMode          ="Explicit"
res2@tmYRValues        = (/".01",".05",".10",".15"/)
res2@tmYRLabels        = (/".01",".05",".10",".15"/)
res2@tmYROn            = True 
res2@tmYRLabelJust     = "CenterLeft"
precipwindow = precip(:,(daynumber(0)-1):(daynumber(85)-1))
plot  = gsn_csm_xy2(wks,days,window(6,:),precipwindow,res,res2) ; create plot

delete(res)

res                    = True
res@gsnFrame      = False
res@tiMainFontHeightF  = 0.018
res@tiXAxisFontHeightF = 0.015
res@tiYAxisFontHeightF = 0.015
res@tmBorderThicknessF = 40
res@tmXBTickSpacingF  = 1.
res@tmXBTickStartF    = 0.
res@tmXMajorGrid      = True
res@tmYLMinorPerMajor = 1
res@tmYMajorGrid      = False
res@trXMinF           = days(0)
res@trXMaxF           = days(85)
res@tmXBTickSpacingF  = 5
res@tmXBMode          = "Explicit"  			; explicitly set X-axis labels
res@tmXBValues        = days(75)
res@tmXBLabels        = "Today"
res@tmXBLabelFontHeightF = .015
res@tmYROn               = False 
res@xyLineColors      = (/"black","black","red","blue","blue","blue","green","green"/)          ; change line color
res@xyLineThicknesses = (/15.0,15.0,8.0,24.0,4.0,4.0,24.0,24.0/)     ; this seems to control the bar chart too!
res@tiYAxisFontHeightF = 0.015
res@trYMinF           =  min(window)-30
res@trYMaxF           =  max(window)+10
res@tmYLPrecision     = 2
res@tmXMajorGridThicknessF = 10

resbar = res
resbar@gsnXYBarChart         = True
resbar@gsnXYBarChartBarWidth = .04      ; Default is 1.0
resbar@gsnYRefLineColor      = -1          ; reference line   
resbar@gsnXYBarChart         = True            ; create bar chart 
resbar@gsnAboveYRefLineColor = "black"           ; above ref line fill red
resbar@gsnBelowYRefLineColor = "black"          ; below ref line fill blue
plot1 = new((/2,days(85)/),graphic)

; get indices of value data, and only plot those
indices=ind(.not.ismissing(official(0,:)))
do xx = 0,(dimsizes(indices)-1)   
     gg = indices(xx)
     resbar@gsnXYBarChartColors = "black"
     resbar@gsnYRefLine         =  (official(0,gg)+official(1,gg))/2
     plot1(0,gg) = gsn_csm_xy(wks,days(gg),official(0,gg),resbar)    ; draw each individual max bar
     plot1(1,gg) = gsn_csm_xy(wks,days(gg),official(1,gg),resbar)    ; draw each individual min bar
end do

res@tmBorderThicknessF = 0
res@xyMarkers         = 16                  ; choose type of marker  
res@xyMarkerColor     = "Black"            ; Marker color
res@xyMarkerSizeF     = 0.015   
res@xyMarkLineMode    = "Markers"             ; choose to use markers
plot  = gsn_csm_xy(wks,days,official(0,:),res)                    ; create plot
plot  = gsn_csm_xy(wks,days,official(1,:),res)                    ; create plot

res@xyMarkerColor     = "Red"            ; Marker color
res@xyMarkerSizeF     = 0.0075                     ; Marker size (default 0.01)

plot  = gsn_csm_xy(wks,days,official(0,:),res)                    ; create plot
res@xyMarkerColor     = "Blue"            ; Marker color
plot  = gsn_csm_xy(wks,days,official(1,:),res)                    ; create plot

; print extremes, etc.

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

label = "Location: " + meta(0)
gsn_text_ndc(wks,label,.09,.11 ,txres)

label = "Station elevation: " + .01*round(100*tofloat(meta(4)),3) + " m"
gsn_text_ndc(wks,label,.09,.09 ,txres)

label = "Station lat/lon: " + .01*round(100*tofloat(meta(2)),3) + "; " + .01*round(100*tofloat(meta(3)),3)
gsn_text_ndc(wks,label,.09,.07 ,txres)

label = "Figure created on " + systemfunc("date -u +%a") + ", " + systemfunc("date -u +%F") + ", at "+systemfunc("date -u +%R") + " UTC"
gsn_text_ndc(wks,label,.09,.05 ,txres)

; get actual availably highs/lows
h = official(0,ind(.not.ismissing(official(0,:))))
l = official(1,ind(.not.ismissing(official(0,:))))

; find expected precip, but first check to make sure they exist

if dimsizes(indices).gt.1
   pexpvalid = new((/dimsizes(indices)-1/), "float")
   pvalid = new((/dimsizes(indices)-1/), "float")
   do xx = 0,(dimsizes(indices)-1) 
       i = 0
       pexpvalid(i) = precipwindow(1,indices(xx))
       pvalid(i) = official(2,indices(xx))
       i = i+1
   end do
else
  pexpvalid = 0
  pvalid    = 0
end if

p =.01*round(1000*sum(pvalid),3)

; get various means
expectedmean = avg(window(6,ind(.not.ismissing(official(0,:)))))
actualmean   = (avg(h)+avg(l))/2
d = .1*round(10*(actualmean - expectedmean),3)
;w =.1*round(10*max(winds(1,:)),3)
;t =.01*round(100*max(cumsum(precip,2)),3)

txres               = True                            
txres@txFontHeightF = .014          ; Set the font height
txres@txJust        = "CenterRight"
label = "Summary Of Past "+dimsizes(indices)+" Days"
gsn_text_ndc(wks,label,.9,.11 ,txres)

label = "Warmest/Coldest:"
gsn_text_ndc(wks,label,.8,.09 ,txres)

label = "Versus Expected:"
gsn_text_ndc(wks,label,.8,.07 ,txres)

label = "(Mean Temperature/Precipitation)"
gsn_text_ndc(wks,label,.9,.05 ,txres)

txres@txFont = "helvetica-bold"
label = .1*round(10*max(official(0,:)),3)+"/"+ .1*round(10*min(official(1,:)),3)
gsn_text_ndc(wks,label,.9,.09 ,txres)

txres@txFont = "helvetica-bold"
txres@txFontColor = "blue"
text = ""
if d.ge.0
   text = "+"
   txres@txFontColor = "red"
end if

label = text(0)+d+"/"+"?"
gsn_text_ndc(wks,label,.9,.07 ,txres)

frame(wks)

end