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

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 vapor
qc(1)=0;                                   % initial cloud water
qr(1)=0;                                   % initial rain water
qi(1)=0;                                   % initial ice

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)=dt*alpha*(qc(i-1)-at);
 Kc(i)=dt*beta*qc(i-1)*qr(i-1);
 Fr(i)=-dt*qr(i-1)*V/b;

 if(T(i-1)>273)
  Gl(i)=0;
 else
  Gl(i)=betai*qr(i-1);
 end

 Fi(i)=-qi(i-1)*Vi/b;

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

qt=qv+qc+qr+qi;

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

figure
subplot(1,2,1) 
%plot(qv,z,'b-',qc,z,'g-',qr,z,'r-',qi,z,'m-',qt,z,'c-')
plot(qc,z,'g-',qr,z,'r-',qi,z,'m-')
%hold on
%set(plot(qt,z),'linewidth',[1.5])
%legend('q_v(vapor)','q_c(cloud water)','q_r(rain)','q_i(ice)','q_t(total)')
legend('q_c(cloud water)','q_r(rain)','q_i(ice)')
axis([0 .003 0 10000])
xlabel('mixing ratio')
ylabel('height (m)')
%title('Water Content')
title(sprintf('Water Content w=%.1f Vf=%.1f Vi=%.1f bi=%.3f\n',W,V,Vi,betai))

subplot(1,2,2)
plot(Cc,z,'b-',Ac,z,'g-',Kc,z,'r-',-Fr,z,'m-',Gl,z,'c-',-Fi,z,'k-')
%hold on
%set(plot(Gl,z,':',-Fi,z),'linewidth',[1.5])
legend('Cc(condensation)','Ac(autoconversion)','Kc(collection)',...
       'Fr(rain fallout)','Gl(glaciation)','Fi(ice fallout)')
axis([0 .000015 0 10000])
xlabel('rate (s^-^1)')
ylabel('height (m)')
%title('Processes')
title(sprintf('Processes w=%.1f Vf=%.1f Vi=%.1f bi=%.3f\n',W,V,Vi,betai))
