% Find the cloud top from the sounding. Start at the inversion base
% (minimum of temperature) and look up in the next 100m for air that when
% mixed with cloud base (LCL) air will result in subsaturation. Test.
% Compare this level with the top from the radar wind profiler.
% SPdeS - 9 January 2007

% preamble
read_parameters;
thermo_constants
%clear t;

% construct title string, e.g. 'STRATUS 2007'
tit = cruise;
lt = length(cruise);
tit(lt+1) = ' ';
tit((lt+2):lt+5) = year; % This year is a string read in from the parameter file

% locate raw sounding files
aa = dir([way_raw_data_balloon '*_LATLON.txt']);
nsonde=length(aa);
% select subset based on yymmddhhMM dates
fst=str2mat(aa.name);
start=0710181500; % the date string as a number for identifying times
endt =0711082559; % include all of Stratus 2007
i1=max([1;find(str2num(fst(:,1:10))>=start,1,'first')]);
i2=min([find(str2num(fst(:,1:10))<=endt,1,'last');length(aa)]);
% convert filename timestamps to time variables
for i=i1:i2;
    balfn=aa(i).name; % simple choice of file for now
    yr(i) = str2double(balfn(1:2))+2000;
    month(i) = str2double(balfn(3:4));
    day (i)  = str2double(balfn(5:6));
    hour(i)  = str2double(balfn(7:8));
    minute(i) = str2double(balfn(9:10));
end

%load([way_raw_data_balloon(1:end-4) 'mat\wtec_edt_interp2p.mat'])
load([way_raw_data_balloon(1:end-4) 'mat\wtec_edt_interp2z.mat'])

la=nanmean(latz); lo=nanmean(lonz);
yday=datenum(yr,month,day,hour,minute,zeros(size(minute)))-datenum(2007,0,0,0,0,0);
x=repmat(la,[3 1]); y=repmat([500;550;NaN],[1 length(la)]);

% identify LCL from dewpoint depression dpd
dpd=tz-tdz;
% there's no cloud if the dewpoint never gets above T-0.5 C.
nocloud=min(dpd)>0.5;
dpda=dpd-repmat(min(dpd),size(tz,1),1);
% by inspection, LCL is lowest point with dpd>1 C below dpd<1
lcl=NaN+zeros(1,size(dpd,2));
m1=(1:length(dpd))'<250;
for sni=1:size(dpd,2)
    ci=find(dpda(:,sni)<=1 & m1,1,'last');
    if ~isempty(ci);
        while dpda(ci,sni)<1
            ci=ci-1;
        end
        lcl(sni)=ci;
    end
end
lcl(2)=lcl(3); % fill in where there's no LCL
thetalcl=diag(thetaz(lcl,:));
qlcl=diag(qz(lcl,:));
hlcl=hz(lcl);

% inversion base height from sonde T minimum
mask=repmat(hz(1:end-1)'>500 & hz(1:end-1)'<2000,[1 size(tz,2)]);
[dTdz,iiz]=max(mask.*diff(tz));
for ni=1:size(tz,2)
    mask=(hz<=hz(iiz(ni)))';
    [tinv(ni),iinvz(ni)]=min(1e12*~mask+tz(:,ni));
end
zinv=hz(iinvz);

