%% Produce figures showing the KHI-billows during PCAPS IOP-1

close all
clear all

%% netcdf handling
%add a bunch of paths to make sure netcdf handling is easier
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/'); %we call some scripts from here
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/mexcdf/mexnc/') % netcdf
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/mexcdf/snctools'); %netcdf

%% define times of interest

for IOP=[1]%3 4 5 6 9]
    if IOP ==1
        starter=datenum(2010,12,2,0,00,00);
        ender=datenum(2010,12,4,0,00,00);
        
    elseif IOP==2
        starter=datenum(2010,12,7,00,00,00);
        ender=datenum(2010,12,11,0,00,00);
    elseif IOP==3
        starter=datenum(2010,12,11,00,00,00);
        ender=datenum(2010,12,15,0,00,00);   
    elseif IOP==4
        starter=datenum(2010,12,23,18,00,00);
        ender=datenum(2010,12,26,19,00,00);   
    elseif IOP==5
        starter=datenum(2010,12,31,18,00,00);
        ender=datenum(2011,1,10,0,00,00);       
    elseif IOP==6
        starter=datenum(2011,1,10,00,00,00);
        ender=datenum(2011,1,18,0,00,00);
    elseif IOP==7
        starter=datenum(2011,1,19,00,00,00);
        ender=datenum(2011,1,23,0,00,00);
    elseif IOP==8
        starter=datenum(2011,1,22,00,00,00);
        ender=datenum(2011,1,27,0,00,00);
    elseif IOP==9
        starter=datenum(2011,1,25,18,00,00);
        ender=datenum(2011,1,31,16,00,00);  
    elseif IOP==10
        starter=datenum(2011,2,1,00,00,00);
        ender=datenum(2011,2,6,0,00,00);
    end
end


%% Ceilometer data

% point to Joe's ceilometer data
addpath('/uufs/chpc.utah.edu/common/home/whiteman-group1/jyoung/hdf5/pcaps/');
fname='ncar_ceil.h5';
path='/uufs/chpc.utah.edu/common/home/whiteman-group1/jyoung/hdf5/pcaps/ncar_ceil.h5';


time = h5read(path,'/time');
tkey=time.key;
time=double(time.time)./86400 + datenum(1970,1,1);
tidx=find(time>=starter & time<=ender);
time=time(tidx);
height=h5read(path,'/height');
height=height+1320;
hidx=find(height < 2500 & height>1340);
height=height(hidx);
bsnow=h5read(path,'/bs',[hidx(1) tidx(1)],[length(hidx) length(tidx)]);
bsnow=double(bsnow);

starter2=datenum(2010,12,3,3,50,0);
%starter2=datenum(2010,12,3,4,40,0);

%now bring in RASS/RWP data
f2='PCAPS_RWP.nc';
nc_dump(f2);
uRWP=nc_varget(f2,'uRWP');
vRWP=nc_varget(f2,'vRWP');
tRWP=nc_varget(f2,'timeRWP');
hRWP=nc_varget(f2,'heightRWP');
[~,rwptidx]=min(abs(tRWP-starter2) );

%% RASS
file='PCAPS_RASS.nc';
nc_dump(file)

thetaRASS=nc_varget(file,'ThetaRASS');
timeRASS=nc_varget(file,'timeRASS');
heightRASS=nc_varget(file,'heightRASS');
[~,rasstidx]=min(abs(timeRASS-starter2) );
uRASS=interp2(tRWP',hRWP,uRWP,timeRASS',heightRASS);
vRASS=interp2(tRWP',hRWP,vRWP,timeRASS',heightRASS);

%%

scrz=get(0,'screensize');
f=figure(1);
clf
set(f,'Renderer','Painters');
set(f,'color','w','Units','pixels','pos',[1200 100 1600 900]);
sp1=subplot(3,6,7);
spos1=get(sp1,'pos');
spos1(3)=1.5*spos1(3);
spos1(1)=spos1(1)-.075;
set(sp1,'pos',spos1);
spos1=get(sp1,'pos');
redge=spos1(1)+spos1(3);
plot(thetaRASS(:,rasstidx),heightRASS,'-k','linewidth',2);
ylim([1330 1700]);
xlim([284 304])
%grid on;
xlabel('Pot. Temp. [K]','fontsize',10,'fontweight','bold');
ylabel('Height [m]','fontsize',10,'fontweight','bold');
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold')
text(284.5,1670,'b)','backgroundcolor',[.7 .7 .7],'edgecolor','k','fontsize',10,'fontweight','bold');

