function [res]=find_Q_nat(cote_bas,ks,D,h,Cd0,S,L,pf,C,Q,sigma)

maxfun=1000;
maxiter=1000;
tolfun=1e-8;
tolx=1e-8;
g=9.81;
kappa=0.41;
U0=Q./L./pf;

Fr=U0./(9.81*pf).^0.5;
Frg=Fr/(1-C^0.5);


if ks==0;ks=1e-34;end

%if Cd0==2
 %   fFr=(min(2.5,Frg^(-4/3)));%
%else
    fFr=(min(1/(1-(Frg.^2)/4),Frg.^(-2/3))).^2;
%end

if Frg>1.5
    fFr=(Frg.^(-2/3)).^2;
end



alpha=1-(1.*C).^0.5-1/2*sigma.*C;
%Cd=Cd0.*(1+0.4./(pf./D).^2);
Cd=Cd0.*(1+1./(pf./D).^2);
%Cd=Cd0.*(0.8-2*C).*(1+0.4./(pf./D).^2).*fFr;
R=(1-sigma*C);%.*(1-C.^0.5).^2;


if pf/h>1; %fFr=1;

    choixturb=1;
    htilde=h./D;
    hstar=pf./h;
    Rh=h.*(hstar-1);
    ustar=(g.*S.*Rh).^0.5;

    Cd1=Cd;
    CdCh=Cd1.*C.*htilde;

 Cf=2./(5.1.*log10(h./ks)+6).^2;
    U0b=(2*g.*S.*R./(Cd1.*C.*h./D+Cf.*R).*h).^0.5;

    % U0=(2*g.*S.*D.*R./(Cd1.*C)).^0.5;

    [P]=fminbnd(@(alphai) resolve_alphaR(alphai,CdCh,R,U0b,hstar,h,C,D,Cd1,ustar,choixturb),1e-5*h,h,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx));

    alpha=P(1);

    beta2=h.*CdCh./alpha./R;
    beta=(beta2).^0.5;
    a1=beta*(hstar-1)/(cosh(beta));
    c=1;

    UhU0=(a1*sinh(beta)+c)^0.5;
    Uh=UhU0*U0b;

    dhp=1-1/kappa*alpha./h.*Uh./ustar;
    z0hp=(1-dhp).*exp(-1*(kappa*Uh./ustar));

    qsup=ustar./kappa.*h.*((hstar-dhp).*(log((hstar-dhp)./z0hp) - 1)-((1-dhp).*(log((1-dhp)./z0hp) - 1)));

        %calcul int�grale dans la canop�e----
        dzinf=0.01;
    Zinf=(0:dzinf:1);
    Uinf=U0b    .*(beta.*Rh./h.*sinh(beta*Zinf)./cosh(beta)+1).^0.5;
    Ub=zeros(size(Uinf));
    Ub(1:end-1)=Uinf(2:end);
    qinf=sum((Uinf(1:end-1)+Ub(1:end-1))/2*dzinf.*h);

    qtot=qinf+qsup;

    PI=0;
    delta=1;

    Umoy=qtot./pf;
    res=abs(U0-Umoy);

else

    hstar=pf/D;
    Re=U0.*pf/1e-6;

    Cd=Cd.*fFr;

    if ks==0
        Cf=0.3164/4.*Re.^(-0.25);
    else
        Cf=2/(5.1*log10(pf/ks)+6)^2;
    end

    N= (1.*Cf)./(pf./D.*Cd.*C);

 U0a=(2*g.*S.*D.*(R)./(Cd.*C.*(1+N))).^0.5;



[u]=fminsearch(@(U0i) find_U0_complet(U0i,pf,C,D,sigma,Cd0,Cf,S),U0a,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx));

    res=abs(U0-u);

end