find_Q_nat.m

function [res]=find_Q_nat(cote_bas,ks,D,h,Cd0,S,L,pf,C,Q,sigma)
fprintf('*************************************\n')
fprintf('find_Q_nat Q=%f\n',Q)

maxfun=5000;
maxiter=5000;
tolfun=1e-16;
tolx=1e-16;
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


coeff_contraction=0.4*Cd0+0.7;

%if Cd0==2
 %   fFr=(min(2.5,Frg^(-4/3)));%
%else
    fFr=(min(coeff_contraction./(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+1./(pf./D).^2).*fFr;

%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.1; %fFr=1;
    
    choixturb=1;
    htilde=h./D;
    hstar=pf./h;
    Rh=h.*(hstar-1);
    ustar=(g.*S.*Rh).^0.5;
    fprintf('ustar=%f\n',ustar)
    
    Cd1=Cd0.*(1+0.4./(pf./D).^2).*fFr;
    CdCh=Cd1.*C.*htilde;
    fprintf('CdCh=%f\n',CdCh)
    
    Cf=2./(5.1.*log10(h./ks)+6).^2;
    
    U0b=(2*g.*S.*R./(Cd1.*C.*h./D+Cf.*R).*h).^0.5;
    fprintf('U0b=%f\n',U0b)
    
    % U0=(2*g.*S.*D.*R./(Cd1.*C)).^0.5;
    
    [P]=fminbnd(@(alphai) resolve_alpha(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);
    fprintf('alpha=%f\n',alpha)
    
    beta2=h.*CdCh./alpha./R;
    beta=(beta2).^0.5;
    a1=beta*(hstar-1)/(cosh(beta));
    c=1;