#This program reads in the min and max raw data files and calculates annual average slope for each station#
#Final output files is "slopescomb0420" (Note:This has been commented out)#
#Other interim output files were written during program checkout, and are now commented out#
#Additional code has been added to take the average slope for each station & get a weighted average for USA lower 48#
#read in data#

loc="9641c_201012_raw.min" 

wd=c(6,1,4,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6)  

USAmin=read.fwf(loc,widths=wd)



ctdata=nrow(USAmin)
k=0
cnt=c()
cnt[k]=0
temp=c()
stano=c()
first=c()
last=c()

#Cnt[k] counts the line where each new station starts#

for (i in 2:(ctdata)) {
stacode1=USAmin[[i-1,1]]
stacode2=USAmin[[i,1]]
if (stacode2 != stacode1) {
k=k+1 
cnt[k]=i}
}
#I think I found the bug-the next to last line was not printing out#
#When i reaches ctdata, it stops and k does not get incremented#
#Therefore I will increment it once more#
#Update:That fixed the problem, in addition printing of first had to be changed to cnt[m-1]#
k=k+1
cnt[k]=i


#Define first line of table of station number, and first and last line#


stano[1] = USAmin[1,1]
first[1] = 1
last[1] = cnt[1]-1

#Define remaining table of  station numbers, and and first and last line#

for (m in 2:k){

stano[m] = USAmin[cnt[m-1],1]
first[m] = cnt[m-1]
last[m] = cnt[m]-1

}

#Define last line of table of  station numbers, and first and last line#


stano[m] = USAmin[ctdata,1]
first[m] = cnt[m-1]
last[m] = ctdata

datamin = data.frame(stano, first, last)
#write.table(datamin, file="dataminnew0420")#

ctdata1=nrow(datamin)

ctdata2=nrow(USAmin)

k=0
cnt[k]=0
m=0
n=0
slope=c()


#Read datamin back in as table#

#stainfo=read.table("datamin",header=FALSE)#
#ctdata1=nrow(stainfo)-1#

for (p in 1:ctdata1) {

#define location of first point for current station and next station#

m=datamin[p,2] #First point ofnrow current station#
n=datamin[p,3] #Last point of current station#

year=c(USAmin[m:n,3])


jan=c(USAmin[m:n,4])
feb=c(USAmin[m:n,6])
mar=c(USAmin[m:n,8])
apr=c(USAmin[m:n,10])
may=c(USAmin[m:n,12])
jun=c(USAmin[m:n,14])
jul=c(USAmin[m:n,16])
aug=c(USAmin[m:n,18])
sep=c(USAmin[m:n,20])
oct=c(USAmin[m:n,22])
nov=c(USAmin[m:n,24])
dec=c(USAmin[m:n,26])
avg=c(USAmin[m:n,28])

stot=0
#cat(mar,file="temp",sep=" ",append=TRUE,"\n")#

temp=data.frame(jan,feb,mar,apr,may,jun,jul,aug,sep,oct,nov,dec,avg)



for (j in 1:13) {
for (i in 1:(n-m+1)) {
if (temp[[i,j]] ==-9999) temp[i,j]<- NA else temp[i,j] = temp[i,j]/10
}
}



k=13

y=temp[[k]]
c=coef(lm(y~year))
slope[p]=c[2]


}


slopes=data.frame(datamin[ ,1],slope)
#write.table(slopes, file="slopes3new0420")#

#read in Max data#

loc="9641c_201012_raw.max" 

wd=c(6,1,4,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6,1,6)  

USAmax=read.fwf(loc,widths=wd)



ctdata=nrow(USAmax)
k=0
cnt=c()
cnt[k]=0
temp=c()
stanomax=c()
firstmax=c()
lastmax=c()


#Cnt[k] counts the line where each new station starts#

for (i in 2:(ctdata-1)) {
stacode1=USAmax[[i-1,1]]
stacode2=USAmax[[i,1]]
if (stacode2 != stacode1) {
k=k+1
cnt[k]=i}
}
l=i

#When i reaches ctdata, it stops and k does not get incremented#
#Therefore I will increment it once more#
k=k+1
cnt[k]=i

#Define first line of table of station number, and first and last line#


stanomax[1] = USAmax[1,1]
firstmax[1] = 1
lastmax[1] = cnt[1]-1



#Write define table of  station numbers, and and first and last line#

