% Homework 4
clear
TE=273+[26 8 -9.5 -31 -31];         % environment sounding
PE=[950 600 400 250 100];          

Po=1e3;                             % surface pressure
g=9.8;
Cp=1004;
L=2.5e6;
Rv=461.5;
RH=.8;
V=4;
beta=2;
at=5e-4;
gamma=6e-3;
R=287;
dt=input('dt = ');

deltaz(1)=R/g*TE(1)*log(1e3/PE(1));
for i=1:4
 H(i)=R/g*(TE(i)+TE(i+1))/2;
 deltaz(i+1)=H(i)*log(PE(i)/PE(i+1));
end
ZE=cumsum(deltaz);

lambda=2e-5;
wc(1)=5;
b=.2/lambda;
Tc(1)=TE(1);
P(1)=PE(1);
%es=6.11*exp(L/Rv*(1/273-1/Tc(1)));
qvc(1)=.622*es(Tc(1))/(P(1)-es(Tc(1)));
qcc(1)=0;
qrc(1)=0;
B1(1)=0;
B2(1)=0;
B3(1)=0;
i=2;
h=ZE(1);
z(1)=h;

while wc(i-1)>=0;
 h=h+wc(i-1)*dt;
 z(i)=h;

 if(z(i)>ZE(4))
  P(i)=PE(4)*exp(-(z(i)-ZE(4))/H(4));
 elseif(z(i)>ZE(3))
  P(i)=PE(3)*exp(-(z(i)-ZE(3))/H(3));
 elseif(z(i)>ZE(2))
  P(i)=PE(2)*exp(-(z(i)-ZE(2))/H(2));
 else
  P(i)=PE(1)*exp(-(z(i)-ZE(1))/H(1));
 end

 Te=interp1(ZE,TE,z(i-1),'linear');
 qve=RH*.622*es(Te)/(P(i-1)-es(Te));

 Tc1=Tc(i-1);
 Tc2=Tc1-.1;
 err=1;

while err>1e-10    % solve for Tc(i) using the secant method
 f1=ff(Tc1,Tc(i-1),wc(i-1),P(i),qvc(i-1),dt);
 f2=ff(Tc2,Tc(i-1),wc(i-1),P(i),qvc(i-1),dt);

 if(abs(f1)>abs(f2))
  temp=Tc2;
  Tc2=Tc1;
  Tc1=temp;
  temp=f2;
  f2=f1;
  f1=temp;
 end

 Tc2=Tc1-f1*(Tc1-Tc2)/(f1-f2);
 err=abs(Tc1-Tc2);
end

 Tc(i)=Tc2;
 qvc(i)=.622*es(Tc(i))/(P(i)-es(Tc(i)));

 Cc(i-1)=(qvc(i-1)-qvc(i))/dt;

% entrainment update
 Tc(i)=Tc(i-1)+lambda*(Te-Tc(i-1)+L/Cp*(qve-qvc(i-1)))*wc(i-1)*dt;
 qvc(i)=qvc(i-1)+lambda*(qve-qvc(i-1))*wc(i-1)*dt;

 if(qcc(i-1)>at)
  alpha=1e-3;
 else
  alpha=0;
 end

 Ac(i-1)=alpha*(qcc(i-1)-at);
 Kc(i-1)=beta*qcc(i-1)*qrc(i-1);
 Fr(i-1)=-qrc(i-1)*V/b;

 qcc(i)=qcc(i-1)+(Cc(i-1)-Ac(i-1)-Kc(i-1))*dt+lambda*(-qcc(i-1))*wc(i-1)*dt;
 qrc(i)=qrc(i-1)+(Ac(i-1)+Kc(i-1)+Fr(i-1))*dt+lambda*(-qrc(i-1))*wc(i-1)*dt; 

 B1(i-1)=g*(Tc(i-1)-Te)/Te;
 B2(i-1)=g*.61*(qvc(i-1)-qve);
 B3(i-1)=-g*(qcc(i-1)+qrc(i-1));
 B(i-1)=B1(i-1)+B2(i-1)+B3(i-1);

 wc(i)=wc(i-1)+(2/3*B(i-1)-lambda*wc(i-1)^2)*dt;

 b=b*(P(i-1)*Tc(i)/P(i)/Tc(i-1))^(1/3);
 lambda=.2/b;
 i=i+1;
end

Cc(i-1)=Cc(i-2);
Ac(i-1)=Ac(i-2);
Kc(i-1)=Kc(i-2);
Fr(i-1)=Fr(i-2);

B1(i-1)=B1(i-2);
B2(i-1)=B2(i-2);
B3(i-1)=B3(i-2);
B(i-1)=B(i-2);

qtc=qvc+qcc+qrc;

