% modications following Ogura and Takahasi 1971
% 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.19;                                    % parameter for Kc was 2
%beta=2.67;
at=5e-4;                                   % parameter for Ac %OK
%W=2;                                       % rising rate
MAXW=input('Max Vertical velocity?? in m/s ');
%TOPKM=input('Top of W profilem?? in km ');
V=6;                                       % fall speed of raindrop, was 4
betai=.02 % 0.04glaciation                  % 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

% changes for parabolic plot
nstep=6000/dz;                             % total time steps
z=1000:dz:7000;

W(1)=0;
for i=2:nstep+1;
W(i)=-1*(MAXW-.5)*(z(i)/1000-1)*(z(i)/1000-7)/9+.5;
if (W(i)<=0) 
W(i)=0.5;
end
dt(i)=dz/W(i);                                % delta t is now an array
end
% changes for parabolic plot

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;  %OK
 else
  alpha=0;
 end

 Cc(i)=(qv(i-1)-qv(i))/dt(i);
 Ac(i)=dt(i)*alpha*(qc(i-1)-at);

%collection of rain
Kc(i)=dt(i)*beta*qc(i-1)*qr(i-1)^(0.875);
% rain fallout
Fr(i)=-dt(i)*qr(i-1)*V/b;

% collection of graupel
if (T(i-1)>273)
 Kci(i)=0;
else
 Kci(i)=dt(i)*1.4*beta*qc(i-1)*qi(i-1)^(0.9125);
end

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

% ice fallout
 Fi(i)=-dt(i)*wqi(i-1)*Vi/b;



 qc(i)=qc(i-1)+(Cc(i)-Ac(i)-Kc(i)-Kci(i))*dt(i);
 qr(i)=qr(i-1)+(Ac(i)+Kc(i)+Fr(i)-Gl(i))*dt(i);
 qi(i)=qi(i-1)+(Gl(i)+Fi(i)+Kci(i))*dt(i);
 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);
Kci(1)=Kci(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)',4)
axis([0 .003 0 7000])
xlabel('mixing ratio (gm/gm)')
ylabel('height (m)')
%title('Water Content')
title(sprintf('Water Content max w=%.1f Vf=%.1f Vi=%.1f bi=%.3f\n',MAXW,V,Vi,betai))

t=now;
text(-10^-3,-1100,datestr(t))

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


