global m n nr co Ir Hr cL
global rho mu grav
global Rey ff
m=21;n=17;nr=3;
co=[1 2; 2 3; 3 4; 3 10; 4 5; 4 7; 5 6; 6 13; 7 6; 7 10;...
    8 9; 9 7; 10 11; 11 12; 11 15; 12 13; 13 14; 14 16;...
    15 17; 16 15; 17 14];
Ir=[1 8 16];Hr=[210 210 200];
cL=-[0 0 0.025 0.06 0.05 0.04 0.03 0 0 0.03...
    0.025 0.03 0.03 0.03 0.03 0 0.02];
horas=[3 6 9 12 19 22 24];
pf=[0.5 0.2 1 1.4 1.1 1.6 1.1];
rho=997;mu=0.890e-3;grav=9.80665;
q=0.01*ones(m,1);h=200*ones(n,1);vR=0*ones(nr,1);
x0=[q;h;vR];

time=0; ## em segundos
timeh=0; ## em horas
dt=1800; ## em segundos
dth=dt/3600; ## em horas
nt=24/dth;
cL0=cL;
k=1;
q=[];h=[];vR=[];tempos=[];hR=[];
Reynolds=[];ffactor=[];
## htanque=Hr(3); ## condicao inicial
for i=1:nt
 if (timeh > horas(k))  k=k+1; endif
 cL=cL0*pf(k); 
 [x fval info]=fsolve(@fhidr2,x0);
 tempos=[tempos timeh];
 q=[q x(1:m)];h=[h x(m+1:m+n)];vR=[vR x(m+n+1:m+n+nr)];
 hR=[hR Hr(3)];
 if (size(Rey,2)>1) Rey=Rey'; endif
 if (size(ff,2)>1) ff=ff'; endif
 Reynolds=[Reynolds Rey];ffactor=[ffactor ff];
 Hr(3)=Hr(3)-dt*x(m+n+3)/(pi*100);
 timeh=timeh+dth; x0=x;
endfor
plot(tempos,hR,"linewidth",2,"-o")