disp('bulk_err')

bet=1.25;
k=0.4;
al=9.8/300;

ue=[.2 .03];
te=[.3 0];
qe=[.3 0];

ecu=.002;
ect=.002;
ecq=.002;

z=10;
Ts=29.2;
qs=24.8;

u=5;
dt=1.6;
dq=6.7;
dq=3.0;
x=[4.5 0 29.2 27.6 24.8 18.1 203 413 0 600 1010 10 10 10];
x=[4.5 0 15.2 13.6 10 7 203 353 0 600 1010 10 10 10];


%dt=-.5:-1:-4;
dt=.5:1:4;

[b,m]=size(dt);

u=.5:.5:15;
[b,n]=size(u);
x(12)=z;x(13)=z;x(14)=z;

for j=1:m,
dth=dt(j)+6.1e-4*(Ts-dt(j)+273.16)*dq;

for i=1:n,

zi=x(10);

x(1)=u(i);
x(3)=Ts;
x(4)=Ts-dt(j);
x(5)=qs;
x(6)=qs-dq;

eu=sqrt(ue(1)^2+(ue(2)*u(i))^2);
et=sqrt(te(1)^2+(te(2)*dt(j))^2);
eq=sqrt(qe(1)^2+(qe(2)*dq)^2);

y=lkb(x);

%y=[hsb hlb tau zo zot zoq L usr tsr qsr dter dqer tkt RF hl_webb Cd Ch Ce Cdn_10 
%    1    2  3   4   5  6  7  8   9   10  11   12   13  14  15    16 17 18  19
%Chn_10 Cen_10];
hs(j,i)=y(1);
hl(j,i)=y(2);
tau(j,i)=y(3);
hb(j,i)=hs(j,i)+6.1e-4*300*hl(j,i)/2.5;
zet=z/y(7);
zt(j,i)=zet;
cunh=k/log(z/y(4));
ctnh=k/log(z/y(5));
cqnh=k/log(z/y(6));

dcu=zet*(psiu(zet)-psiu(.99*zet))/.01/zet;
dct=zet*(psit(zet)-psit(.99*zet))/.01/zet;
dcq=dct;

usr=y(8);
ust(j,i)=usr;
tsr=y(9);
qsr=y(10);
ws=(al*usr*abs(tsr+6.1e-4*300*qsr)*zi)^.333;

su=(1-cunh/k*psiu(zet));
st=(1-ctnh/k*psit(zet));
sq=(1-cqnh/k*psit(zet));
uw=sqrt(u(i)*u(i)+(bet*ws)^2);
a=((bet*ws)/uw)^2;

fs=cunh/k*dcu/su+a/3;
fu=cunh/k*dcu/su;
ft=ctnh/k*dct/st;
fq=cqnh/k*dcq/sq;

xs1=ecu/su/ctnh;
xt1=ect/st/ctnh;
xq1=ecq/sq/cqnh;
xu1=ecu/su/cunh;

xs2=eu/u(i)*(u(i)/uw)^2;
xt2=et/dt(j);
xq2=eq/dq;
xu2=eu/u(i);

d(j,i)=(1-ft)*(1-a)+2*fs;

dus(j,i)=sqrt((1-ft)^2*(xs1^2+xs2^2)+fs^2*(xt1^2+xt2^2))/d(j,i);

dzet(j,i)=sqrt((1-a)^2*(xt1^2+xt2^2)+4*(xs1^2+xs2^2))/d(j,i);

dwth(j,i)=sqrt(((1-a)+3*fs)^2*(xt1^2+xt2^2)+(xs1^2+xs2^2)*(1-3*ft)^2)/d(j,i);

dwq(j,i)=sqrt(xq1^2+xq2^2+((fq*(1-a)+fs)^2*(xt1^2+xt2^2)+((1-ft)-2*fq)^2*(xs1^2+xs2^2))/d(j,i).^2);

dwu(j,i)= sqrt((xu1*((1-ft)*(2-a)+.67*a))^2+(xu2*((1-ft)*(1+(u(i)/uw)^2-a)+2*fu*(1-(u(i)/uw)^2)+.67*a))^2+(xt1^2+xt2^2)*(fu*(2-a)+.33*a)^2)/d(j,i);


end;%for i
end;%for j

figure;
semilogy(u,abs(zt));
%axis([0 u(n) 0 1]);
xlabel('Wind speed (m/s)');
ylabel('abs({\zeta}) ');
text(10,9,'{\Delta \Theta}=.5, 1.5, 2.5, 3.5');
title('COARE conditions, unstable');
%hold on;


figure;
plot(u,dzet);
axis([0 u(n) 0 1]);
xlabel('Wind speed (m/s)');
ylabel('{\delta\zeta}/{\zeta} ');
text(10,.9,'{\Delta \Theta}=.5, 1.5, 2.5, 3.5');
title('COARE conditions, unstable');
%hold on;

figure;
plot(u,ust);
axis([0 u(n) 0 .6]);
xlabel('Wind speed (m/s)');
ylabel('u_{*t} (m/s)');
text(5,.55,'{\Delta \Theta}=.5, 1.5, 2.5, 3.5');
title('COARE conditions, unstable');


figure;
plot(u,dus);
axis([0 u(n) 0 .4]);
xlabel('Wind speed (m/s)');
ylabel('{\delta}u_{*t}/u_{*t} ');
text(10,.38,'{\Delta \Theta}=.5, 1.5, 2.5, 3.5');
title('COARE conditions, unstable');
%hold on;

figure;
plot(u,hl);
axis([0 u(n) 0 300]);
xlabel('Wind speed (m/s)');
ylabel('H_l (W/m^2) ');
text(5,275,'{\Delta\Theta}_v=.5, 1.5, 2.5, 3.5');
title('COARE conditions, unstable');
%hold on;


figure;
plot(u,dwq);
axis([0 u(n) 0 .4]);
xlabel('Wind speed (m/s)');
ylabel('{\delta}H_l/H_l ');
text(10,.38,'{\Delta\Theta}=.5, 1.5, 2.5, 3.5');
title('COARE conditions, unstable');
%hold on;

figure;
plot(u,tau);
axis([0 u(n) 0 .4]);
xlabel('Wind speed (m/s)');
ylabel('Tau (N/m^2)');
text(10,.38,'{\Delta\Theta}=.5, 1.5, 2.5, 3.5');
title('COARE conditions, unstable');
%hold on;

figure;
plot(u,dwu);
axis([0 u(n) 0 .4]);
xlabel('Wind speed (m/s)');
ylabel('{\delta}Tau/Tau ');
text(10,.38,'{\Delta\Theta}=.5, 1.5, 2.5, 3.5');
title('COARE conditions, unstable');
%hold on;

figure;
plot(u,hs);
axis([0 u(n) 0 70]);
xlabel('Wind speed (m/s)');
ylabel('H_{s} (W/m^2) ');
text(5,65,'{\Delta\Theta}=.5, 1.5, 2.5, 3.5');
title('COARE conditions, unstable');
%hold on;


figure;
plot(u,dwth);
axis([0 u(n) 0 1]);
xlabel('Wind speed (m/s)');
ylabel('{\delta}H_{sb}/H_{sb} ');
text(10,.90,'{\Delta\Theta}=.5, 1.5, 2.5, 3.5');
title('COARE conditions, unstable');
%hold on;
