function [us_ib,ts_ib,qs_ib]=gasid(x,y,C);

disp('gasid.m');

vkon=0.4;
alphau=0.52;
g=9.81;
Cp=1004;
Cle=2500-2.274*x(4);
T=273.16+x(4);
Tv=T*(1+.61*x(6));
rho=1013.25/287/Tv;
z_u=x(12);
z_q=x(12)-1.7;
L=y(7);
zet_u=z_u/L;
zet_q=z_q/L;
Cu2=C(1);Ct2=C(3);Cq2=C(5);

us_ib=sqrt(Cu2*(z_u^(2/3))/psi_fu(zet_u));
ts_ib=sign(zet_u)*sqrt(Ct2*(z_u^(2/3))/psi_ft(zet_u));
qs_ib=-sqrt(Cq2*(z_q^(2/3))/psi_ft(zet_q))/1000;
ts_ib=ts_ib-0.00051*T*qs_ib/rho;

us0=us_ib;
ts0=ts_ib;
qs0=qs_ib;

%us0=y(8);
us0=us_ib;
ts0=y(9);
qs0=y(10);

for iii=1:10,
   zet=vkon*z_u*g/Tv*(ts0+0.00061*T*qs0/rho)/(us0*us0);
   us0=sqrt(Cu2*(z_u^(2/3))/psi_fu(zet));
   %disp(['zet = ',num2str(zet),'  us0 = ',num2str(us0)]);
   %ts0=sign(zet)*sqrt(Ct2*(z^(2/3))/psi_ft(zet));
   %qs0=-sqrt(Cq2*(z^(2/3))/psi_ft(zet))/1000;
   %ts0=ts0-0.00051*T*qs_ib/rho;
end;
disp(['z/L = ',num2str(zet),'     id = ',num2str(us0),'     bulk = ',num2str(y(8))])

us_ib=us0;
