function y=cor26_ks(x)
%version with shortened iteration
%x=[1 0 30. 29.8 24.2 18.5 141 419 0 600 1010 15 15 15];
u=x(1);
us=x(2);
ts=x(3);
t=x(4);
Qs=x(5)/1000;
Q=x(6)/1000;
Rs=x(7);
Rl=x(8);
rain=x(9);
zi=x(10);
P=x(11);
zu=x(12);
zt=x(13);
zq=x(14);

     %***********   set constants *************
     Beta=1.25;
     von=0.4;
     fdg=1.00;
     tdk=273.16;
     grav=9.82;
     %*************  air constants ************
     Rgas=287.1;
     Le=(2.501-.00237*ts)*1e6;
     cpa=1004.67;
     cpv=cpa*(1+0.84*Q);
     rhoa=P*100/(Rgas*(t+tdk)*(1+0.61*Q));
     visa=1.325e-5*(1+6.542e-3*t+8.301e-6*t*t-4.8e-9*t*t*t);
     %************  cool skin constants  *******
     Al=2.1e-5*(ts+3.2)^0.79;
     be=0.026;
     cpw=4000;
     rhow=1022;
     visw=1e-6;
     tcw=0.6;
     bigc=16*grav*cpw*(rhow*visw)^3/(tcw*tcw*rhoa*rhoa);
     wetc=0.622*Le*Qs/(Rgas*(ts+tdk)^2);
     
     %**************  compute aux stuff *******
     Rns=Rs*.945;
     Rnl=0.97*(5.56e-8*(ts+tdk)^4-Rl);
     
     %***************   Begin bulk loop *******
     
     %***************  first guess ************
     du=u-us;
     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=0.011*usr*usr/grav+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-dter)+.61*ta*dq)/ut^2;
	nits=3;
	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-dter)*von*fdg/(log(zt/zot10)-psit_26(zt/L10));
     qsr=-(dq-wetc*dter)*von*fdg/(log(zq/zot10)-psit_26(zq/L10));

     tkt=.001;
	
     
     %disp(usr)
     
     %***************  bulk loop ************
     for i=1:nits;
     
     zet=von*grav*zu/ta*(tsr+.61*ta*qsr)/(usr*usr);
      %disp(usr)
      %disp(zet);
     zo=0.011*usr*usr/grav+0.11*visa/usr;
     rr=zo*usr/visa;
     L=zu/zet;
     zoq=5e-5/rr^.6;
     zot=zoq;
     usr=ut*von/(log(zu/zo)-psiu_26(zu/L));
     tsr=-(dt-dter)*von*fdg/(log(zt/zot)-psit_26(zt/L));
     qsr=-(dq-wetc*dter)*von*fdg/(log(zq/zoq)-psit_26(zq/L));
     Bf=-grav/ta*usr*(tsr+.61*ta*qsr);
     if Bf>0
     ug=Beta*(Bf*zi)^.333;
     else
     ug=.2;
     end;
     ut=sqrt(du*du+ug*ug);
     hsb=-rhoa*cpa*usr*tsr;
     hlb=-rhoa*Le*usr*qsr;
     qout=Rnl+hsb+hlb;
     dels=Rns*(.065+11*tkt-6.6e-5/tkt*(1-exp(-tkt/8.0e-4))); 	% Eq.16 Shortwave
     qcol=qout-dels;
     alq=Al*qcol+be*hlb*cpw/Le;					% Eq. 7 Buoy flux water

     if alq>0;
     		xlamx=6/(1+(bigc*alq/usr^4)^.75)^.333;				% Eq 13 Saunders
            tkt=xlamx*visw/(sqrt(rhoa/rhow)*usr);			%Eq.11 Sub. thk

     else
            xlamx=6.0;
            tkt=min(.01,xlamx*visw/(sqrt(rhoa/rhow)*usr));			%Eq.11 Sub. thk
     end;
      dter=qcol*tkt/tcw;%  Eq.12 Cool skin
      dqer=wetc*dter;
     
     end;%bulk iter loop
     tau=rhoa*usr*usr*du/ut;                %stress
     
     %****************   rain heat flux ********
     
      dwat=2.11e-5*((t+tdk)/tdk)^1.94;%! water vapour diffusivity
      dtmp=(1.+3.309e-3*t-1.44e-6*t*t)*0.02411/(rhoa*cpa); 	%!heat diffusivity
      dqs_dt=Q*Le/(Rgas*(t+tdk)^2);                        	%!Clausius-Clapeyron
      alfac= 1/(1+0.622*(dqs_dt*Le*dwat)/(cpa*dtmp));      	%! wet bulb factor
      RF= rain*alfac*cpw*((ts-t-dter)+(Qs-Q-dqer)*Le/cpa)/3600;
     %****************   Webb et al. correection  ************
     wbar=1.61*hlb/rhoa/Le+(1+1.61*Q)*hsb/rhoa/cpa/ta;
     hl_webb=wbar*Q*Le;
     %**************   compute transfer coeffs relative to du @meas. ht **********
     Cd=tau/rhoa/du^2;
     Ch=-usr*tsr/du/(dt-dter);
     Ce=-usr*qsr/(dq-dqer)/du;
     %************  10-m neutral coeff realtive to ut ********
     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);
     
lam=13.3;%original Soloviev value
A=.625;%as per Hare et al. 2004
%CO2 variables
s=solco2(ts);
sol=s(1);%mole/m^3/atm
alc=s(2);%dimensionless
scwc=s(3);%schmidt #, water
scac=0.9;%schmidt #, air

%Woolf bubble parameterization
f=3.8e-6*u.^3.4;%whitecap fraction
fd=alc*(1+1./(14*alc/scwc^.5)^(1/1.2))^1.2;%Woolf equation
kbb=4900*f/fd/3600/100;%m/s units

%Fairall et al. 1999 parameterization
ha=lam;%neglect air side sublayer buoyancy effects
hw=lam/A/6*xlamx;%includes water side buoyancy effect
vs1=sqrt(rhow/rhoa)*(hw*sqrt(scwc)+(log(.5./tkt)/.4));%water side
vs2=ha*sqrt(scac)+1./sqrt(Cd)-5+.5*log(scac)/.4; %air side
vtc=usr/(vs1+alc*vs2);%non-bubble xfer velocity
vtcT=vtc+kbb;

alc_oz=1e3;%ozone dimensionless solubility
vtc_oz=usr/(vs1+alc_oz*vs2);%non-bubble ozone xfer velocity

    k_co2=1.85*usr/13.3/sqrt(660)*6/xlamx*sqrt(rhoa/rhow);  %C0@ xfer velocity at Sc=660, no bubble component

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 vtcT vtcT*sol vtc_oz*alc_oz  vs2];
%   1  2    3   4  5   6  7   8  9   10  11   12   13 14    15   16 17 18   19    20      21    22     23        24         25
%hsb=sensible heat flux, W/m^2
%hlb=latent heat flux, W/m^2
%tau=stress, N/m^2
%vtcT=total transfer velocity for CO2 of the form
%Flux = vtcT*(Concwater-alpha*Concair) where alpha is dimensionless solubility
%Flux=sol*vtcT*dPCO2, where sol is solubility in X/atm and dP the CO2 partial pressure difference
%23 collumn=flux/dPCO2, mole/m^2/s/atm
%24 collumn is transfer velocity for ozone, prior to ocean chemical reaction
%25 collumn is airside transfer velocity