% Homework 2 (Warm Cloud Water Continuity)
clear
gamma=6e-3;                                % lapse rate 
To=295;                                    % surface temperature
Po=1e3;                                    % surface pressure
b=1000;                                    % initial radius 
beta=2;                                    % parameter for Kc 
at=5e-4;                                   % parameter for Ac
V=4;                                       % fall speed
W=8;                                       % rising rate 
L=2.5e6;                                   % latent heat of vaporation
Rv=461.5;                                  % gas constant for water vapor
R=287;                                     % gass constant for dry air
dz=8;                                      % delta z
dt=dz/W;                                   % delta t 
nstep=9000/dz;                             % total time steps
g=9.8;                                     % gravitional acceleration

z=1000:dz:10000;
T=To-gamma*z;                              % temperature profile 
P=Po*(1-gamma*z/To).^(g/R/gamma);          % pressure profile
es=6.11.*exp(L./Rv.*(1/273-1./T));         % saturation pressure
qv=.622*es./(P-es);                        % mixing ratio of water vapor
qc(1)=0;                                   % initial cloud water 
qr(1)=0;                                   % intial rain water

for i=2:nstep+1;                           % Kessler's equations
 if(qc(i-1)>at)
  alpha=1e-3;
 else
  alpha=0;
 end

 Cc(i)=(qv(i-1)-qv(i))/dt;
 Ac(i)=alpha*(qc(i-1)-at);
 Kc(i)=beta*qc(i-1)*qr(i-1);
 Fr(i)=-qr(i-1)*V/b;

 qc(i)=qc(i-1)+(Cc(i)-Ac(i)-Kc(i))*dt;
 qr(i)=qr(i-1)+(Ac(i)+Kc(i)+Fr(i))*dt;
 b=b*(T(i)*P(i-1)/(T(i-1)*P(i)))^(1/3);
end

qt=qv+qc+qr;

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

subplot(1,2,1) 
plot(qv,z,'b-',qc,z,'g-',qr,z,'r-',qt,z,'m-')
legend('q_v(vapor)','q_c(cloud water)','q_r(rain)','q_t(total)')
axis([0 .02 0 10000])
xlabel('mixing ratio')
ylabel('height (m)')
title('Water Content')

subplot(1,2,2)
plot(Cc,z,'b-',Ac,z,'g-',Kc,z,'r-',-Fr,z,'m-')
legend('Cc(condensation)','Ac(autoconversion)','Kc(collection)','Fr(fallout)')
axis([0 .00003 0 10000])
xlabel('rate (s^-^1)')
ylabel('height (m)')
title('Processes')
