; 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"
;load "/ct12/abrammer/ncl6.3.0/lib/ncarg/nclscripts/csm/contributed.ncl"
begin

;set history count here, for number of days in period:
history = 365 ; 365 is the current maximum value. 

mindays = round((history*.94),3) ;; minimum number of days required to compute climate info

;; check to see if it's worth running whole script - no longer necessary
;recorddata = asciiread("records.ascii", -1, "string")
;sitedays = stringtofloat(recorddata(52))

;if sitedays.lt.mindays
;   exit()
;end if

    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


;;;;;;;;;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:5))

;;;;;; read the water data, and calculate distance from site if needed

delim = ","

if metadata(3).eq.0
    waterfile  = "water.dat"
    waterdata  = asciiread(waterfile,-1,"string")
    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("Minimum distance to water for "+meta(1)+":"+min(distance))
else
    distance = metadata(3)
end if

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

;sdtmax(0-4),sdtmin(5-9),tavg(10-14),tmax(15-19),tmin(20-24),ytdp(25-28/9)

temp=new((/8,365/), "float")
temp(0,:)=data(15,:)+metadata(2)*data(16,:)+metadata(0)*data(17,:)+metadata(1)*data(18,:)+min(distance)*data(19,:)      ; tmax
temp(1,:)=temp(0,:)+(data(0,:)+metadata(2)*data(1,:)+metadata(0)*data(2,:)+metadata(1)*data(3,:))+min(distance)*data(4,:)   ; tmax+
temp(2,:)=temp(0,:)-(data(0,:)+metadata(2)*data(1,:)+metadata(0)*data(2,:)+metadata(1)*data(3,:))+min(distance)*data(4,:)   ; tmax-
temp(3,:)=data(20,:)+metadata(2)*data(21,:)+metadata(0)*data(22,:)+metadata(1)*data(23,:)+min(distance)*data(24,:)      ; tmin
temp(4,:)=temp(3,:)+(data(5,:)+metadata(2)*data(6,:)+metadata(0)*data(7,:)+metadata(1)*data(8,:))+min(distance)*data(9,:)   ; tmin+
temp(5,:)=temp(3,:)-(data(5,:)+metadata(2)*data(6,:)+metadata(0)*data(7,:)+metadata(1)*data(8,:))+min(distance)*data(9,:)   ; tmin-
temp(6,:)=data(10,:)+metadata(2)*data(11,:)+metadata(0)*data(12,:)+metadata(1)*data(13,:)+min(distance)*data(14,:)        ; tavg
temp(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((/(history+12)/), "string")
yearstring=new((/(history+12)/), "string")
daynumber=new((/(history+12)/), "integer")
official=new((/4,(history+12)/), "float")
do i=((history*-1)-1),10,1
   daystring(i+(history+1)) = tostring(systemfunc("date +%j -d "+i+"days"))
   yearstring(i+(history+1))= tostring(systemfunc("date +%y -d "+i+"days"))
   daynumber(i+(history+1)) = tointeger(systemfunc("date +%j -d "+i+"days"))
end do

;;;;;; determine how many days are in this year
offset=365  ; normal year
first = ((history*-1)-1) ; set as minimum i above
if mod(tointeger(systemfunc("date +%y -d "+first+"days")),4.0).eq.0
    offset = 366
end if

;;;;;;reading in station data
official(:,:) = -996
do i=((history*-1)+1),0,1
   official(:,i+history) = asciiread(meta(1)+"/"+meta(1)+yearstring(i+history)+daystring(i+history)+".txt", -1, "float")
end do

;;;;;; getting right climo data
do i=0,(history),1
    if daynumber(i).gt.daynumber((history))
       daynumber(i)=daynumber(i)-offset
    end if
end do

;;;;;; getting right climo data
days = daynumber+offset

temps = table_attach_columns(temp,temp,0) ; make two years worth of data to handle transition across year
temps(7,365:729) = temps(7,365:729)+ temp(7,364) ; make sure precip is continuous

window = temps(:,(days(0)-1):(days((history+11))-1))

official@_FillValue = -996

overlap=new((/2,(history+12)/), "float")
print(dimsizes(overlap))
print(dimsizes(window))
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,:))

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

