function y=cor_ice(x)
%version with shortened iteration
%modified for ice or water
%the 5th variable is a switch, x(5)<3 means ice, >3 means water
%Qs is calculated using Buck's equs for ice and water
%uses LKB Rt and Rq for seawater and Andreas 1987 for ice
%presently has fixed zo=3e-4 for ice and Smith for water
%x=[4.5 0 -10 -5 1 2 203 250 0 600 1010 15 15 15];
u=x(1);%wind speed
us=x(2);%surface current
ts=x(3);%bulk water/ice surface temperature
t=x(4);%air temperature
ice=x(5);%ice indicator; 1=snow covered; 
Q=x(6)/1000;%air moisture mixing ratio
Rs=x(7);%downward solar flux
Rl=x(8);%downward longwave flux
rain=x(9);%rainrate (mm/hr)
zi=x(10);%inversion height
P=x(11);%surface pressurre (mb)
zu=x(12);%height of anemometer
zt=x(13);%hieght of thermometer
zq=x(14);%height of hygrometer

%***********   set constants *************
Beta=1.25;%gustiness coeff
von=0.4;%von Karman constant
fdg=1.00;%ratio of thermal to wind von Karman
tdk=273.16;
grav=9.82;%accel of earth grav
%*************  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));							%density of air
visa=1.325e-5*(1+6.542e-3*t+8.301e-6*t*t-4.8e-9*t*t*t);	%kinematic viscosity
     %************  cool skin constants  *******
Al=2.1e-5*(ts+3.2)^0.79;			%water thermal expans coeff
be=0.026;								%water salinity expans coef
cpw=4000;								%water specific heat			
rhow=1022;								%water density
visw=1e-6;								%water kinematic visconsity
tcw=0.6;									%water thermal diffusivity
bigc=16*grav*cpw*(rhow*visw)^3/(tcw*tcw*rhoa*rhoa);
     
%**************  compute aux stuff *******
Rns=Rs*.945;							%net solar flux
Rnl=0.97*(5.56e-8*(ts+tdk)^4-Rl);%net longwave flux
     
     %***************   Begin bulk loop *******
     
     %***************  first guess ************
if ice<4     
	if ts<=0 
   	Qs=qice([ts,P])/1000;
	else
   	Qs=qwat([ts,P])/1000;
	end;
else
   Qs=qsee([ts,P])/1000;
end;
wetc=0.622*Le*Qs/(Rgas*(ts+tdk)^2);


        
        
du=u-us;
dt=ts-t-.0098*zt;
dq=Qs-Q;
ta=t+tdk;
ug=.5;
if ice<3
   dter=0;
else
   dter=0.3; 
end;

ut=sqrt(du*du+ug*ug);
zogs=1e-4;
if ice<3
   zogs=3e-4;%snow covered ice roughness guess, Cdn_10=0.00147
end;

     %***************  neutral coeff  **********
u10=ut*log(10/zogs)/log(zu/zogs);
cdhg=von/log(10/zogs);     
usr=cdhg*u10;
if ice<4
   zo10=zogs;
else
   zo10=0.011*usr*usr/grav+0.11*visa/usr;
end;

Cd10=(von/log(10/zo10))^2;
if ice<3
   Ch10=0.0015;
else
   Ch10=0.00115;
end;

Ct10=Ch10/sqrt(Cd10);
zot10=10/exp(von/Ct10);
Cd=(von/log(zu/zo10))^2;
     %***************  Grachev and Fairall (JAM, 1997) *********
Ct=von/log(zt/zot10);		%temp transfer coeff
CC=von*Ct/Cd;					%z/L vs Rib linear coefficient
Ribcu=-zu/zi/.004/Beta^3;	%Saturation Rib
Ribu=-grav*zu/ta*((dt-dter)+.61*ta*dq)/ut^2;
nits=3;							%Number of iterations
if Ribu<0;
	zetu=CC*Ribu/(1+Ribu/Ribcu);		%Unstable, G&F
else;
	zetu=CC*Ribu*(1+27/9*Ribu/CC);	%Stable, I forget where I got this
