function y=corland26(x)
%version with shortened iteration; modified Rt and Rq
%x=[5 30 30  18.5 180 600 1010 10 2 2];%sample data line
u=x(1);%wind speed (m/s)
ts=x(2);%surface temperature (C)
t=x(3);%air temperature (C)
Q=x(4)/1000;%air humidity (input g/kg  converted to g/g)
dir=x(5);%wind direction for future ref to zo
zi=x(6);%PBL height
P=x(7);%air pressure (mb)
zu=x(8);%wind speed measurement ht
zt=x(9);%air temp measurement ht
zq=x(10);%humidity measurement ht
jd=x(11);%julian day
lat=38;%latitude for nashville
h=1.00;%stomatal resistance factor (1.00 implies no resistance), no longer used
Qs=qgnd([ts P h])/1000;
if isnan(P);P=1000;end;

     %***********   set constants *************
     Beta=1.2;%gustiness factor
     von=0.4;%von karman
     fdg=1.00;%lkb factor
     tdk=273.16;
     grav=grv(lat);%9.82;
     %*************  air constants ************
     Rgas=287.1;%gas constant
     Le=(2.501-.00237*ts)*1e6;%latent heat of vaporization of water
     cpa=1004.67;%specific heat of air, constant p
     cpv=cpa*(1+0.84*Q);%inlcuding moisture
     rhoa=P*100/(Rgas*(t+tdk)*(1+0.61*Q));%air density
     visa=1.326e-5*(1+6.542e-3*t+8.301e-6*t*t-4.84e-9*t*t*t);%kinematic viscosity
     %**************  compute aux stuff *******
     %***************   Begin bulk loop *******
     
     %***************  first guess ************
     du=u;
     dt=ts-t-.0098*zt;
     dq=Qs-Q;
     ta=t+tdk;
     ug=.5;
     dter=0.3; 
     ut=sqrt(du*du+ug*ug);
	  u10=ut*log(10/1e-4)/log(zu/1e-4);
     usr=.035*u10;
	zo10=.01+0.11*visa/usr;%
	Cd10=(von/log(10/zo10))^2;
	Ch10=0.00115;
	Ct10=Ch10/sqrt(Cd10);
	zot10=10/exp(von/Ct10);
	Cd=(von/log(zu/zo10))^2;
	Ct=von/log(zt/zot10);
	CC=von*Ct/Cd;
	Ribcu=-zu/zi/.004/Beta^3;
	Ribu=-grav*zu/ta*((dt)+.61*ta*dq)/ut^2;
	nits=6;
	if Ribu<0;
		zetu=CC*Ribu/(1+Ribu/Ribcu);
	else;
		zetu=CC*Ribu*(1+27/9*Ribu/CC);
		end;		
	L10=zu/zetu;
	if zetu>50;
		nits=1;
	end;
     usr=ut*von/(log(zu/zo10)-psiu_26(zu/L10));
     tsr=-(dt)*von*fdg/(log(zt/zot10)-psit_26(zt/L10));
     qsr=-(dq)*von*fdg/(log(zq/zot10)-psit_26(zq/L10));

     tkt=.001;
	
   charn=0.011;
   if ut>10
      charn=0.011+(ut-10)/(18-10)*(0.018-0.011);
   end;
   if ut>18
      charn=0.018;
   end;
   
     %disp(usr)
     
     %***************  bulk loop ************
     for i=1:nits;
     
     zet=von*grav*zu/ta*(tsr+.61*ta*qsr)/(usr*usr);
      %disp(usr)
      %disp(zet);
      
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %specify surface roughness
      zo=.2+0.11*visa/usr;%area average  roughness
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
      %%%%   allow zo to depend on wind direction, not used
      %if isnan(dir);dir=0;end;
      %zo=interp1([0 45 135 225 315 360],[.2 .12 .36 .4 .89 .2],dir);
      %zo=interp1([0 22 67 112 157 202 247 292 337 360],[.21 .073 .051 .18 .21 .28 .57 .7 .64 .21],dir);
     rr=zo*usr/visa; %roughness reynolds #
     L=zu/zet;
     
     %%%%%%%%%%%   specify temperature roughness length  *********
     bhm=min(6*rr^.25-5,2/.4);%k u* Rb from Garratt
     zot=zo/exp(.4*bhm);
     %kurc=interp1([165 168 170 183 200],[15 20 10 10 10],jd);
     
     %%%%%%%%%%  interpolate time series of specified moisture resistance, kapp u* Rc
   	jdx=165:1:200;
		kurcx=[25 25 25 30 35 25 20 20 20 20 15 9  9  9  9  9  12 15 12  9  8  7  5  8 8  15  9  9 12 10 10 10 10  7  7  7];
       %     65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 00 

      kurc=interp1(jdx,kurcx,jd);
     
     zoq=zot/exp(kurc);%k u* Rc = 10
     %zoq=zot;
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
     
     usr=ut*von/(log(zu/zo)-psiu_26(zu/L));%friction velocity
     tsr=-(dt)*von*fdg/(log(zt/zot)-psit_26(zt/L));%t scaling parameter
     qsr=-(dq)*von*fdg/(log(zq/zoq)-psit_26(zq/L));%q scaling parameter
     Bf=-grav/ta*usr*(tsr+.61*ta*qsr);%buoyancy flux
     if Bf>0
     ug=Beta*(Bf*zi)^.333;%gustiness
     else
     ug=.2;
     end;
     ut=sqrt(du*du+ug*ug);
     hsb=-rhoa*cpa*usr*tsr;%sensible heat flux
     hlb=-rhoa*Le*usr*qsr;%latent heat flux
     
     end;%bulk iter loop
     tau=rhoa*usr*usr*du/ut;%stress
     
     %****************   rain heat flux ********
     
     
     %****************   Webb et al. correection  ************
     wbar=1.61*hlb/Le/(1+1.61*Q)/rhoa+hsb/rhoa/cpa/ta;%formulation in hlb already includes webb
     %wbar=1.61*hlb/Le/rhoa+(1+1.61*Q)*hsb/rhoa/cpa/ta;
     hl_webb=rhoa*wbar*Q*Le;
     %**************   compute transfer coeffs relative to du @meas. ht **********
     Cd=tau/rhoa/du^2;
     Ch=-usr*tsr/du/(dt);
     Ce=-usr*qsr/(dq)/du;
     %************  10-m neutral coeff realtive to 10-m ********
     Cdn_10=von*von/log(10/zo)/log(10/zo);
     Chn_10=von*von*fdg/log(10/zo)/log(10/zot);
     Cen_10=von*von*fdg/log(10/zo)/log(10/zoq);
   
    
     y=[hsb hlb tau zo zot zoq L usr tsr qsr wbar Cd Ch Ce Cdn_10 Chn_10 Cen_10 ];