%% FROPA FIGURES

%% 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
addpath('/uufs/chpc.utah.edu/common/home/u0559900/RESEARCH/PCAPS_SURF/SFC_PRESSURE_GRADIENT/');

%% define times of interest

for IOP=[1]%3 4 5 6 9]
    if IOP ==1
        starter=datenum(2010,12,2,18,00,00);
        ender=datenum(2010,12,3,18,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>1300);
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);

%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);

%% Overview
starter2=datenum(2010,12,3,2,0,0);
ender2=starter2+(8/24);
f=figure(10);clf; set(f,'color','w','units','pixels','pos',[1200 100 1600 900]);
%tsp=tight_subplot(4,1);
%axes(tsp(1:2));
subplot(4,1,1:2)
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 interp;
%pcolor(1:length(bssmooth),height,bssmooth);shading interp;
%contourf(1:length(tidx2),height,bsnow(:,tidx2),[-10:.005:-6.5],'linestyle','none');
cbh=colorbar;
ylabel(cbh,'Backscatter','fontsize',12,'fontweight','bold');
axis tight;
ylim([1350 1700]);
xticker=get(gca,'xtick');
xlab=tmat(1,xticker);
set(gca,'xtick',[],'xticklabel',[])
ylabel('Height [m]','fontsize',12,'fontweight','bold');
set(gca,'fontsize',10,'fontweight','bold','box','on','layer','top','linewidth',2);
spos=get(gca,'pos');
bot=spos(2);
text(30,1685,'a) Ceilometer','backgroundcolor',[.7 .7 .7],'edgecolor','k','fontsize',12,'fontweight','bold');
%% load ISS surface data

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
ISStemp=data(:,7)+273.15;
%temp=interp1(tvec,temp,SFCtimes);
ISSpres=data(:,11);
%pres=interp1(tvec,pres,SFCtimes);

ISStheta=ISStemp.*(1000./ISSpres).^(2/7);
%%

ISSu = - data(:,end-6) .* sin ( (pi/180) .* data(:,end-5) );
ISSv = - 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);
sp3=subplot(4,1,3);
spos3=get(sp3,'pos');
spos3(2)=bot-spos3(4)-.001;
set(sp3,'pos',spos3);
plot(1:length(tidx),ISStheta(tidx),'r','linewidth',2);axis tight;
hold on;
grid on; box on;
set(gca,'yaxislocation','right','fontsize',12,'fontweight','bold');
ylabel('[K]','fontsize',12,'fontweight','bold');
set(gca,'xticklabel',[]);
set(gca,'fontsize',12,'fontweight','bold','box','on','linewidth',2);
cbh=colorbar; set(cbh,'visible','off');
text(5,292,'b) Surface Pot. Temp.','backgroundcolor',[.7 .7 .7],'edgecolor','k','fontsize',12,'fontweight','bold');

%% 
sp4=subplot(4,1,4); hold on;
plot(1:length(tidx),ISSv(tidx),'-r','linewidth',2);
axis tight;
spos4=get(sp4,'pos');
spos4(2)=spos3(2)-spos4(4);
set(sp4,'pos',spos4);
xticker=get(gca,'xtick');
xlab=datestr(SFCtimes(tidx(xticker+1)),'HH:MM');
set(gca,'xticklabel',{xlab});
grid on; box on; 
plot([1 length(tidx)],[0 0],'--k','linewidth',1);
grid on;
box on;
set(gca,'fontsize',12,'fontweight','bold');
xlabel('Time [UTC]','fontsize',10,'fontweight','bold');
ylabel('[m/s]','fontsize',10,'fontweight','bold');
set(gca,'fontsize',12,'fontweight','bold','box','on','linewidth',2);
cbh=colorbar; set(cbh,'visible','off');
text(5,7,'c) Meridional Wind','backgroundcolor',[.7 .7 .7],'edgecolor','k','fontsize',12,'fontweight','bold');

