hymod=function(par) { areab=1314000000 tstep=3600 ndeltat=6960 #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/(kslow*fatconv) wquick=0 er1=rep(0,ndeltat) er2=rep(0,ndeltat) er=rep(0,ndeltat) e=rep(0,ndeltat) qt=rep(0,ndeltat) # computation loop # for (i in 1: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 # wslow=(1-kslow)*wslow+(1-kslow)*uslow qslow=(kslow/(1-kslow))*wslow # Quick flow # wquick=rep(0,3) qquick=0 for (j in 1:3) { wquick[j]=(1-kquick)*wquick[j]+(1-kquick)*uquick qquick=(kquick/(1-kquick))*wquick[j] uquick=qquick } # Total flow # qt[i]=(qslow+qquick)*fatconv } # Drawing the plot plot(qt,type="l",ylim=c(0,max(c(max(qt),max(qoss)))),col="red",lwd=2) lines(qoss) #cat(par,fill=T) #return(qt) return(sum((qt-qoss)^2)) }