ad4.m

fla=1; %1=AED 2=old Bauru style

flag=0;
dp=-10; %pressure step

clear p T th thcs r theta

A=2.53e11;
epsilon=0.622;
B=5.42e3;
k=0.286;
L=2.5e6;
cp=1005;


%read in sounding
if fla==2
    FileName='c:/documents and settings/user/my documents/cprogs/0302052040UTC.TXT';
else
    FileName='c:/documents and settings/user/my documents/HIBISCUS/baurufield/soundings/ascii/1/04_02_06_12AED.txt';
end

fid=fopen(FileName,'rb');

if fla==2
    txt=fscanf(fid,'%s',71);
    dat=fscanf(fid,'%g',[12,108]);
    Tdat=dat(6,:);
    Pdat=dat(5,:);
    rhstart=dat(7,1); %relative humidity at ground (%)
else
    dat=fscanf(fid,'%g',[9,437]);
    Tdat=dat(5,:);
    Pdat=dat(4,:);
    rhstart=dat(6,1); %relative humidity at ground (%)
end



THdat=(Tdat+273.13).*(1e3./Pdat).^k;
pstart=Pdat(1)*100;
tstart=Tdat(1)+273.13;
pend=2000;


% pstart=100000;
% tstart=20+273;

T(1)=tstart;
p(1)=pstart;
thstart=T(1)*(1e5/p(1))^k;
iend=-fix( (pstart-pend)/dp )
esg=A*exp(-B/T(1)); %es at ground
wg=rhstart/100*esg*epsilon/p(1); %mixing ratio at ground=constant

psc=0;

for i=1:iend
dwsdp=-epsilon*A*exp(-B/T(i))/p(i)/p(i);
dwsdt=B*epsilon*A*exp(-B/T(i))/T(i)/T(i)/p(i);
es=A*exp(-B/T(i)); %saturation vapour pressure at T
e=p(i)*wg/epsilon;  %vapour pressure at T

if e>es    %if saturation reached
    psc=psc+1;
    if psc==1
        psati=i;
    end
    dTdp=( (k/p(i)) - (L*dwsdp/T(i)/cp) ) / ( (1/T(i)) + (L*dwsdt/T(i)/cp) ); %psuedo adiabat
    T(i+1)=T(i)+dTdp*dp;
else
    T(i+1)=thstart*(p(i)/1e5)^k;  %dry adiabat
    i
end

p(i+1)=p(i)+dp;

end

th=T.*(1e5./p).^k;

hold on;

figure;
pcs=[1000 800 700 600 500 400 300 200 100 50]; %mb
ths=273:20:413;
thsfine=260:0.1:600;
ts=30:-20:-90;

minT=min(T);
maxT=max(T);
minTH=log10(thsfine(1));
maxTH=log10(thsfine(end));
Y=1/1.2*(maxT-minT); %ratio of page height to width for page transform


size(pcs)
for i=1:ans(2)
    
thcs(i,:)=T.*(1e3./pcs(i)).^k;

hold on;
size(thcs);
siz=ans(2);

y1d=Y*(log10(thcs(1,1))-minTH)/(maxTH-minTH); %in page co-ordinates
y2d=Y*(log10(thcs(1,siz))-minTH)/(maxTH-minTH);

thi=atan( (y2d-y1d) / ( T(siz)-T(1) ) );
%thi=0;


[xx,yy]=rotatedan2(T,log(th)/log(10),thi,Y,minT,minTH,maxTH); 
h=plot(xx,(yy),'r'); %pseudo-adiabat
set(h,'linewidth',2);

[x,y]=rotatedan2(T,log(thcs(i,:))/log(10),thi,Y,minT,minTH,maxTH);
h=plot( x,(y) ); %plot isobars

tol=0;
ii=[];
while size(ii,2)<1
    tol=tol+2
    ii=find(abs(x-xx(3))<tol); %find x value equal to xx(3)
end
h=text(xx(3),y(ii(end)),num2str(pcs(i)));


end




clear x y
size(ths);
for i=1:ans(2)
   [x,y]=rotatedan2(T,log(ths(i))/log(10),thi,Y,minT,minTH,maxTH);
   h=plot( x,(y) );
   h=text(x(fix(siz*15/20)),y(fix(siz*15/20)),num2str(ths(i)),'rotation',-19);
end



clear x y
size(thsfine);
siz=ans(2);
size(ts);
for i=1:ans(2)
   [x,y]=rotatedan2(273.13+ts(i),log(thsfine)/log(10),thi,Y,minT,minTH,maxTH);
   h=plot( x,(y) ); %plot isotherms
   h=text(x(fix(siz*6/20)),y(fix(siz*6/20)),num2str(ts(i)));

end


axis( [xx(end)/1.1 xx(1)*1.1 yy(1)*1.1 yy(end)*1.1 ] ); %zoom into area of interest
set(gca,'XTicklabel','');
set(gca,'YTicklabel',''); %remove axis tick lables




[x,y]=rotatedan2(273.13+Tdat,log(THdat)/log(10),thi,Y,minT,minTH,maxTH);
h=plot( x,(y),'kx'); %plot sounding


[a b]=size(th);
[a c]=size(THdat);
sati=find(abs(Pdat-850)<2);

clear tdiff;
for i=100:b %start at low pressure to avoid getting highest pressure value as cross-over point
    [a2 ay]=min(abs(th(i)-THdat));    
    tdiff(i-99)=abs(T(i)-273-Tdat(ay));
%     if abs(tdiff(i)-0)<1e-10 
%         tdiff(i)=NaN;
%     end
end

[a b]=min(tdiff);
b=b+99;
fprintf(1,'Cross-over pressure = %g mb\n',p(b)/100);
[c d]=min(abs(Pdat.*100-p(b)));
fprintf(1,'Cross-over height = %g m',dat(4,d));