function zd=duct(g)
%von Karman constant
kon=0.4;
%g=[7.6 16 14.5 70 1010 600 250 400 15 15 15]
%Sample bulk meteorological input data
u=g(1);% 			%m/s    %wind speed
Ts=g(2);%    		%C, water surface T  
T=g(3);%				%C, air temperature
H=g(4);%				%%, air relative humidity
P=g(5);%				%mb, atmosperhic pressure
zi=g(6);%			%m, estimate of atmospheric mixed-layer depth
Rs=g(7);%			%W/m^2, estimate of downward solar radiative flux
Rl=g(8);%			%E/^2, estimate of downward longwave radiative flux

%Measurement heights for u, T, and RH
zu=g(9);%			%m, height of wind speed measurement
zT=g(10);%			%m, height of mean temperature measurement
zh=g(11);%			%m, height of mean humidity measurement

%compute humidity in g/kg variables

esats=6.112*exp(17.502*Ts/(241+Ts))*.98;   %sat vapor pressure for seawater at temp=Ts
qs=0.622*esats/(P-0.378*esats)*1000;		%surface humidity, g/kg
es=6.112*exp(17.502*T/(241+T))*H/100;  		%vapor pressure at temp=T and relhum=H
q=0.622*es/(P-0.378*es)*1000;					%air humidity, g/kg


x=[u 0 Ts T qs q Rs Rl 0 zi P zu zT zh];

y=coul_lkb(x);			%compute bulk fluxes and scaling parameters with Toga Coare version of LKB algortihm  (Fairall et al., JGR, vol. 101, 3747-3764, 1996)
                     
%set MO scaling paras
L=y(7);
tsr=y(9);
qsr=y(10)*1000;
dqer=y(12)*1000;
zoq=y(6);
%first guess
zd=5;
deduct=1;
dzet=0.01;
zdoo=8*(7.2*qsr-1.3*tsr)/kon;
zet=zdoo/L;
phi=1-zet*(psit(zet+dzet/2)-psit(zet-dzet/2))/dzet;
zdo=zdoo*phi;
zd=zdo;


%set critical gradient condition
while abs(deduct)>.01,
   zet=zd/L;
   phi=1-zet*(psit(zet+dzet/2)-psit(zet-dzet/2))/dzet;
	deduct=0.157-0.032+(7.2*qsr-1.3*tsr)*phi/zd/kon;

   zd=zd-0.25*deduct;
   if zd<0
      zd=0.01;
      break;
   end;
	%compute q at zd
	qd=qs-dqer+qsr/kon*(log(zd/zoq)-psit(zet));
	if q<0
   	break;
	end;
end;

