#Uploading rainfall, evapotranspiration and river flow data
rainevapo=scan("rain-evaposynt.txt",skip=1)
rainevapo=array(rainevapo,dim=c(2,8760))
p=rainevapo[1,]
ep=rainevapo[2,]
qoss=scan("dischargesynt.txt",skip=1)

#Function to apply hymod

hymod=function(par)
{
areab=1214000000
tstep=3600
ndeltat=8760 #be careful to adjust it as needed
qi=15
fatconv=1/1000/tstep*areab

# parameters #

cmax=par[1]
beta=par[2]
alfa=par[3]
kslow=par[4]
kquick=par[5]

# initialisation of variables #

w2=0
c1=0
wslow=qi/(1/kslow*fatconv)
wquick=rep(0,3)
er1=rep(0,ndeltat)
er2=rep(0,ndeltat)
er=rep(0,ndeltat)
e=rep(0,ndeltat)
qt=rep(0,ndeltat)

# computation loop #
sk=1 #Set initial time step of simulation

for (i in sk:ndeltat)

  {
  
  w1=w2
  c1=cmax*(1-((1-((beta+1)*w1/cmax))^(1/(beta+1))))   
  c2=min(c1+p[i],cmax)
  er1[i]=max((p[i]-cmax+c1),0)
  w2=(cmax/(beta+1))*(1-((1-c2/cmax)^(beta+1)))  
  er2[i]=max((c2-c1)-(w2-w1),0)
  e[i]=(1-(((cmax-c2)/(beta+1))/(cmax/(beta+1))))*ep[i]
  w2=max(w2-e[i],0)
  er[i]=er1[i]+er2[i]
  
# Subdivision of the surface runoff #  
  
  uquick=alfa*er[i]
  uslow=(1-alfa)*er[i]
  
# Slow flow #

# The linear reservoir is resolved with the following relationship:
# w(t)=w(t-1)+p(t) - (w(t-1)+p(t))/k*Dt = w(t-1)+p(t) - (w(t-1)+p(t))/(k/Dt)
# If k is measured in units of Dt then Dt=1 and one obtains: 
# w(t) = (w(t-1)+p(t)) - (w(t-1)+p(t))/k = (w(t-i)+p(t))(1-1/k), then:
# w(t) = (1-1/k) w(t-1) + (1-1/k)p(t), which is the equation below.

  wslow=(1-1/kslow)*wslow+(1-1/kslow)*uslow

# Then, from the above proof it follows that q(t) = (w(t-1)+p(t))/k
# Namely: 1/k/(1-1/k) w(t), which is the equation below.
  
  qslow=(1/kslow/(1-1/kslow))*wslow
  
# Quick flow # We repeat the above computation three times.

  qquick=0  
  for (j in 1:3)
    {
      wquick[j]=(1-1/kquick)*wquick[j]+(1-1/kquick)*uquick
      qquick=(1/kquick/(1-1/kquick))*wquick[j]
      uquick=qquick    }

# Total flow #  
  
  qt[i]=(qslow+qquick)*fatconv
  
  }

# Drawing the plot
i=sk
if (max(qt)<1000000){
plot(qt[i:ndeltat],type="l",ylim=c(0,max(c(max(qt[i:ndeltat]),max(qoss[i:ndeltat])))),col="red",lwd=2)
lines(qoss[i:ndeltat])}
#cat(par,fill=T)
#return(qt)
return(sum((qt[i:ndeltat]-qoss[i:ndeltat])^2))
}

#Calibration of parameters on the whole time series
#Initial parameter values are already quite good. Pay attention to factr: the 
#higher its value the faster the convergence. Default would be 1e8

#To make it faster one can remove the instruction for plotting the graph at 
#each iteration

pr1=optim(c(400,0.5,0.8,50,10),hymod,method="L-BFGS-B",lower=c(10,0.1,0,1,1),upper=c(400,10,0.9,1000,1000),control=list(factr=1e10))