% Homework 4
clear
TE=273+[26 8 -9.5 -31 -31];             % environment sounding temperature
PE=[950 600 400 250 100];               % environment sounding pressure
m=5;                                    % number of data points

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

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);                      % height at each data point

lambda=5e-5;                            % entrainment rate
wc(1)=5;                                % intial vertical velocity
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
i=2;                                    % index for time steps 
z(1)=ZE(1);                             % initial height for the parcel

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

 for j=m-1:-1:1
  if(z(i)>ZE(j))                        % interpolate to get pressure
   P(i)=PE(j)*exp(-(z(i)-ZE(j))/H(j));  % for current time
  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

 Tc1=Tc(i-1);                           % two guesses for solving new Tc
 Tc2=Tc1-.1;                            % using the secant method
 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);
 f2=f_conv(Tc2,Tc(i-1),wc(i-1),P(i),qvc(i-1),Te,qve,dt,lambda);

 if(abs(f1)>abs(f2))                    % exchange two guesses
  temp=Tc2;                             % keep Tc1, Tc2 will be updated
  Tc2=Tc1;
  Tc1=temp;
  temp=f2;
  f2=f1;
  f1=temp;
 end

 Tc2=Tc1-f1*(Tc1-Tc2)/(f1-f2);          % update Tc2
 err=abs(Tc1-Tc2);                      % difference in Tc1 and Tc2
end

 Tc(i)=Tc2;                             % new Tc 
 qvc(i)=.622*es(Tc(i))/(P(i)-es(Tc(i)));% new qvc
                                        
                                        % get Cc using eqn(7.24) 
 Cc(i-1)=lambda*(qve-qvc(i-1))*wc(i-1)-(qvc(i)-qvc(i-1))/dt; 

 if(qcc(i-1)>at)                        % Kessler's equation starts
  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;             % buoyancy terms
 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);
                                        % get wc using eqn(7.23)
 wc(i)=wc(i-1)+(2/3*B(i-1)-lambda*wc(i-1)^2)*dt;
 
 i=i+1;                                 % update index
end


for j=1:i-1                             % interpolate to get environment
 Te(j)=interp1(ZE,TE,z(j),'linear');    % temperature at each time step
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;

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

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

figure
plot(qvc,z,'-.',qcc,z,qrc,z,':',qtc,z,'--')
title('Mixing ratios')
legend('water vapor','cloud water','rain water','total water')
xlabel('mixing ratios')
ylabel('height (m)')
