AESINDAT='/home/jallen/aes_data/1154410.001.out'
REANINDAT='/home/jallen/reanal_data/daily/1154410.001'

AESINDAT='/home/jallen/aes_data/1096450.001.out'
REANINDAT='/home/jallen/rdata_modified/1096450.001'

PLOTOUT='PG_JA.ps.test'
STATOUT='PG_JA.stat.test'
DATAOUT="PG_JA.dat.test"
HREANLAB='Reanalysis Temperature at 00Z(oC)'
HAESLAB='Station Maximum Daily Temperature (oC)'
HTITLE="Comparison of Reanalysis and Station Data, Maximum Daily Temperature, PG"
BYLAB='Temperature (oC)'
BXLAB='Station Data; Reanalysis Data'
BTITLE="Station and Reanalysis Data, Maximum Daily Temperature, PG"

#collumn contents, if different from above*

#---------------------Start of Program--------------------------------

#Read the files AESINDAT and REANINDAT, which contain station (aes) and
#reanalysis data for the station and element indicated; assign this information 
#to the arrays 'aes' and 'rean'.

outa<-"Reading data files.";print(outa,quote=F)

aes<-read.table(AESINDAT, header=F)
rean<-read.table(REANINDAT, header=F)


#initialize counter variables
i<-1;k<-1;s<-0;numdat<-0;nummissing<-0;totaldat<-0

outb<-"Searching for dates in July and August.";print(outb,quote=F)
#This loop locates the relevant data and writes it to the file DATAOUT

i<-1;k<-1;s<-0;numdat<-0;nummissing<-0;totaldat<-0

while(i<=length(aes[,1])){
 if(aes[i,2]==7||aes[i,2]==8){
  if(aes[i,4]*0.1!=-99.9){
  print(i)
   while(s!=1){
#      print(k)
      if(rean[k,1]==aes[i,1]){
if(rean[k,2]==aes[i,2]){
if(rean[k,3]==aes[i,3]){
out<-array(c(aes[i,1],aes[i,2],aes[i,3],aes[i,4]*0.1, (rean[k,5]-273.1)), dim=c(1,5))
write(out, file=DATAOUT, append=T)
numdat<-numdat+1
s<-1
}
}
}
k<-k+1 
if(k==length(rean[,1])){

s<-1
k=1
} 
}
s<-0
}else{nummissing<-nummissing+1}
  totaldat<-totaldat+1
 }
 i<-i+1
}


#Read the file DATAOUT; store this information in the array 'res'.
res<-read.table(DATAOUT, header=F)

#Create a scatter plot of reanalysis data as a function of station
#data, and boxplots of station and reanalysis data.  Save the
#resultant graphs to the file PLOTOUT, in postscript format.
 
postscript(PLOTOUT)
plot(res[,4],res[,5], ylab=HREANLAB, xlab=HAESLAB, main=HTITLE)
boxplot(res[,4], res[,5], ylab=BYLAB, xlab=BXLAB, main=BTITLE)
dev.off()

#Perform a linear regression on reanalysis data as a function of 
#station data.  Save a summary of the results to the file STATOUT.

sink(STATOUT)
out1<-"The number of non-missing data points is"
out2<-"The number of missing data points is"
print(out1,quote=F);print(numdat)
print(out2,quote=F);print(nummissing)

comp.lm<-lm(res[,5] ~ res[,4])
sum<-summary(comp.lm)
print(sum)
sink()

CompJ:

ISOLATE='inv_rean.dates'
SDATA='inv_rean.date'
OUTDAT='compositePrec.dat'
OUTTALLY="/home/jallen/compare/PrinceGeorge/compositePrec.tal"
OUTPLOT="/home/jallen/compare/PrinceGeorge/compositePrec.ps"
TITLE='Time series Composite for Daily Precipitation, Prince George'
YLAB="Precipitation (mm)"
XLAB='D-3            D-2            D-1            D            D+1            D+2            D+3'
ycol<-3
mcol<-4
dcol<-5
#---------------------Start of Program--------------------------------

#read in a table containing the desired dates, in PST.  The year, month, and
#day are found in collumns ycol, mcol, and dcol, respectively; this
#information is saved in the array 'iso'. 

iso<-read.table(ISOLATE, header=FALSE)
out1<-"Reading input data";print(out1, quote=F)


#read station data file; this data is saved as the array sdata.
sdata<-read.table(SDATA, header=FALSE)

#set counter variables
i<-1;j<-1;numdata<-0;ok<-1

#find the first entry in iso for which station data is present
while(iso[j,ycol]<sdata[1,1]){j<-j+1}
out2<-"Searching for matching dates";print(out2, quote=F)

#the following loop isolates the dates indicated in iso, 
#extracts the relevant data from sdata, and writes the results
#to the file OUTDATA

while(j<=length(iso[,ycol])){
 if(iso[j,ycol]==sdata[i,1]){
  if(iso[j,mcol]==sdata[i,2]){
   if(iso[j,dcol]==sdata[i,3]){
    for(k in -3:3){
     if(sdata[i+k,4]==-999){ok<-0;k<-4}
    }
    if(ok==1){
     out<-array(c(sdata[i+1,1], sdata[i+1,2], sdata[i+1,3], sdata[i-3,4]/10, sdata[i-2,4]/10,sdata[i-1,4]/10,sdata[i,4]/10,sdata[i+1,4]/10, sdata[i+2,4]/10, sdata[i+3,4]/10),dim=c(1,10))
#     write(out, file=OUTDAT, ncolumns=10, append=T)
     print(out)
     for(k in -2:1){
      if(sdata[i+k,4]!=0){
       tally<-array(c(sdata[i+k+1,1], sdata[i+k+1,2], sdata[i+k+1,3], k,sdata[i+k,4]/10), dim=c(1,5))
       write(tally, file=OUTTALLY, ncolumns=5, append=T)
      }
     }
     numdata<-numdata+1
    };
   ok<-1;
   j<-j+1
   }
  }
 };
 if(i==length(sdata[,1])){
  j<-length(iso[,ycol])+1;
 }else{i<-i+1}
}

outc<-"Matches were found for ";outd<-" of the requested dates."
print(outc, quote=F);print(numdata);print(outd,quote=F)

#read the output file just created
res<-read.table(OUTDAT, header=FALSE)
Dm3<-mean(res[,4])
Dm2<-mean(res[,5])
Dm1<-mean(res[,6])
D<-mean(res[,7])
Dp1<-mean(res[,8])
Dp2<-mean(res[,9])
Dp3<-mean(res[,10])

assign("Mean", c(mean(res[,4]), mean(res[,5]),mean(res[,6]),mean(res[,7]),mean(res[,8]),mean(res[,9]),mean(res[,10])))

res2<-read.table(OUTTALLY, header=F)


#create boxplots of the data for each composite day, and save the resultant
#graph to the file OUTPLOT, in postcript format

postscript(OUTPLOT)
boxplot(res[,4],res[,5],res[,6],res[,7],res[,8],res[,9],res[,10], main=TITLE, xlab=XLAB, boxwex=0.7, ylab=YLAB)
add=T
points(Mean, type="l", col=4, lty="dashed")
legend(1,15,legend<-c("mean"), col=4,lty="dashed")
dev.off()

plot(res2[,5])
for(h in 1:10){
 m<-1;lowCount<-0
 while(m<=length(res2[,1])){
 if(res2[m,5]<=(0.1*h)){
   lowCount<-lowCount+1
  };
  m<-m+1
 }
 print(0.1*h);print(lowCount/length(res2[,1]))
}