#Be careful of the path in the line below

riverflow=scan("po-pontelagoscuro.txt")
duration=seq(1,32850)
duration=duration/90
riverflow1=sort(riverflow,decreasing=T)
plot(duration,riverflow1,type="l")

plot(duration,riverflow1,type="l",main="Flow duration curve for the Po River at Pontelagoscuro (1920-2009)",lwd=4, col="red",xlab="Duration (days)",ylab="River flow (m3/s)")

riverflow2=array(riverflow,dim=c(365,90))
for (i in 1:90)
riverflow2[,i]=sort(riverflow2[,i],decreasing=T)

#matplot(riverflow2,type="l")

avg=rep(0,365)
for (i in 1:365)
avg[i]=mean(riverflow2[i,])

#plot(duration,riverflow1,type="l",main="Flow duration curve for the Po River at Pontelagoscuro (1920-2009)",lwd=4, col="red",xlab="Duration (days)",ylab="River flow (m3/s)")
lines(avg,lwd=4,col="blue")

standev=rep(0,365)
for (i in 1:365)
standev[i]=sd(riverflow2[i,])

matplot(riverflow2,type="l")
lines(avg,lwd=6,col="red")
lines(avg+standev,lwd=6,col="blue")
lines(avg-standev,lwd=6,col="blue")