

h=nanmedian(htt);
%h=cldtop(:,2);
h=nanmedian(cldtop(:,2));
kk=find(ht<h);
topi=max(kk);
cldthk=nanmean(cldtop(:,8))/12600*25/1e3;
kk=find(ht<h-cldthk);
boti=max(kk);
dkcld=min(25,topi-boti+2);


j=1;
dbb=-35:0;
jj=find(Ref(topi-2,:)<-20);
clear Pi ygg;
for dk=3:dkcld;
    jj=find(Ref(topi-dk,:)>-35);

%[P,S] = polyfit(Ref(topi-dk,jj),DV(topi-dk,jj)-DV(topi-2,jj),2);
[P,S] = polyfit(Ref(topi-dk,jj),DV(topi-dk,jj),2);

yg=polyval(P,dbb);
ygg(j,:)=yg;
Pi(j,:)=P;
j=j+1;
end
figure;plot(dbb,mean(ygg),dbb,polyval(mean(Pi),dbb),'o');

j=1;
velk=[];
for dk=1:dkcld;
vel = polyval(mean(Pi),Ref(topi-dk,:));
velk(j,:)=vel;
j=j+1;
end


vc=median(DV(topi-dkcld+1:topi,:)-velk);

%dk=dkcld;
Pg=[];
[Pg,Sg] = polyfit(mean(Ref(topi-dkcld:topi-floor(dkcld/2),:)),mean(DV(topi-dkcld:topi-floor(dkcld/2),:))-vc,2);
dbt=[ones(size(dbb')) exp(.06*dbb')];
tg=nanmean(Ref(topi-dkcld:topi-floor(dkcld/2),:))';vg=nanmean(DV(topi-dkcld:topi-floor(dkcld/2),:))'-vc';xp=[ones(size(tg)) exp(.06*tg)];aa=xp\vg;ygb=dbt*aa;
figure;plot(dbb,mean(ygg),dbb,polyval(mean(Pi),dbb),'o',dbb,polyval(Pg,dbb),'--',dbb,ygb,'-x',tg,vg,'.');xlabel('dBZ');ylabel('Vertical Motion (m/s)');

velk2=[];velk3=[];
j=1;
for dk=1:dkcld;
vel = polyval(Pg,Ref(topi-dk,:));
velk2(j,:)=vel;
tgg=Ref(topi-dk,:)';xpg=[ones(size(tgg)) exp(.06*tgg)];
velg=xpg*aa;
velk3(j,:)=velg';
j=j+1;
end

dk=5;vc2=mean(DV(topi-dkcld:topi-dk,:)-velk2(dk:dkcld,:));vc3=mean(DV(topi-dkcld:topi-dk,:)-velk3(dk:dkcld,:));

wp=[.02];ws= [.01];rp=3;rs=15;[n,wmn]=cheb1ord(wp,ws,rp,rs);
[b,a]=cheby1(n,rp,wmn);
vcfL=filter(b,a,vc-mean(vc));%\
vcfH=vc-vcfL-mean(vc);
vcfL3=filter(b,a,vc3-mean(vc3));%\
vcfH3=vc3-vcfL-mean(vc3);

read_Wband_motion
DV_mot=interp1(hhk,vkons(9,:),hh);
ii=find(isnan(DV_mot));DV_mot(ii)=0;

figure;
km=floor(dkcld/2);km=10;
ii=find(Ref(topi-km,:)>-50);
subplot(3,1,1);plot(hh(ii),DV(topi-km,ii)-velk2(km,ii),hh(ii),DV(topi-km,ii));axis([hh(100) hh(1200) -2 2]);ylabel('Wvelocity (m/s)');
subplot(3,1,2);plot(hh(ii),DV_mot(ii));axis([hh(100) hh(1200) -2 2]);ylabel('Wmotion (m/s)');
subplot(3,1,3);plot(hh(ii),DV(topi-km,ii)-velk2(km,ii)+DV_mot(ii));hold on;
plot(hh(ii),DV(topi-km-1,ii)-velk2(km-1,ii)+DV_mot(ii),'r');
plot(hh(ii),DV(topi-km-2,ii)-velk2(km-2,ii)+DV_mot(ii),'g');
axis([hh(100) hh(1200) -2 2]);ylabel('Wvelocity (m/s)');xlabel('Time (hrs)');
hold off;

fw=1/.3;
[S_w,F]=psd2(detrend(DV(topi-km,:)),length(hh),fw);
[S_ws,Fsl]=specsmoo(S_w,fw);
[S_w2,F]=psd2(detrend(DV(topi-km,ii)-velk2(km,ii)+DV_mot(ii)),length(hh),fw);
[S_ws2,Fsl]=specsmoo(S_w2,fw);

lss=length(Fsl);
figure;loglog(Fsl,S_ws.*Fsl,Fsl,S_ws2.*Fsl,Fsl(lss-41:lss),((Fsl(lss-41:lss)).^(-2/3))./(Fsl(lss-25).^(-2/3)).*S_ws2(lss-25).*Fsl(lss-25));
title([ 'VOCALS Wband Radar Spectrum:   Day ',num2str(ddd),'   Hour ',num2str(hh(1))]);
xlabel('Frequency in Hz');
ylabel('W-Spectrum * Frequency');

figure;plot(hh,vcfL);

j=1;Swi=[];kss=25;eps=[];Swif=[];
dkcld2=50;
for dk=1:dkcld2;
    veldbz = polyval(mean(Pg),Ref(topi-dk,:));
    velf=DV(topi-dk,:)+DV_mot-veldbz;%-velk2(dk,:)
    jj=find(Ref(topi-dk,:)<-35);
    velf(jj)=0;
    [S_w,F]=psd2(detrend(velf),length(hh),fw);
    [S_ws,Fsl]=specsmoo(S_w,fw);
    Swi(j,:)=S_ws;
    Swif(j,:)=S_ws.*Fsl;
    vls=.75/.54*(2*pi/12)^.667*Fsl(lss-25).^(5/3).*S_ws(lss-25); 
    eps(j)=vls^1.5;
    j=j+1;
end
figure;semilogx(eps,ht(topi-dkcld2:topi-1));;axis([ 1e-5 1e-3 0 1.5]);title([ 'VOCALS Wband Radar TKE Dissipation Rate:   Day ',num2str(ddd),'   Hour ',num2str(hh(1))]);
xlabel('\epsilon (m^2/s^3)');
ylabel('ALtitude (km)');

figure;loglog(Fsl,Swif);axis([1e-4 10 1e-4 1]);title([ 'VOCALS Wband Radar Spectrum:   Day ',num2str(ddd),'   Hour ',num2str(hh(1))]);
xlabel('Frequency in Hz');
ylabel('W-Spectrum * Frequency');


