%Cloud Entrainment Model

clear
ET=[299 281 263.5 242 242];               %environment temp from sounding
EP=[950 600 400 250 100];                 %environment press from sounding
m=5;                                      %# of sounding pts.
R=287;                                    %Gas constant for dry air
Rv=461.5;                                 %Gas constant for water vapor
Cp=1004;                                  %Specific Heat for dry air
L=2.5e6;                                  %Latent heat for evaporation
g=9.8;                                    %gravity
RH=.8;                                    %relative humidity of environ.
V=4;                                      %Rain water fall speed
Ps=1e3;                                   %Surface pressure
beta=2;                                   %Collection Parameter
at=5e-4;                                  %Accumulation Parameter
dt=input('dt= ');                         %Time step
 
%calulate the height of the sounding

dz(1)=R/g*ET(1)*log(1e3/EP(1));           %initialize dz
for i=1:m-1
 H(i)=R/g*(ET(i)+ET(i+1))/2;              %get H from ave. of last step
 dz(i+1)=H(i)*log(EP(i)/EP(i+1));         %get dz for each layer
end
EZ=cumsum(dz);                            %accumulate the sum to get z

lambda=5e-5;                              %entrainment rate
wc(1)=5;                                  %beginning vertical velocity
b=.2/lambda;                              %initial radius
P(1)=EP(1);                               %initial cloud pressure
Tc(1)=ET(1);                              %initial cloud temperature
qvc(1)=.622*es(Tc(1))/(P(1)-es(Tc(1)));   %initial vapor mixing ratio
qrc(1)=0;                                 %init. rain water mix. ratio
qcc(1)=0;                                 %init. cloud water mixing ratio

i=2;                                      %beg.index at 2/need bottom step
z(1)=EZ(1);                               %initialize z


%%%%%%
%BEGIN HEART OF PROGRAM
%%%%%%

while wc(i-1)>=0;                         %while we stil have vert. vel.
 z(i)=z(i-1) + wc(i-1)*dt;                %height equals previous height
 %where is the parcel??

  for j=m-1:-1:1;
   if(z(i)>EZ(j))                         %interpolate the pressure field
    P(i)=EP(j)*exp(-(z(i)-EZ(j))/H(j));   %P(i)=press where the parcel is
   end
  end

  Te=interp1(EZ,ET,z(i-1),'linear');     %calc. envrion. temp.
  qve=RH*.622*es(Te)/(P(i-1)-es(Te));     %this is a function to get qvc
                                          %in the envir. as a f(es)
  %We need the next Tc
  Tc1=Tc(i-1);                            %these are our first two guesses
  Tc2=Tc1-.1;
  err=1;                                  %initialize the error
 
  %Use the SECANT method to get the cloud temp at the first timestep
  %This will solve the problem of the fact that the equation for temp
  %is implicit

 while err>1e-10
  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);

 %loop to switch Tc1 and Tc2 if abs(f1)>abs(f2)
  if(abs(f1)>abs(f2))
   temp=Tc2;                              %temp is a temporary variable
   Tc2=Tc1;
   Tc1=temp;
   temp=f2;                               %now switch f2 and f1
   f2=f1;
   f1=temp;
  end
  
  Tc2=Tc1-f1*(Tc1-Tc2)/(f1-f2);           %update the Tc2
  err=abs(Tc1-Tc2);                       %test weth. to do loop again
 end

%Out of the evil loop! Now we have a pretty good Tc2)

 Tc(i)=Tc2;
 qvc(i)=.622*es(Tc(i))/(P(i)-es(Tc(i)));

%Houze, eq.7.24:
 Cc(i-1)=lambda*(qve-qvc(i-1))*wc(i-1)-(qvc(i)-qvc(i-1))/dt;

%Ready for Kessler's Eqns.
 if(qcc(i-1)>at)                          %check alpha
  alpha=1e-3;
 else
  alpha=0;
 end
 
 Ac(i-1)=alpha*(qcc(i-1)-at);             %Processes
 Kc(i-1)=beta*qcc(i-1)*qrc(i-1);
 Fr(i-1)=-qrc(i-1)*V/b;

%Houze, eq.7.25 qcc,qrc from Kessler's eq. and entrainment
 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;

%Calculate the Buoyancy terms
 B1(i-1)=g*(Tc(i-1)-Te)/Te;      %entrainment of env. air at diff. temp
 B2(i-1)=g*.61*(qvc(i-1)-qve);   %entrainment of less humid air
 B3(i-1)=-g*(qcc(i-1)+qrc(i-1)); %dec. in buoyancy from inc. qcc and qrc
 B(i-1)=B1(i-1)+B2(i-1)+B3(i-1);

%Houze eqn.7.23 Update the vertical velocity
%The 2/3 is an estimate of the effect of the drag (Houze, p.235)
 wc(i)=wc(i-1)+(2/3*B(i-1)-lambda*wc(i-1)^2)*dt;
 i=i+1;                          %update index
end

%go back to the first while and start a new time step. Do this as long as
%the parcel is still rising

%Now update processes and buoyancy
Cc(i-1)=Cc(i-2);
Ac(i-1)=Ac(i-2);
Fr(i-1)=Fr(i-2);
Kc(i-1)=Kc(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;                           %total cloud water

%Done with the heart of the program
figure
plot(wc,z)
title('Updraft Velocity, W_c')

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')
grid

figure
plot(qvc,z,'--',qcc,z,qrc,z,':',qtc,z,'-.')
title('Cloud Mixing Ratios')
legend('cloud water vapor','cloud water','rain water','total water')
%%%%%%%%%%%%%%%