% candidate top parcels 0-250m above inversion base
ican=repmat(iinvz,[26 1])+repmat((0:25)',[1 length(iinvz)]);
pcan=zeros(size(ican));
thetacan=pcan;
qcan=pcan;
for i=1:length(ican)
    thetacan(:,i)=thetaz(ican(:,i),i);
    qcan(:,i)=qz(ican(:,i),i);
    pcan(:,i)=pz(ican(:,i),i);
end
% mixtures of top and LCL parcels
% thetal is conserved in the cloud; at LCL and top theta=thetal
frac=0.25; % fractional dilution of the updraft
thetamix=frac*thetacan+(1-frac)*repmat(thetalcl',[size(ican,1) 1]);
tmix=temp(pcan*1e2,thetamix)-273.15; % degree C
qmix=frac*qcan+(1-frac)*repmat(qlcl',[size(ican,1) 1]);
qsmix=reshape(qsat([tmix(:),pcan(:)])',size(qmix));

% find first level where mixed parcel evaporates
elect=zeros([1 length(qmix)]);
for i=1:length(qmix)
    e=find(qmix(:,i)<qsmix(:,i),1,'first');
    if ~isempty(e)
        elect(i)=e;
    end
end

% intersection where parcel evaporates
thetainv=NaN+zeros(1,length(elect));
qinv=thetainv;
pinv=thetainv;
hinv=thetainv;
isinv=elect>0;
thetainv(isinv)=diag(thetacan(elect(isinv),isinv));
qinv(isinv)=diag(qcan(elect(isinv),isinv));
pinv(isinv)=diag(pcan(elect(isinv),isinv));
hinv(isinv)=hz(iinvz(isinv)+elect(isinv)-1);
mask=all(qmix<qsmix); % all mixtures unsaturated
qinv(mask)=0;
hinv(mask)=0;
pinv(mask)=0;
thetainv(mask)=0;
qevap(isinv)=diag(qmix(elect(isinv),isinv));
thetaevap(isinv)=diag(thetamix(elect(isinv),isinv));
qevap(~isinv)=NaN; thetaevap(~isinv)=NaN;
thetaevap(mask)=NaN;
qevap(mask)=NaN;
if 0 % plot intersection
    for np=1:12
        subplot(4,3,np,'align')
        plot(thetamix(:,(np-1)*6+(1:6)),qmix(:,(np-1)*6+(1:6)))
        hold on
        plot(thetamix(:,(np-1)*6+(1:6)),qsmix(:,(np-1)*6+(1:6)))
        plot(thetaevap((np-1)*6+(1:6)),qevap((np-1)*6+(1:6)),'.');
    end
end

% load profiler PBL top
top=load([way_proc_data_profiler cruise year '_top_all.txt']);
% load ceilo cloud base
base=load([way_proc_data_ceilo cruise year 'ceilo_10min_a.txt']);
% load surface met values to compute Tlcl from TRH at 15.5m
sfc=load([way_proc_data_flux 'flux_5hf_STRATUS07.txt']);
tlclsfc=Tlcl(101500,sfc(:,7)+273.15,sfc(:,9)/1e3);
zlclsfc=(sfc(:,7)+273.15-tlclsfc)/9.8; % km
% heights from 296K and 4 g/kg contours
%hth296=contourc(yday(i1:i2),hz(iz)/1e3,thetaz(iz,i1:i2),[296 296]);
%hqv4=contourc(yday(i1:i2),hz(iz)/1e3,qz(iz,i1:i2),[4 4]);
for i=i1:i2;
    [junk,iun]=unique(thetaz(:,i));
    isf=find(isfinite(thetaz(:,i)) & hz'<2100);
    ii=sort(iun(ismember(iun,isf)));
    hth296(i)=interp1(thetaz(ii,i),hz(ii),296);
    [junk,iun]=unique(qz(:,i));
    isf=find(isfinite(qz(:,i)) & hz'<2100);
    ii=sort(iun(ismember(iun,isf)));
    hqv4(i)=interp1(qz(ii,i),hz(ii),4);
end

% plot time-height series for the whole cruise
iz=1:251;
clf
b2rcolormap(24);
subplot(2,1,1,'align')
[c,h]=contourf(yday(i1:i2),hz(iz)/1e3,thetaz(iz,i1:i2),284:2:312);
set(h,'edgecolor','none')
set(gca,'ylim',[0 2.5])
hold on
plot(top(:,1)+top(:,2)/86400,top(:,3),'m.','markersize',2)
plot(yday(i1:i2),hinv(i1:i2)/1e3,'r.')
plot(base(:,1),base(:,13)/1e3,'g.','markersize',2)
plot(yday(i1:i2),zinv(i1:i2)/1e3,'wx')
plot(yday(i1:i2),hlcl(i1:i2)/1e3,'r.')
plot(sfc(:,1),zlclsfc,'y')
contour(yday(i1:i2),hz(iz)/1e3,qz(iz,i1:i2),[4 4],'k')
contour(yday(i1:i2),hz(iz)/1e3,thetaz(iz,i1:i2),[296 296],'g')
%plot(yday,hth296/1e3,'k.')
caxis([284 312])
colorbar
title({'STRATUS 2007 time-height section' 'potential temperature (K)'})
subplot(2,1,2,'align')
[c,h]=contourf(yday(i1:i2),hz(iz)/1e3,qz(iz,i1:i2),0:1:13);
set(h,'edgecolor','none')
set(gca,'ylim',[0 2.5])
hold on
plot(top(:,1)+top(:,2)/86400,top(:,3),'m.','markersize',2)
plot(yday(i1:i2),hinv(i1:i2)/1e3,'r.')
plot(base(:,1),base(:,13)/1e3,'g.','markersize',2)
plot(yday(i1:i2),zinv(i1:i2)/1e3,'kx')
plot(yday(i1:i2),hlcl(i1:i2)/1e3,'r.')
plot(sfc(:,1),zlclsfc,'b')
contour(yday(i1:i2),hz(iz)/1e3,qz(iz,i1:i2),[4 4],'k')
contour(yday(i1:i2),hz(iz)/1e3,thetaz(iz,i1:i2),[296 296],'g')
%plot(yday,hqv4/1e3,'k.')
caxis([0 11])
colorbar
title('specific humidity (g/kg)');
xlabel('year day (UTC, Oct 18-Nov 3)')
orient landscape
%print('-depsc2',[way_proc_images_balloon 'compare_sounding_profilertop_ceilobase.eps'])
%print('-djpeg',[way_proc_images_balloon 'compare_sounding_profilertop_ceilobase.jpg'])

% time series of cloud/MABL top and base within 1 hour of soundings.
for ii=1:length(yday)
    mask=abs(yday(ii)-top(:,1)+top(:,2)/86400)<=1/24;
    prof_top(ii)=nanmean(top(mask,3));
    mask=abs(yday(ii)-sfc(:,1))<=1/24;
    lcl_base(ii)=nanmean(zlclsfc(mask));
    mask=abs(yday(ii)-base(:,1))<=1/24 & base(:,13)<2e3;
    ceil_base(ii)=nanmean(base(mask,13))/1e3;
end
[zi,si]=sort(zinv);

% correlations of different measures of cloud top/inv. height
Mt=[hqv4(yday>293.5)',hth296(yday>293.5)',zinv(yday>293.5)',prof_top(yday>293.5)'*1e3];
Ct=corr(Mt);
% correlations of different measures of cloud base
Mb=[lcl_base(yday>293.5)' hlcl(yday>293.5)'/1e3 ceil_base(yday>293.5)'];
Cb=corr(Mb);
% results copied to My Documents\Reports\cloud_base+top.xcl
% Comparison summarized in My Documents\Reports\Vertical cloud...
% distribution in STRATUS 2007.doc
P=polyfit(zinv(zinv>700),prof_top(zinv>700)*1e3,1);

plot(zinv,prof_top*1e3,'.')
hold on
plot([200,2000],[200,2000],'k-')
ylabel('height of max. profiler SNR (m)')
xlabel('height of min. sonde temperature (m)')
rectangle('position',[200 700 500 1700])
plot([200 2000],P(1)*[200; 2000]+P(2),'r')
text(1200,600,sprintf('y=%2.3g x + %3.3g',P),'color','r')
text(1200,700,'mean offset=121 m')
axis([200 2000 200 2000])
axis square
%print('-depsc2',[way_proc_images_balloon 'sonde_prof_zinv.eps'])
%print('-djpeg',[way_proc_images_balloon 'sonde_prof_zinv.jpeg'])