%% Add ISS3 to the figure
fname='ISS3.dat';
data=importdata(fname,',',6);
hdr=data.textdata;
hdr=(hdr(1)); 

ST=regexp(hdr,' ','split');
ST=ST{1};
elev=str2double(ST(5));
lat3=str2double(ST(3));
lon3=str2double(ST(4));

data=data.data;
yr=data(:,3);
mn=data(:,1);
dy=data(:,2);
hr=data(:,4);
mins=data(:,5);
times=datenum(yr,mn,dy,hr,mins, zeros(length(mins),1));
datestr(times(1:5))
temp3=data(:,7)+273.15;
pres3=data(:,11);
spd3=data(:,9);
drct3=data(:,10);
u3 = - spd3 .* sin ( (pi/180) .* drct3);
v3 = - spd3 .* cos ( (pi/180) .* drct3 );

iss3tidx=find(times>=starter2 & times<= ender2);
theta3=temp3.*(1000./pres3).^(2/7);

tnow=interp1(times(iss3tidx),theta3(iss3tidx),SFCtimes(tidx));
vnow=interp1(times(iss3tidx),v3(iss3tidx),SFCtimes(tidx));
axes(sp3)
plot( 1:length(tidx),tnow,'-k','linewidth',2);axis tight;

legend('ISS-S','ISS-N')
axes(sp4)
hold on;
plot(1:length(tidx),vnow,'-k','linewidth',2);

return;
%% Timing of warm and cold fronts
figure(51);clf;
plot(SFCtimes(tidx),ISStheta(tidx),'linewidth',2);
hold on;
plot(SFCtimes(tidx),tnow,'k','linewidth',2);
xlim([datenum(2010,12,3,5,0,0) datenum(2010,12,3,6,45,0)]);
grid on; box on;
datetick('x','keeplimits')

%distance between sites
between=deg2km(lat3-lat);

figure(52);clf;
plot(SFCtimes(tidx),ISStheta(tidx),'linewidth',2);
hold on;
plot(SFCtimes(tidx),tnow,'k','linewidth',2);
xlim([datenum(2010,12,3,7,0,0) datenum(2010,12,3,7,45,0)]);
grid on; box on;
datetick('x','keeplimits');

%% CAP structure between the fronts

[~,rassidx]=min(abs(timeRASS- datenum(2010,12,3,6,40,0)));
figure(53);clf
plot(thetaRASS(:,[(rassidx):(rassidx+1)]),heightRASS);
ylim([1300 1800]);
%legend(datestr(timeRASS(rassidx)))
legend

%% Density Current Cold FROPA

[~,rassidx]=min(abs(timeRASS- datenum(2010,12,3,7,05,0)));
[~,rassidx2]=min(abs(timeRASS- datenum(2010,12,3,7,35,0)));

%[left bottom width height]


starter2=datenum(2010,12,03,7,03,0);
ender2=starter2+(.8/24);
f=figure(100);clf; set(f,'color','w','renderer','painters');%,'units','pixels','pos',[100 100 1000 200]);

sp1=subplot(1,8,1);
spos1=get(sp1,'pos');
spos1(1)=spos(1)-.05;
spos1(3)=1.4*spos1(3);
set(sp1,'pos',spos1);
spos1=get(sp1,'pos');
redge=spos1(1)+spos1(3);%left + width;
plot(thetaRASS(:,rassidx), heightRASS,'-k','linewidth',2);
%plot(vRASS(:,rassidx),heightRASS,'-k','linewidth',2);
grid on;
ylim([1310 1700]);
set(gca,'fontsize',10,'fontweight','bold','box','on','layer','top','linewidth',1);
xlim([291 301])
ylabel('Height [m]','fontsize',10,'fontweight','bold');
xlabel('Pot. Temp. [K]','fontsize',10,'fontweight','bold');
text(291.5,1670,'a)Pre-frontal: 0705 UTC','backgroundcolor',[.7 .7 .7],'edgecolor','k','fontsize',10,'fontweight','bold');

