% Homework 4
clear
TE=273+[26 -31 -31];                    % environment sounding temperature
PE=[950 250 100];                       % environment sounding pressure
m=3;                                    % number of sounding points
wc(1)=5;                                % intial vertical velocity

Po=1e3;                                 % surface pressure
g=9.8;                                  % gravitational acceleration
Cp=1004;                                % specific heat for dry air 
L=2.5e6;                                % latent heat for vaporation
Rv=461.5;                               % gas constant for vapor
RH=.8;                                  % relative humidity for environment
V=4;                                    % fall speed for rain water
beta=2;                                 % parameter for Fr
at=5e-4;                                % parameter for Ac
R=287;                                  % gas constant for dry air
dt=input('dt = ');                      % delta t
index=input('ice included?? 1 = yes, 0 = no  ');
if(index == 1)
 betai=.003;                            % parameter for Gl
 Vi=1.5;                                % ice fall speed
else
 betai=0;
 Vi=0;
end

deltaz(1)=R/g*TE(1)*log(1e3/PE(1));     % calculate height for sounding
for i=1:m-1
 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=5e-5;                            % entrainment rate
b=.2/lambda;                            % initial radius
Tc(1)=TE(1);                            % initial cloud temperature
P(1)=PE(1);                             % initial cloud pressure
qvc(1)=.622*es(Tc(1))/(P(1)-es(Tc(1))); % initial vapor mixing ratio
qcc(1)=0;                               % initial cloud water mixing ratio
qrc(1)=0;                               % initial rain water mixing ratio
qic(1)=0;
i=2;
z(1)=ZE(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%               INTEGRATION STARTS                              %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while wc(i-1)>=0;                       %
 z(i)=z(i-1)+wc(i-1)*dt;

 for j=m-1:-1:1
  if(z(i)>ZE(j))                        % interpolating pressure field
   P(i)=PE(j)*exp(-(z(i)-ZE(j))/H(j));
  end
 end
 
 Te=interp1(ZE,TE,z(i-1),'linear');     % environment temperature
 qve=RH*.622*es(Te)/(P(i-1)-es(Te));    % environment vapor mixing ratio

 if(Tc(i-1)>273)                        % check if glaciation occurs
  Gl(i-1)=0;
 else
  Gl(i-1)=betai*qrc(i-1);
 end

 Fi(i-1)=-qic(i-1)*Vi/b;                % ice fallout
 qic(i)=qic(i-1)+(Gl(i-1)+Fi(i-1))*dt+lambda*(-qic(i-1))*wc(i-1)*dt;

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

while err>1e-10                         % the secant method
 f1=f_conv(Tc1,Tc(i-1),wc(i-1),P(i),qvc(i-1),Te,qve,dt,lambda,qic(i),qic(i-1));
 f2=f_conv(Tc2,Tc(i-1),wc(i-1),P(i),qvc(i-1),Te,qve,dt,lambda,qic(i),qic(i-1));

 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;                             % new Tc 
 qvc(i)=.622*es(Tc(i))/(P(i)-es(Tc(i)));

 Cc(i-1)=lambda*(qve-qvc(i-1))*wc(i-1)-(qvc(i)-qvc(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)-Gl(i-1))*dt+...
        lambda*(-qrc(i-1))*wc(i-1)*dt; 

 B1(i-1)=g*(Tc(i-1)-Te)/Te;             % buoyancy terms
 B2(i-1)=g*.61*(qvc(i-1)-qve);
 B3(i-1)=-g*(qcc(i-1)+qrc(i-1)+qic(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+qic;

figure
plot(wc,z)
title('W_c')
xlabel('W_c (m/s)')
ylabel('height (m)')

figure
plot(B1,z,B2,z,'-.',B3,z,'--')
title('Buoyancy')
legend('B_1 = g*(T_c-T_e)/T_e','B_2 = .61*g*dq_v','B_3 = -g*q_H')
xlabel('buoyancy (m s^-^2)')
ylabel('height (m)')
grid

figure
plot(qvc,z,'-.',qcc,z,qrc,z,':',qic,z,'--')
hold on
set(plot(qtc,z),'linewidth',[1.5])
title('Mixing ratios')
legend('water vapor','cloud water','rain water','ice','total water')
ylabel('height (m)')
xlabel('mixing ratios')
