clear;

%ncload /backup/mirrored_data/zhup/data/ARMDATA/twpsondewnpnC3.b1.20060123.231500.cdf;

%save sounding_dat2.mat alt pres tdry dp wspd deg;

load data/sounding_dat2.mat;

g=9.8;
cp=1004;
xlv=2450000;
rv=461;
rd=287;
t0=273.16;

zlcl=124*(tdry(1)-dp(1));

tp(1)=tdry(1);

rs(1)=g/cp;
for k=1:length(tdry)-1
if alt(k+1)<=zlcl
rs(k+1)=g/cp;
tp(k+1)=tp(k)-rs(k+1)*(alt(k+1)-alt(k));
else
es=6.11*exp(xlv/rv*tp(k)/t0/(t0+tp(k)));
qs=0.622*es/(pres(k)-es);
tv=(t0+tp(k))*(1+0.608*qs);
rs(k+1)=g/cp-g*qs*xlv*(rd*xlv/cp/rv/(t0+tp(k))-1)/tv/rd/(cp+xlv^2*qs/rv/tv/(t0+tp(k)));
tp(k+1)=tp(k)-rs(k+1)*(alt(k+1)-alt(k));
end
end

figure(1);
plot(tdry,alt/1000);
hold on;
plot(tp,alt/1000,'r');
plot(dp,alt/1000,'g');

hold off;

 
figure(2)
plot(rs*100,alt);


