disp('idiss_u');
kon=0.4;
g=9.82;

x=[6.0 0 28.7 26 24.2 18.5 141 419 0 600 1010 10 10 10];
s=num2str(x(1));
T=273.16+x(4);
z=x(12);
y=cool_lkb(x);
L=y(7);
zet=z/L;
dzetn=.4;

%y=[hsb hlb tau zo zot zoq L usr tsr qsr dter dqer tkt RF hl_webb Cd Ch Ce Cdn_10 Chn_10 Cen_10];
%eps=1e-4;
%ct=1e-3;
tsrv=y(9)+0.61*T*y(10);
ct=tsrv*tsrv*z^(-.667)*psi_ft(zet);
usrv=y(8);
eps=(usrv^3)/(kon*z)*psi_eps(zet);

	us_q=(eps*kon*z/psi_eps(zet))^(1/3);
   ts_q=-(ct*z^.667/psi_ft(zet))^0.5;
   
hsc0=z*(ct/5.9)^.75*(8*kon*g/T)^.5;
hs0=-y(8)*tsrv;
us0=y(8);
gt=-2/3/(1-1/(8*zet));
ge=(20/3/(1-20*zet)^(4/3)-1-1/(7-zet)^2)*zet/psi_eps(zet);
bet=1+.5*gt-.667*ge;
spamex=70;
dspm=100;
dspt=100000;
conv=.01;
for j=0:spamex;
   i=j+1;
   n(i)=j/min(dspm,dspt);
   epsp=eps*(1+j/dspm);
   ctp=ct*(1+j/dspt);
	us_ip(i)=(epsp*kon*z/psi_eps(zet*(1-dzetn)))^(1/3);
	ts_ip(i)=-sqrt(ctp*z^.667/psi_ft(zet*(1-dzetn)));
	hsp(i)=-us_ip(i)*ts_ip(i);
   hscp(i)=z*(ctp/5.9)^.75*(8*kon*g/T)^.5;
   usx=(epsp*kon*z/psi_eps(zet))^(1/3);
   tsx=-sqrt(ctp*z^.667/psi_ft(zet));
   contp(i)=1;
   jam=1;
   zetx=zet;
   while jam>conv,
      zetx0=kon*g*z/T*tsx/usx^2;
      usx=(epsp*kon*z/psi_eps(zetx0))^(1/3);
      tsx=-sqrt(ctp*z^.667/psi_ft(zetx0));
      jam=abs((zetx0-zetx)/max(abs(zetx),.001));
      zetx=zetx0;
      contp(i)=contp(i)+1;  
      zetp(i)=zetx;
   end;   
   us_jp(i)=usx;
   ts_jp(i)=tsx;
   hsjp(i)=-usx*tsx;
   
   epsm=eps*(1-j/dspm);
   ctm=ct*(1-j/dspt)  ;
	us_im(i)=(epsm*kon*z/psi_eps(zet*(1+dzetn)))^(1/3);
	ts_im(i)=-sqrt(ctm*z^.667/psi_ft(zet*(1+dzetn)));
	hsm(i)=-us_im(i)*ts_im(i);
   hscm(i)=z*(ctm/5.9)^.75*(8*kon*g/T)^.5;
   usx=(epsm*kon*z/psi_eps(zet))^(1/3);
   tsx=-sqrt(ctm*z^.667/psi_ft(zet));
   contm(i)=1;
   jam=1;
   zetx=zet;
   while jam>conv,
      zetx0=kon*g*z/T*tsx/usx^2;
      usx=(epsm*kon*z/psi_eps(zetx0))^(1/3);
      tsx=-sqrt(ctm*z^.667/psi_ft(zetx0));
      jam=abs((zetx0-zetx)/max(abs(zetx),.001));
      zetx=zetx0;
      contm(i)=contm(i)+1;  
      zetm(i)=zetx;
      if contm(i)>1000
         break;
      end;
   end;   
   us_jm(i)=usx;
   ts_jm(i)=tsx;
   hsjm(i)=-usx*tsx;
	end;
   
   %plot(n(1:spamex),hsp(1:spamex)/hs0,'^',n(1:spamex),hsm(1:spamex)/hs0,'v',n(1:spamex),hsjp(1:spamex)/hs0,n(1:spamex),hsjm(1:spamex)/hs0,'--');
   plot(n(1:spamex),us_ip(1:spamex)/us0,'^',n(1:spamex),us_im(1:spamex)/us0,'v',n(1:spamex),us_jp(1:spamex)/us0,n(1:spamex),us_jm(1:spamex)/us0,'--');
   set(get(gca,'Xlabel'),'String','\delta\epsilon/\epsilon');
  % set(get(gca,'Ylabel'),'String','H_s/H_{s0}');
  set(get(gca,'Ylabel'),'String','u_*/u_{*0}');
   title(strcat('u = ',s,' m/s'));
   axis([0 spamex/dspm 0 2]);