end;		
L10=zu/zetu;								%MO length
if zetu>50;
	nits=1;									%cutoff iteration if too stable
end;
   %**************  first guess stability dependent scaling params  ****
     usr=ut*von/(log(zu/zo10)-psiu(zu/L10));
     tsr=-(dt-dter)*von*fdg/(log(zt/zot10)-psit(zt/L10));
     qsr=-(dq-wetc*dter)*von*fdg/(log(zq/zot10)-psit(zq/L10));
   	tkt=.001;%cool skin thickness, first guess
		    
     %disp(usr)
     
     %***************  bulk loop ************
for i=1:nits;
 	zet=von*grav*zu/ta*(tsr+.61*ta*qsr)/(usr*usr);
      %disp(usr)
      %disp(zet)
   if ice<3
      zo=zogs;
      %zo=0.018*usr*usr/grav+0.135*visa/usr+1.13e-3*exp(-((usr-0.25)/0.052)^2)/20;
      %zo=0.135*(visa/usr)+1.13e-3*exp(-((usr - 0.25)/0.052)^2)+0.065*(usr^2/grav);
   else
   	zo=0.011*usr*usr/grav+0.11*visa/usr;
   end;
         
   rr=zo*usr/visa;
   if ice<3				%Andreas 1987 for snow/ice
      if rr<=.135,
   		rt=rr*exp(1.250);
     		rq=rr*exp(1.610);
   	elseif rr<=2.5,
     		rt=rr*exp(0.149-.55*log(rr));
     		rq=rr*exp(0.351-0.628*log(rr));
   	elseif rr<=1000,
     		rt=rr*exp(0.317-0.565*log(rr)-0.183*log(rr)*log(rr));
      	rq=rr*exp(0.396-0.512*log(rr)-0.180*log(rr)*log(rr));
   	end;
   
	else				%Liu et al. for ocean
        
   	if rr<=.11,
   		rt=.177;
     		rq=.292;
   	elseif rr<=.8,
     		rt=1.376*rr^.929;
     		rq=1.808*rr^.826;
   	elseif rr<=3.0,
     		rt=1.026*rr^(-.599);
     		rq=1.393*rr^(-.528);
   	elseif rr<=10.0,
     		rt=1.625*rr^(-1.018);
     		rq=1.956*rr^(-.870);
   	elseif rr<=30.0,
     		rt=4.661*rr^(-1.475);
     		rq=4.994*rr^(-1.297);
   	elseif rr<=100.0,
     		rt=34.904*rr^(-2.067);
     		rq=30.709*rr^(-1.845);
   	elseif rr<=300.0,
     		rt=1667.19*rr^(-2.907);
     		rq=1448.68*rr^(-2.682);
   	elseif rr<=1000.0,
     		rt=5.88e5*rr^(-3.935);
     		rq=2.98e5*rr^(-3.616);
   	end;
   end;  
   L=zu/zet;
   zot=rt*visa/usr;
   zoq=rq*visa/usr;
   usr=ut*von/(log(zu/zo)-psiug(zu/L));
   tsr=-(dt-dter)*von*fdg/(log(zt/zot)-psitg(zt/L));
   qsr=-(dq-wetc*dter)*von*fdg/(log(zq/zoq)-psitg(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*(.137+11*tkt-6.6e-5/tkt*(1-exp(-tkt/8.0e-4))); 	% Eq.16 Shortwave
   qcol=qout-dels;
   if qcol>0;
   	alq=Al*qcol+be*hlb*cpw/Le;					% Eq. 7 Buoy flux water
    	xlamx=6/(1+(bigc*alq/usr^4)^.75)^.333;			% Eq 13 Saunders
    	tkt=xlamx*visw/(sqrt(rhoa/rhow)*usr);			%Eq.11 Sub. thk
     	dter=qcol*tkt/tcw;%  Eq.12 Cool skin
   else
   	dter=0;
   end;
   if ice<3
      dter=0;
   end;
   
   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);
     
     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];