grid on;
%% 
ender2=starter2+(.5/24);
sp2=subplot(3,6,8:10);
spos2=get(sp2,'pos');
spos2(1)=redge+.01;
spos2(3)=1.1*spos2(3);
set(sp2,'pos',spos2);
spos2=get(sp2,'pos');
redge=spos2(1)+spos2(3);
cmap=[gray(8);flipud(gray(8))];
cmap(1:8,3)=1;
cmap(9:end,1)=1;
colormap(cmap);
colormap((jet(24)))
colormap(flipud(pink(24)))

tidx2=find(time>=starter2 & time<=ender2);

caxis([-7.75 -6.3]);
%caxis([-7.75 -7]);
hold on
bssmooth=interp2(bsnow(:,tidx2),2);
[tmat1,hmat1]=meshgrid(time(tidx2),height);
tmat=interp2(tmat1,2);
hmat=interp2(hmat1,2);
pcolor(1:size(bssmooth,2),(hmat(:,1))',bssmooth);shading flat;
%pcolor(tmat1,hmat1,bsnow(:,tidx2));shading interp;
axis tight;
ylim([1330 1700]);

xticker=get(gca,'xtick');
xlab=tmat(1,xticker);
set(gca,'xticklabel',{datestr(xlab,'HH:MM')});
xlabel('Time HH:MM [UTC]','fontsize',10,'fontweight','bold');
set(gca,'yticklabel',[]);
box on;
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold')
cbh=colorbar('north');cpos=get(cbh,'pos');cpos(3)=.2.*cpos(3);set(cbh,'pos',cpos);

%%

sp3=subplot(3,6,11); hold on;
spos3=get(sp3,'pos');
spos3(1)=redge+.01;
set(sp3,'pos',spos3);
spos3=get(sp3,'pos');
redge=spos3(1)+spos3(3);
plot(uRASS(:,rasstidx),heightRASS,'-k','linewidth',2);
plot(vRASS(:,rasstidx),heightRASS,'r','linewidth',2);
ylim([1330 1700]);
xlim([-5 17])
%grid on;
set(gca,'yticklabel',[],'xtick',[-5:5:20]);
legend('u','v');
xlabel('Wind Speed [m/s]','fontsize',10,'fontweight','bold');
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold')
grid on;
box on;
%% RI

dudz=(uRASS(:,rasstidx)-circshift(uRASS(:,rasstidx),1))./60;
dvdz=(vRASS(:,rasstidx)-circshift(vRASS(:,rasstidx),1))./60;
dtdz=(thetaRASS(:,rasstidx)-circshift(thetaRASS(:,rasstidx),1))./60;
RIG=(9.81./(thetaRASS(:,rasstidx))).*(dtdz)./(dudz.^2 + dvdz.^2);
sp6=subplot(3,6,12);
spos6=get(sp6,'pos');
spos6(1)=redge+.01;
set(sp6,'pos',spos6);
plot(RIG,heightRASS,'-k','linewidth',2)
ylim([1330 1700]);
xlim([-.1 4]);
%grid on;
set(gca,'yaxislocation','right','xtick',[0 1 2 3 4]);
xlabel('Ri','fontsize',10,'fontweight','bold');
ylabel('Height [m]','fontsize',10,'fontweight','bold');
box on;
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold')
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SAME thing but using sounding data near 12 UTC 3 December

file='/uufs/chpc.utah.edu/common/home/horel-group/pcaps.utah.edu/data/SONDE_PROC/20101203/NCAR_QC_2010120312.dat';


[b,c]=butter(5,1/18/.5,'low');
[b2,c2]=butter(5,1/3/.5);

lehi=importdata(file,' ',5);
sounding=lehi.data;
sounding(sounding == -999)= nan;
hgts=sounding(4:end,1);
theta=sounding(4:end,9);
theta=filtfilt(b2,c2,theta);
dtdz=gradient(theta,10);
dtdz=filtfilt(b2,c2,dtdz);
spd=sounding(4:end,12);
u=sounding(4:end,14);
u=inpaint_nans(u);
usm=filtfilt(b,c,u);
dudz=gradient(usm,10);

v=sounding(4:end,15);
v=inpaint_nans(v);

vsm=filtfilt(b,c,v);

dvdz=gradient(vsm,10);
%dudz=(u-circshift(u,1))./10;
%dvdz=(v-circshift(v,1))./10;
%dtdz=(theta-circshift(theta,1))./10;
dtdz=gradient(theta,10);
nsqr=(9.81./(theta)).*(dtdz);
RIG=nsqr./(dudz.^2 + dvdz.^2);
RIG(isnan(RIG))=0;
% 



RICH=sounding(4:end,end);RICH(RICH>100)=nan;

%richardson number tests
% figure(88);clf
% axx=tight_subplot(2,3);
% axes(axx(1))
% plot(theta,hgts); ylim([1300 2500]);
% axes(axx(4));
% plot(nsqr,hgts); ylim([1300 2500]);
% axes(axx(2))
% hold on;
% plot(vsm,hgts,'k','linewidth',2)
% plot(v,hgts,'k','linewidth',1);
% plot(usm,hgts,'r','linewidth',2)
% plot(u,hgts,'r'); ylim([1300 2500]);
% axes(axx(5));
% plot((dudz.^2 + dvdz.^2),hgts);ylim([1300 2500]);
% axes(axx(6));
% hold on;
% plot(RIG,hgts,'r'); ylim([1300 2500]);
% plot(RICH,hgts,'k'); xlim([-.1 5]);
% axes(axx(3));
% RInorm=RIG./1; 
% shearnorm=(dudz.^2 + dvdz.^2)./.02;
% nnorm=nsqr./(35*10^-4);
% plot(RInorm,hgts);
% hold on;
% plot(shearnorm,hgts,'r');
% plot(nnorm,hgts,'k');
% ylim([1300 2500]);
% xlim([0 5]);
%%

sp1=subplot(3,6,13);
spos1=get(sp1,'pos');
spos1(3)=1.5*spos1(3);
spos1(1)=spos1(1)-.075;
set(sp1,'pos',spos1);
spos1=get(sp1,'pos');
redge=spos1(1)+spos1(3);
plot(theta,hgts,'-k','linewidth',2)
ylim([1330 1700]);
xlim([284 303])
grid on;
box on;
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold');
xlabel('Pot. Temp. [K]','fontsize',10,'fontweight','bold');
text(284.5,1670,'c)','backgroundcolor',[.7 .7 .7],'edgecolor','k','fontsize',10,'fontweight','bold');
ylabel('Height [m]','fontsize',10,'fontweight','bold');


%% loop to look for good KH waves
starter2=datenum(2010,12,3,10,50,0);

ender2=starter2+(.5/24);
figure(1)

sp2=subplot(3,6,14:16);
spos2=get(sp2,'pos');
spos2(1)=redge+.01;
spos2(3)=spos2(3)*1.1;
set(sp2,'pos',spos2);
spos2=get(sp2,'pos');
redge=spos2(1)+spos2(3);
cmap=[gray(8);flipud(gray(8))];
cmap(1:8,3)=1;
cmap(9:end,1)=1;
colormap(cmap);
colormap((jet(24)))
colormap(flipud(pink(24)))

tidx2=find(time>=starter2 & time<=ender2);
caxis([-7.75 -6.3]);
hold on
bssmooth=interp2(bsnow(:,tidx2),2);
[tmat1,hmat1]=meshgrid(time(tidx2),height);
tmat=interp2(tmat1,2);
hmat=interp2(hmat1,2);
colormap(pink(48));
pcolor(1:size(bssmooth,2),(hmat(:,1))',bssmooth);shading interp;
axis tight;
ylim([1330 1700]);
xticker=get(gca,'xtick');
xlab=tmat(1,xticker);
set(gca,'xticklabel',{datestr(xlab,'HH:MM')});
xlabel('Time HH:MM [UTC]','fontsize',10,'fontweight','bold');
%cbh=colorbar('north');cpos=get(cbh,'pos');cpos(3)=.2.*cpos(3);set(cbh,'pos',cpos);

set(gca,'fontsize',10,'fontweight','bold','box','on','layer','top','linewidth',1);
set(gca,'yticklabel',[]);

%%
sp3=subplot(3,6,17); hold on;
spos3=get(sp3,'pos');
spos3(1)=redge+.01;
set(sp3,'pos',spos3);
spos3=get(sp3,'pos');
redge=spos3(1)+spos3(3);
plot(usm,hgts,'-k','linewidth',2)
plot(vsm,hgts,'-r','linewidth',2)
ylim([1330 1700]);
xlim([-5 17])
set(gca,'yticklabel',[],'xtick',[-5:5:20]);
legend('u','v');
xlabel('Wind Speed [m/s]','fontsize',10,'fontweight','bold');
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold')
grid on;
box on;
%% RI


sp6=subplot(3,6,18);
spos6=get(sp6,'pos');
spos6(1)=redge+.01;
set(sp6,'pos',spos6); hold on;
%plot(RICH,hgts,'linewidth',2)
plot(RIG,hgts,'-k','linewidth',2);
ylim([1330 1700]);
xlim([-.1 4]);
set(gca,'yaxislocation','right','xtick',[0 1 2 3 4]);
xlabel('Ri','fontsize',10,'fontweight','bold');
ylabel('Height [m]','fontsize',10,'fontweight','bold');
box on;
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold')
grid on;
%export_fig khi_2010_12_3_11.png -m2 -painters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SAME thing but using sounding data near 12 UTC 2 December

file='/uufs/chpc.utah.edu/common/home/horel-group/pcaps.utah.edu/data/SONDE_PROC/20101202/cpin_KSLC_2010120212_smth.dat';
%file='/uufs/chpc.utah.edu/common/home/horel-group/pcaps.utah.edu/data/SONDE_PROC/20101202/cpin_NCAR_QC_2010120206.dat';


lehi=importdata(file,' ',5);
sounding=lehi.data;
sounding(sounding == -999)= nan;
hgts=sounding(:,1);
theta=sounding(:,9);
dtdz=gradient(theta,10);
spd=sounding(:,12);
u=sounding(:,14);
dudz=gradient(u,10);
v=sounding(:,15);
dvdz=gradient(v,10);
%dudz=(u-circshift(u,1))./10;
%dvdz=(v-circshift(v,1))./10;
%dtdz=(theta-circshift(theta,1))./10;
dtdz=gradient(theta,10);
Nsqr=(9.81./(theta)).*(dtdz);
shear=(dudz.^2 + dvdz.^2);
RIG=Nsqr./shear;
RICH=sounding(:,end-2);


%%


% scrz=get(0,'screensize');
% f=figure(3);
% clf
% set(f,'color','w','Units','pixels','pos',[1200 100 1600 300],'render','painters');
sp1=subplot(3,6,1);
spos1=get(sp1,'pos');
spos1(3)=1.5*spos1(3);
spos1(1)=spos1(1)-.075;
set(sp1,'pos',spos1);
spos1=get(sp1,'pos');
redge=spos1(1)+spos1(3);
plot(theta,hgts,'-k','linewidth',2)
ylim([1330 1700]);
xlim([280 293])
grid on;
box on;
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold');
xlabel('Pot. Temp. [K]','fontsize',10,'fontweight','bold');
text(280.5,1670,'a)','backgroundcolor',[.7 .7 .7],'edgecolor','k','fontsize',10,'fontweight','bold');


%% loop to look for good KH waves
starter2=datenum(2010,12,2,12,0,0);

ender2=starter2+(.5/24);
%figure(1)

sp2=subplot(3,6,2:4);
spos2=get(sp2,'pos');
spos2(1)=redge+.01;
spos2(3)=spos2(3)*1.1;
set(sp2,'pos',spos2);
spos2=get(sp2,'pos');
redge=spos2(1)+spos2(3);
cmap=[gray(8);flipud(gray(8))];
cmap(1:8,3)=1;
cmap(9:end,1)=1;
colormap(cmap);
colormap((jet(24)))
colormap(flipud(pink(24)))
tidx2=find(time>=starter2 & time<=ender2);

caxis([-7.75 -6.3]);
hold on
bssmooth=interp2(bsnow(:,tidx2),2);
[tmat1,hmat1]=meshgrid(time(tidx2),height);
tmat=interp2(tmat1,2);
hmat=interp2(hmat1,2);
pcolor(1:size(bssmooth,2),(hmat(:,1))',bssmooth);shading interp;
axis tight;
ylim([1330 1700]);
xticker=get(gca,'xtick');
xlab=tmat(1,xticker);
set(gca,'xticklabel',{datestr(xlab,'HH:MM')});
xlabel('Time HH:MM [UTC]','fontsize',10,'fontweight','bold');
set(gca,'fontsize',10,'fontweight','bold','box','on','layer','top','linewidth',1);
set(gca,'yticklabel',[]);

%%
sp3=subplot(3,6,5); hold on;
spos3=get(sp3,'pos');
spos3(1)=redge+.01;
set(sp3,'pos',spos3);
spos3=get(sp3,'pos');
redge=spos3(1)+spos3(3);
plot(u,hgts,'-k','linewidth',2)
plot(v,hgts,'-r','linewidth',2)
ylim([1330 1700]);
xlim([-5 17])
set(gca,'yticklabel',[],'xtick',[-5:5:20]);
legend('u','v');
xlabel('Wind Speed [m/s]','fontsize',10,'fontweight','bold');
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold')
grid on; box on; 
%% RI
sp6=subplot(3,6,6);
spos6=get(sp6,'pos');
spos6(1)=redge+.01;
set(sp6,'pos',spos6);
plot(RIG,hgts,'-k','linewidth',2)
hold on;
%plot(RICH,hgts,'--k','linewidth',2);
ylim([1330 1700]);
xlim([-.1 4]);
set(gca,'yaxislocation','right','xtick',[0 1 2 3 4]);
xlabel('Ri','fontsize',10,'fontweight','bold');
ylabel('Height [m]','fontsize',10,'fontweight','bold');
box on;
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold')
grid on;
return
%export_fig khi_2010_12_2_12.png -m2 -painters

%return;

%%


file='/uufs/chpc.utah.edu/common/home/horel-group/pcaps.utah.edu/data/SONDE_PROC/20101203/NCAR_QC_2010120306.dat';

lehi=importdata(file,' ',5);
sounding=lehi.data;
sounding(sounding == -999)= nan;
hgts=sounding(:,1);
theta=sounding(:,9);
dtdz=gradient(theta,10);
spd=sounding(:,12);
u=sounding(:,14);
dudz=gradient(u,10);
v=sounding(:,15);
%dvdz=gradient(v,10);
dudz=(u-circshift(u,1))./10;
dvdz=(v-circshift(v,1))./10;
dtdz=(theta-circshift(theta,1))./10;
RIG=(9.81./nanmean(theta)).*(dtdz)./(dudz.^2 + dvdz.^2);

RICH=sounding(:,end);



%%


scrz=get(0,'screensize');
f=figure(5);
clf
set(f,'color','w');%,'Units','pixels','pos',[1200 100 1600 400]);
sp1=subplot(1,6,1);
spos1=get(sp1,'pos');
spos1(3)=1.2*spos1(3);
spos1(1)=spos1(1)-.05;
set(sp1,'pos',spos1);
spos1=get(sp1,'pos');
redge=spos1(1)+spos1(3);
plot(theta,hgts,'-k','linewidth',2)
ylim([1320 1600]);
xlim([289 303])
grid on;
box on;
xlabel('Pot. Temp. [K]','fontsize',10,'fontweight','bold');
set(sp1,'box','on','layer','top','linewidth',1,'fontsize',10,'fontweight','bold');


%% loop to look for good KH waves
starter2=datenum(2010,12,3,03,50,0);
for i=1;%1:10;
ender2=starter2+(1.5/24);
figure(5)

sp2=subplot(1,6,2:4);
spos2=get(sp2,'pos');
spos2(1)=redge;
spos2(3)=1.2*spos2(3);
set(sp2,'pos',spos2);
spos2=get(sp2,'pos');
redge=spos2(1)+spos2(3);
cmap=[gray(8);flipud(gray(8))];
cmap(1:8,3)=1;
cmap(9:end,1)=1;
colormap(cmap);
colormap((jet(24)))
tidx2=find(time>=starter2 & time<=ender2);
%pcolor(1:length(tidx2),height,bsnow(:,tidx2));shading interp; 

caxis([-7.75 -6.3]);
hold on
bssmooth=interp2(bsnow(:,tidx2),1);
[tmat,hmat]=meshgrid(time(tidx2),height);
tmat=interp2(tmat,1);
hmat=interp2(hmat,1);
pcolor(1:size(bssmooth,2),(hmat(:,1))',bssmooth);shading flat;%interp;
axis tight;
ylim([1320 1600]);
xticker=get(gca,'xtick');
xlab=tmat(1,xticker);%time(tidx2(xticker));
set(gca,'xticklabel',{datestr(xlab,'HH:MM')});
xlabel('Time HH:MM [UTC]','fontsize',10,'fontweight','bold');
set(gca,'fontsize',10,'fontweight','bold','box','on','layer','top','linewidth',1);
set(gca,'yticklabel',[]);


end

%%
sp3=subplot(1,6,5); hold on;
spos3=get(sp3,'pos');
spos3(1)=redge;
set(sp3,'pos',spos3);
spos3=get(sp3,'pos');
redge=spos3(1)+spos3(3);
plot(u,hgts,'-k','linewidth',2)
plot(v,hgts,'-r','linewidth',2)
ylim([1320 1600]);
xlim([-4 20])
grid on;
box on;
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold');
set(gca,'yticklabel',[],'xaxislocation','top');
legend('u','v');
xlabel('Wind Speed [m/s]','fontsize',10,'fontweight','bold');
%% RI


sp6=subplot(1,6,6);
spos6=get(sp6,'pos');
spos6(1)=redge;
set(sp6,'pos',spos6);
plot(RICH,hgts,'linewidth',2)
ylim([1320 1600]);
xlim([-.1 4]);
%set(gca,'xscale','log')
grid on;
box on;
set(gca,'layer','top','linewidth',1,'fontsize',10,'fontweight','bold')
set(gca,'yaxislocation','right');
xlabel('Ri','fontsize',10,'fontweight','bold');

%% Overview
starter2=datenum(2010,12,2,11,00,0);
for i=1;%1:10;
ender2=starter2+(2/24);
f=figure(10);clf; set(f,'color','w','units','pixels','pos',[1200 100 1600 400]);
subplot(2,1,1)

cmap=[gray(8);flipud(gray(8))];
cmap(1:8,3)=1;
cmap(9:end,1)=1;
colormap(cmap);
colormap((jet(24)))
tidx2=find(time>=starter2 & time<=ender2);
%pcolor(1:length(tidx2),height,bsnow(:,tidx2));shading interp; 

caxis([-7.75 -6.3]);
hold on
bssmooth=interp2(bsnow(:,tidx2),1);
[tmat,hmat]=meshgrid(time(tidx2),height);
tmat=interp2(tmat,1);
hmat=interp2(hmat,1);
pcolor(1:size(bssmooth,2),(hmat(:,1))',bssmooth);shading flat;%interp;
%pcolor(1:length(bssmooth),height,bssmooth);shading interp;
%contourf(1:length(tidx2),height,bsnow(:,tidx2),[-10:.005:-6.5],'linestyle','none');
axis tight;
ylim([1320 1900]);
xticker=get(gca,'xtick');
xlab=tmat(1,xticker);%time(tidx2(xticker));
set(gca,'xticklabel',{datestr(xlab,'HH:MM')});
xlabel('Time HH:MM [UTC]','fontsize',10,'fontweight','bold');
%ylabel('Height [m]','fontsize',10,'fontweight','bold');
set(gca,'fontsize',10,'fontweight','bold','box','on','layer','top','linewidth',1);
set(gca,'yticklabel',[]);


end


%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%% 6 hr increments
%starters=[datenum(2010,12,2,12,10,0)  datenum(2010,12,3,3,50,0) datenum(2010,12,3,10,45,0)]
starters=[datenum(2010,12,2,12,0,0)  datenum(2010,12,3,3,50,0) datenum(2010,12,3,10,45,0)]


figure(55); clf

%ax=tight_subplot(3,1)  
for i=1:3;%1:10;
starter2=starters(i);%starter2+6/24;
ender2=starter2+(.5/24);

%axes(ax(i))
subplot(3,1,i)
cmap=[gray(8);flipud(gray(8))];
cmap(1:8,3)=1;
cmap(9:end,1)=1;
colormap(cmap);
colormap((jet(24)))
colormap(flipud(pink(24)))
tidx2=find(time>=starter2 & time<=ender2);
%pcolor(1:length(tidx2),height,bsnow(:,tidx2));shading interp; 

caxis([-7.75 -6.3]);
hold on
bssmooth=interp2(bsnow(:,tidx2),2);
[tmat,hmat]=meshgrid(time(tidx2),height);
tmat=interp2(tmat,2);
hmat=interp2(hmat,2);
pcolor(1:size(bssmooth,2),(hmat(:,1))',bssmooth);shading flat;%interp;
axis tight;
ylim([1320 1650]);
xticker=get(gca,'xtick');
xlab=tmat(1,xticker);%time(tidx2(xticker));
set(gca,'xticklabel',{datestr(xlab,'HH:MM')});
xlabel('Time HH:MM [UTC]','fontsize',10,'fontweight','bold');
set(gca,'fontsize',10,'fontweight','bold','box','on','layer','top','linewidth',1);
set(gca,'yticklabel',[]);


end

%% profiles for the above

starters=[datenum(2010,12,2,12,10,0)  datenum(2010,12,3,3,50,0) datenum(2010,12,3,10,45,0)]



f=figure(1000);
clf
axx=tight_subplot(3,1)  
clrs={'-.k','--k','-k'};
clrs2={'-.b','--b','-b'};
clrs3={'-.r','--r','-r'};
for i=1:3;%1:10;
starter2=starters(i);
[~,rasstidx]=min(abs(timeRASS-starter2) );
ender2=starter2+(.5/24);

axes(axx(1))
hold on;
plot(thetaRASS(:,rasstidx),heightRASS,clrs{i},'linewidth',2);
axis tight
ylim([1330 2200]);
xlabel('Pot. Temp. [K]','fontsize',10,'fontweight','bold');
ylabel('Height [m]','fontsize',10,'fontweight','bold');
grid on;
xlim([275 305])

axes(axx(2))
hold on;
plot(uRASS(:,rasstidx),heightRASS,clrs2{i},'linewidth',2);
plot(vRASS(:,rasstidx),heightRASS,clrs3{i},'linewidth',2);
ylim([1330 2200]);
xlim([-5 17])
grid on;

dudz=(uRASS(:,rasstidx)-circshift(uRASS(:,rasstidx),1))./60;
dvdz=(vRASS(:,rasstidx)-circshift(vRASS(:,rasstidx),1))./60;
dtdz=(thetaRASS(:,rasstidx)-circshift(thetaRASS(:,rasstidx),1))./60;
RIG=(9.81./(thetaRASS(:,rasstidx))).*(dtdz)./(dudz.^2 + dvdz.^2);

axes(axx(3))
hold on;
plot(RIG,heightRASS,clrs{i},'linewidth',2)
ylim([1330 2200]);
xlim([-.1 10]);
grid on;



end
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%

return
%% 
R=287;
g=9.81;

fname='ISSSF_PCAPS.txt';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 
ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(5));

lat=str2double(ST(3));
lon=str2double(ST(4));
data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
tvec=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));

% since there are some missing times lets make a complete time vector
SFCtimes=tvec(1):(1/60/24):tvec(end);

[~,midx]=setdiff(SFCtimes,tvec);


%interp data to fill the voids
temp=data(:,7)+273.15;
temp=interp1(tvec,temp,SFCtimes);
pres=data(:,11);
pres=interp1(tvec,pres,SFCtimes);


u = -1 .* data(:,end-6) .* sin ( (pi/180) .* data(:,end-5) );
v = -1.* data(:,end-6) .* cos ( (pi/180) .* data(:,end-5) );
u=interp1(tvec,u,SFCtimes);
v=interp1(tvec,v,SFCtimes);

u=inpaint_nans(u);
v=inpaint_nans(v);

tidx=find(SFCtimes>=starter2 & SFCtimes <=ender2);

%high pass filter pressure and temp
nhr=7/60;
srate=60;
w=1./(nhr.*srate);
wn=w./.5;
[b,c]=butter(5,wn,'high');
pprime=filtfilt(b,c,pres(tidx));

%low pass
% nhr=3/60;
% w=1./(nhr.*srate);
% wn=w./.5;
% [b2,c2]=butter(5,wn,'low');
% pprime=filtfilt(b2,c2,pprime);


figure(10);
subplot(2,1,2)
plot(pprime);
axis tight