sp2=subplot(1,8,2:7);
spos2=get(sp2,'pos');
spos2(1)=redge;
set(sp2,'pos',spos2);
cmap=[gray(12);flipud(gray(12))];
cmap(1:12,3)=1;
cmap(13:end,1)=1;
colormap(cmap);
colormap((jet(24)))
tidx2=find(time>=starter2 & time<=ender2);
caxis([-7.75 -6.8]);
hold on
bssmooth=interp2(bsnow(:,tidx2),2);
[tmat,hmat]=meshgrid(time(tidx2),height);
tmat=interp2(tmat,2);
hmat=interp2(hmat,2);
%bb=ones(2,2)/4;
%bssmooth(hmat<1330)=nan;%filter2(bb,bssmooth(hmat<1350));%1:size(bssmooth,2)
pcolor(1:size(bssmooth,2),(hmat(:,1))',bssmooth);shading flat;
hold on;
vnow=vRASS(:,rassidx);
hnow=interp(heightRASS,4);
vnow=interp(vnow,4);
vnow2=vRASS(:,rassidx2);
vnow2=interp(vnow2,4);

quiver(100,hnow,vnow,0*vnow,.43,'k','linewidth',2)
quiver(600,hnow,vnow2,0*vnow2,.43,'k','linewidth',2)
axis tight;
%datetick('x','HH:MM:SS','keeplimits');
cbh=colorbar('north'); cpos=get(cbh,'pos'); cpos(3)=.3* cpos(3); cpos(1)=cpos(1)+.2; set(cbh,'pos',cpos);
set(cbh,'fontsize',10,'fontweight','bold','box','on','layer','top','linewidth',1);

ylim([1310 1700]);
set(gca,'yticklabel',[]);
spos2=get(sp2,'pos');
redge=spos2(1)+spos2(3);
xticker=get(gca,'xtick');
xlab=datestr(tmat(1,xticker),'HH:MM');%time(tidx2(xticker));
set(gca,'xticklabel',{xlab});
%set(gca,'xtick',[],'xticklabel',[])
%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);
text(10,1670,'b) Gravity Current','backgroundcolor',[.7 .7 .7],'edgecolor','k','fontsize',10,'fontweight','bold');

sp3=subplot(1,8,8);
spos3=get(sp3,'pos');
spos3(1)=redge;
spos3(3)=1.4*spos3(3);
set(sp3,'pos',spos3);
plot(thetaRASS(:,rassidx2), heightRASS,'-k','linewidth',2);
%plot(vRASS(:,rassidx2), heightRASS,'-k','linewidth',2);

ylim([1310 1700]);
grid on;
set(gca,'fontsize',10,'fontweight','bold','box','on','layer','top','linewidth',1);
set(gca,'yaxislocation','right')
ylabel('Height [m]','fontsize',10,'fontweight','bold');
xlim([285 302])

text(285.5,1670,'c) Post-frontal: 0740 UTC','backgroundcolor',[.7 .7 .7],'edgecolor','k','fontsize',10,'fontweight','bold');
xlabel('Pot. Temp. [K]','fontsize',10,'fontweight','bold');

return


%% Warm front sounding ISS 06 UTC

file='/uufs/chpc.utah.edu/common/home/horel-group/pcaps.utah.edu/data/SONDE_PROC/20101203/NCAR_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);

%%
f=figure(81); clf; plot(theta,hgts,'linewidth',2)
tidxnew=find(SFCtimes == datenum(2010,12,3,5,14,0));
hold on;
plot(ISStheta(tidxnew),1320,'*r','linewidth',2);
plot(293,1320,'*r','linewidth',2)
ylim([1310 1800]);
grid on;

[~,rassidx]=min(abs(timeRASS- datenum(2010,12,3,5,14,0)));
hold on;
plot(thetaRASS(:,[(rassidx-1):(rassidx+1)]),heightRASS);