precipwindow = precip(:,(offset+daynumber(0)-1):(offset+daynumber((history+11))-1))
precipreal = precipwindow
if all(ismissing(official(2,:)))
    precipreal(1,:)=0
else
    precipreal(1,:)=cumsum(official(2,:),2)   ; plot only real data
end if

h = official(0,ind(.not.ismissing(official(0,:))))
l = official(1,ind(.not.ismissing(official(0,:))))

; find expected historical precip, but first check to make sure they exist
pindices=ind(.not.ismissing(official(2,:)))
if dimsizes(pindices).gt.1
   pexpvalid = new((/dimsizes(pindices)/), "float")
   pvalid = new((/dimsizes(pindices)/), "float")
   do xx = 0,(dimsizes(pindices)-1) 
       pexpvalid(xx) = precipwindow(1,pindices(xx))
       pvalid(xx) = official(2,pindices(xx))
   end do
   pflag = 0
else
  pexpvalid = 0.00
  pvalid    = 0.00
  pflag     = 1
end if

p =sum(pvalid)-sum(pexpvalid)
p =.01*round((100*p(0)),3)

; get various means
expectedhigh = avg(window(0,ind(.not.ismissing(official(0,:)))))
expectedlow  = avg(window(3,ind(.not.ismissing(official(0,:)))))
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)


indices=ind(.not.ismissing(official(0,:)))


;write out historical anomaly, if it exists
if dimsizes(indices).ge.round((history*.94),3)
    relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(actualmean-expectedmean)
    quote = inttochar(34) 
    system("echo "+quote+relative(0)+quote+" >> specified.dat")
    relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+actualmean
    system("echo "+quote+relative(0)+quote+" >> specifiedmean.dat")
    relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ avg(h)
    system("echo "+quote+relative(0)+quote+" >> specifiedhigh.dat")
    relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ avg(l)
    system("echo "+quote+relative(0)+quote+" >> specifiedlow.dat")

    relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ max(h)
    system("echo "+quote+relative(0)+quote+" >> specifiedmaximum.dat")
    relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ min(l)
    system("echo "+quote+relative(0)+quote+" >> specifiedminimum.dat")

    relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ (avg(h)-expectedhigh)
    system("echo "+quote+relative(0)+quote+" >> specifiedhighanomaly.dat")
    relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ (avg(l)-expectedlow)
    system("echo "+quote+relative(0)+quote+" >> specifiedlowanomaly.dat")
    monthhc = (h + l)/2
    heating = 0.0
    cooling = 0.0
    do i=0,(dimsizes(indices)-1)
        if (monthhc(i).gt.65)
           cooling=cooling+(monthhc(i)-65)
        else if (monthhc(i).lt.65)
           heating=heating+(monthhc(i)-65)
        end if
        end if
    end do
    monthd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+heating*-1
    system("echo "+quote+monthd(0)+quote+" >> specifiedheating.dat")
    monthd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+cooling
    system("echo "+quote+monthd(0)+quote+" >> specifiedcooling.dat")
end if


;;;;; precip output 
;write out precip historical anomaly and total, if it exists
if dimsizes(pindices).ge.round((history*.94),3)
    relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(p)
    mrelative =""+tofloat(meta(2))+","+tofloat(meta(3))+","+ sum(pvalid)
    quote = inttochar(34) 
    system("echo "+quote+relative(0)+quote+" >> precipspecified.dat")
    system("echo "+quote+mrelative(0)+quote+" >> preciptotalspecified.dat")
end if

;;;;; write out solar data
; write out historical average, if it exists
sindices=ind(.not.ismissing(official(3,:)))
solar = ((sum(official(3,:))*300)/1000000)/dimsizes(sindices)
if dimsizes(sindices).ge.round((history*.94),3)
    relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(solar)
    quote = inttochar(34) 
    system("echo "+quote+relative(0)+quote+" >> solarspecified.dat")
end if

end