floodfreq=function(filed,monthmean=8,meanpref=c(0,0.5,1,1.5))
{
#Be careful: the time series must include complete calendar years of data (start date maust be Jan 01 - End date must be Dec 31st)
#Create the vector of dates
startdate=scan(file=filed, what="character", n=1, skip=4)
enddate=scan(file=filed, what="character", n=1, skip=6)
delbis=scan(file=filed, what="character", n=1, skip=8)
numonth=scan(file=filed, n=1, skip=10)
floodseas=scan(file=filed, n=numonth, skip=12)
flows=scan(file=filed, skip=14)

pr1=seq(as.Date(startdate),as.Date(enddate),1)
if(delbis==F)
{
#Remove 29s Februaries if delbis=T
pr3=flows[strftime(pr1,"%m-%d")!="02-29"]
pr1=pr1[strftime(pr1,"%m-%d")!="02-29"]
}
else
{
pr1=pr1[strftime(pr1,"%m-%d")!="02-29"]
}
#Compute number of years of the time series
iiii=length(pr1[strftime(pr1,"%m-%d")=="01-01"])
#Start the loop to extract the annual maxima in the flood season and the mean flow in each year in the pre-flood season
pr10=0
pr11=0
#Computation of startyear
startyear=strftime(pr1[1],"%Y")
#Extraction of the annual maxima in the flood season
for (i in as.numeric(startyear):(as.numeric(startyear)+iiii-1))
{
pr10=c(pr10,max(flows[as.numeric(strftime(pr1,"%Y"))==i & (as.numeric(strftime(pr1,"%m")) %in% floodseas)==T]))
}
pr10=pr10[2:length(pr10)]
#Extraction of the mean flow in the pre-flood season
for (i in as.numeric(startyear):(as.numeric(startyear)+iiii-1))
{
pr11=c(pr11,mean(flows[as.numeric(strftime(pr1,"%Y"))==i & (as.numeric(strftime(pr1,"%m")) %in% monthmean)==T]))
}
pr11=pr11[2:length(pr11)]
#Transformation of the two vectors to the normal space
pr20=qqnorm(pr10,plot=F)$x
pr21=qqnorm(pr11,plot=F)$x
#Computation of the correlation between annual maxima and mean flow
corcoef=cor(pr20,pr21)

#Computation conditional mean and variance
meancond=corcoef*meanpref
stdcond=(1-corcoef^2)^0.5

#Plot annual maxima distribution
alfa=1.283/sd(pr10)
u=mean(pr10)-0.45*sd(pr10)
quant=seq(min(flows)/2, max(flows)*1.5,length.out=1000)
proba=exp(-exp(-alfa*(quant-u)))

#Computation return period
trit=1/(1-proba)
#plot(quant,trit,log="y",type="l",lwd=2,col="red")

pdf("output.pdf")
plot(quant,-log(-log(proba)),type="l",lwd=2,col="red",axes=F,ylim=c(-log(-log(0.2)),-log(-log(0.995))),ylab="Return Period (years)", xlab=expression('Peak river flow (m'^3*'/s)'))
axis(1)
axis(2,at=-log(-log(c(0.2,0.5,0.8,0.9,0.95,0.98,0.99,0.995))),labels=c("","2","5","10","20","50","100","200"))
segments(y0=c(-log(-log(c(0.2,0.5,0.8,0.9,0.95,0.98,0.99,0.995)))),x0=c(0,0,0,0,0,0,0,0),y1=c(-log(-log(c(0.2,0.5,0.8,0.9,0.95,0.98,0.99,0.995)))),x1=c(15000,15000,15000,15000,15000,15000,15000,15000),lty="dotted")
segments(y0=rep(-log(-log(0.2)),15),x0=seq(0,15000,1000),y1=rep(-log(-log(0.995)),15),x1=seq(0,15000,1000),lty="dotted")

#Generation conditional distributions
qp=0
for (i in 1:length(meanpref))
{
pr100=rnorm(10000,meancond[i],stdcond)
pr200=approx(pr20,pr10,xout=pr100)$y
pr200=pr200[is.na(pr200)==F]
#Estimation of the Gumbel distribution parameters
alfa=1.283/sd(pr200)
u=mean(pr200)-0.45*sd(pr200)
proba=exp(-exp(-alfa*(quant-u)))
quant1=u-1/alfa*log(-log(0.995))
qp[i]=quant1
lines(quant,-log(-log(proba)),col="black")
#lines(quant,trit,log="y",type="l",lwd=2,col="black")
}
dev.off()
#return(corcoef)
#return(qp)
}