% Homework 6, sounding 1 & 2 clear sounding=input('sounding 1 or 2'); % choose sounding if(sounding==1) TE=273+[26 8 -9.5 -31 -31]; % environment sounding temperature PE=[950 600 400 250 100]; % environment sounding pressure m=5; % data point else TE=273+[26 -31 -31]; % environment sounding temperature PE=[950 250 100]; % environment sounding pressure m=3; % data point end 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 % ice or no ice ice=input('ice included?? 1 = yes, 0 = no '); if(ice==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); % 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 qic(1)=0; % initial ice 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)); % at current position 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); % two guesses for the secant method 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)) % exchange Tc1 and Tc2 temp=Tc2; % if |f1| > |f2| Tc2=Tc1; % we'll update Tc2 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 by 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) % parameter for autoconversion 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; % Kessler's equations 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); % update wc by 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 Te Te(j)=interp1(ZE,TE,z(j),'linear'); 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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Eqn (7.22) including ice function y=f_conv(Tc_new,Tc_old,wc,P,qvc_old,Te,qve,dt,lambda,qic_new,qic_old) L=2.5e6; Lf=3e5; Cp=1004; g=9.8; qvc_new=.622*es(Tc_new)/(P-es(Tc_new)); y=Tc_new-Tc_old+g/Cp*wc*dt+L/Cp*(qvc_new-qvc_old)-... lambda*(Te-Tc_old+L/Cp*(qve-qvc_old))*wc*dt-... Lf/Cp*(qic_new-qic_old); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function es(T) function y=es(T) a=17.27; b=35.86; y=3.8/0.62197*exp(a*((T-273)/(T-b)));