for (m in 2:k) {
if (cnt[m]-cnt[m-1] >=20) {

stanomax[m] = USAmax[cnt[m-1],1]
firstmax[m] = cnt[m-1]
lastmax[m] = cnt[m]-1
}
}

#Define the last line of file of  station number, and first and last line#


stanomax[m] = USAmax[ctdata,1]
firstmax[m] = cnt[m-1]
lastmax[m] = ctdata
datamax = data.frame(stanomax, firstmax, lastmax)
#write.table(datamax, file="datamaxnew0420")#


k=0
cnt=c()
cnt[k]=0
tempmx=c()

slopemax=c()


#Read datamax back in as table#

#stainfo=read.table("datamax",header=FALSE)#
ctdata1=nrow(datamax)

for (p in 1:ctdata1) {

#define location of first point for current station and next station#

m=datamax[p,2] #First point ofnrow current station#
n=datamax[p,3] #Last point of current station#

year=c(USAmax[m:n,3])


janmax=c(USAmax[m:n,4])
febmax=c(USAmax[m:n,6])
marmax=c(USAmax[m:n,8])
aprmax=c(USAmax[m:n,10])
maymax=c(USAmax[m:n,12])
junmax=c(USAmax[m:n,14])
julmax=c(USAmax[m:n,16])
augmax=c(USAmax[m:n,18])
sepmax=c(USAmax[m:n,20])
octmax=c(USAmax[m:n,22])
novmax=c(USAmax[m:n,24])
decmax=c(USAmax[m:n,26])
avgmax=c(USAmax[m:n,28])

tempmx=data.frame(janmax,febmax,marmax,aprmax,maymax,junmax,julmax,augmax,sepmax,octmax,novmax,decmax,avgmax)

#cat(mar,file="tempmax",sep=" ",append=TRUE,"\n")#





for (j in 1:13) {
for (i in 1:(n-m+1)) {
if (tempmx[[i,j]] ==-9999)
tempmx[i,j]<- NA else tempmx[i,j] = tempmx[i,j]/10

}
}


k=13
y=tempmx[[k]]
c=coef(lm(y~year))

slopemax[p]=c[2]

}


slopesmax=data.frame(datamax[ ,1],slopemax)
#write.table(slopesmax, file="slopesmax3new0420")#

slopescomb=data.frame(datamax[ ,1],slope, slopemax)
#write.table(slopescomb, file="slopescomb0420")#

loc="ushcn-stations.txt"
wd=c(6,9,10,8,2,31)
Stations=read.fwf(loc,widths=wd)

ctlines=0
state1=c()
state2=c()


ctlines=nrow(Stations)
k=0
cnt=c()
first=c()
last=c()


#cnt() counts the line where each state starts#

for (i in 2:(ctlines)) {
state1=Stations[[i-1, 5]]
state2=Stations[[i,5]]
if (state2 != state1) {
k=k+1
cnt[k]=i}
}

#For the last line#
k=k+1
cnt[k]=i


#define the first line where the first state starts#

state = Stations[1,5]
first[1]=1
last[1]=cnt[1]-1


#Define the remaining list of states and first and last lines#

for (m in 2:k){
state[m]=Stations[cnt[m-1], 5]
first[m]= cnt[m-1]
last[m]=cnt[m]-1
}
#Define the last state and first and last lines#

state[m] = Stations[ctlines,5]
first[m] = cnt[m-1]
last[m]= ctlines

statedata = data.frame(state, first, last)
#write.table(statedata, file = "statedata0906")#



#The stations will be averaged by states#

summin=0
summax=0
p=0
r=0
s=0
t=0
x=0
avgmin=c()
avgmax=c()
 
for (r in 1:48){
s=first[r]
t=last[r]
for (p in s:t){
summin=summin+slope[p]
summax=summax+slopemax[p]
}
x=x+1
avgmin[x]=summin/(t-s+1)
avgmax[x]=summax/(t-s+1)
summin=0
summax=0
}

stateavg=data.frame(state,avgmin, avgmax)
#write.table(stateavg,file="stateaverages")#

starea=c()
starea=read.table(file="StateAreas.txt")

#Calculate area weighted averages#

x=0
y=0
summinxa=0
summaxxa=0
suma=0

y=nrow(stateavg)

for (x in 1:y){
summinxa=summinxa+(stateavg[x,2])*(starea[x,2])
summaxxa=summaxxa+(stateavg[x,3])*(starea[x,2])
suma=suma+starea[x,2]
}
wtmin=summinxa/suma
wtmax=summaxxa/suma