%% 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,15,00,00);
        ender=datenum(2010,12,3,7,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=starter;%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=starter;%datenum(2010,12,02,14,0,0);
[~,rasstidx]=min(abs(timeRASS-starter2) );
[~,rassend]=min(abs(timeRASS-ender));%(starter2+(15/24))));
ender2=ender;%starter2+(18/24);
f=figure(10);clf; set(f,'color','w','units','pixels','pos',[1200 100 1600 900]);

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



%% RASS FROUDE
hidx=find(heightRASS==1972);
def=thetaRASS(hidx,rasstidx)-thetaRASS;def(def<0)=0;
cdef=9.81.*def./thetaRASS(hidx,rasstidx);
BH=trapz(cdef*60,1);
spd=sqrt(uRASS.^2 + vRASS.^2);
spd(def>1)=nan;
umax=max(spd,[],1);
%umax=nanmean(spd,1);
Fr=umax./sqrt(BH);

tspan=(rasstidx:rassend);




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(tmat(1,:),hmat(:,1)-1350,bssmooth);shading interp;
pcolor(timeRASS(tspan),heightRASS-1350,def(:,tspan));caxis([0 15]);shading interp;
%contourf(1:length(tidx2),height,bsnow(:,tidx2),[-10:.005:-6.5],'linestyle','none');
cbh=colorbar;
ylabel(cbh,'Temperature Deficit','fontsize',12,'fontweight','bold');
axis tight;
hold on;
contour(timeRASS(rasstidx:rassend),heightRASS-1350,thetaRASS(:,rasstidx:rassend),0:5:500,'k','linewidth',2);
contour(timeRASS(rasstidx:rassend),heightRASS-1350,thetaRASS(:,rasstidx:rassend),0:1:500,'k','linewidth',1);

ylim([1350 1900]-1350);
%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);



% close all
% figure;
% subplot(3,1,1:2);
% hold on;
% pcolor(timeRASS(tspan),heightRASS,def(:,tspan));caxis([0 20]);shading interp;
% contour(timeRASS(tspan),heightRASS,thetaRASS(:,tspan),0:5:500,'k','linewidth',2);
% contour(timeRASS(tspan),heightRASS,thetaRASS(:,tspan),0:1:500,'k','linewidth',1);
% datetick('x');
% axis tight
% 
% subplot(3,1,3);
% plot(timeRASS(tspan),Fr(tspan));
% axis tight;
% ylim([0 5]);
% set(gca,'ytick',[0:1:5]);
% datetick('x');
% axis tight
% 
% grid on;
%%
% [~,fidx1]=min(abs(Fr(tspan)-1.2));
% 
% [~,fidx2]=min(abs(Fr(tspan)-2.1));
% 
% bstvec=tmat(1,:);
% [~,tidx1]=min(abs(timeRASS(tspan(fidx1))-bstvec));
% [~,tidx2]=min(abs(timeRASS(tspan(fidx2))-bstvec));
% hh=hmat(:,1);
% f=figure(77);clf;set(f,'color','w','render','painters');
% subplot(2,1,1)
% pcolor((tidx1-200):(tidx1+400),hh(1:100),bssmooth(1:100,tidx1-200:tidx1+400));shading interp;
% caxis([-7.75 -6.3]);
% box on; 
% set(gca,'layer','top','linewidth',2,'fontsize',10,'fontweight','bold')
% xticker=datestr(bstvec((tidx1-200):100:tidx1+400),'HH:MM');
% set(gca,'xticklabel',{xticker});
% text(tidx1-195,1750,'a) Fr=1','fontsize',14,'fontweight','bold');
% ylabel('Height [m]','fontsize',10,'fontweight','bold')
% xlabel('Time [HH:MM]','fontsize',10,'fontweight','bold')
% 
% subplot(2,1,2)
% pcolor((tidx2-200):(tidx2+400),hh(1:100),bssmooth(1:100,tidx2-300:tidx2+300));shading interp;
% caxis([-7.75 -6.3]);
% box on; 
% xticker=datestr(bstvec((tidx2-200):100:tidx2+400),'HH:MM');
% set(gca,'xticklabel',{xticker});
% set(gca,'layer','top','linewidth',2,'fontsize',10,'fontweight','bold')
% ylabel('Height [m]','fontsize',10,'fontweight','bold')
% xlabel('Time [HH:MM]','fontsize',10,'fontweight','bold')
% text(tidx2-195,1750,'b) Fr=2','fontsize',14,'fontweight','bold');
% cbh=colorbar;
% cpos=get(cbh,'pos');
% cpos(1)=cpos(1)+.1;
% cpos(3)=.5.*cpos(3);
% set(cbh,'pos',cpos);
%% 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);

plot(timeRASS(tspan),Fr(tspan),'-r','linewidth',2);
axis tight;
spos4=get(sp4,'pos');
spos4(2)=spos3(2)-spos4(4);
set(sp4,'pos',spos4);
set(gca,'ytick',0:1:5);
%xticker=get(gca,'xtick');
%xlab=datestr(SFCtimes(tidx(xticker+1)),'HH:MM');
%set(gca,'xticklabel',{xlab});
datetick('x');
axis tight;
grid on; box on; 
ylim([0 5])

%plot([1 length(tidx)],[0 0],'--k','linewidth',1);

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');
%export_fig pcaps_iop1.png -m2
return
%% now compare with wrf data
addpath('/uufs/chpc.utah.edu/common/home/u0559900/RESEARCH/WRF/matlab_codes/post_process/');
load('wrf_sounding.mat','tint','uint','wint','tbar','ubar','wbar','tvec','zvec','Rtime');
    
%def
hidx=find(zvec==501);
def=-(tint-tint(1,hidx));
def(def<0)=0;
cdef=9.81.*def./(tint(1,hidx));
umax=max(uint,[],2);

BH=trapz(cdef*10,2);
Fr=umax./sqrt(BH);

%%
f=figure(7);clf; set(f,'color','w','units','pixels','pos',[1200 100 1600 900]);

    subplot(4,1,1:2)
    cmap=jet(24);
    colormap((cmap));
    pcolor(tvec/3600,zvec,def');shading flat; ch=colorbar;
    caxis([0 10]);
    colorbar;
    hold on;
    contour(Rtime/60,zvec,tbar',[0:5:500],'-k','linewidth',2);
    contour(Rtime/60,zvec,tbar',[0:1:500],'-k','linewidth',1);

    xlim([0 15]);
    box on;
    set(gca,'layer','top','linewidth',2)
    set(gca,'fontsize',12,'fontweight','bold','box','on','linewidth',2);
ylabel('Height [m]','fontsize',12,'fontweight','bold');

    sp3=subplot(4,1,3);
    spos3=get(sp3,'pos');
    spos3(2)=bot-spos3(4)-.001;
    set(sp3,'pos',spos3);
    hold on;
    plot(tvec/3600,tint(:,2),'r','linewidth',2);
    grid on;
    xlim([0 15]);
    ylim([263 275]);
        box on; set(gca,'layer','top')
    cbh=colorbar;
    set(cbh,'visible','off')
    set(gca,'fontsize',12,'fontweight','bold','box','on','linewidth',2);
ylabel('Pot. Temp. [K]','fontsize',10,'fontweight','bold');
    
    sp4=subplot(4,1,4); hold on;
    plot(tvec/3600,Fr,'r','linewidth',2);
    xlim([0 15]);
    %ylim([-5 10]);
    ylim([0 5]);
    box on; set(gca,'layer','top','ytick',0:1:5)
    spos4=get(sp4,'pos');
    spos4(2)=spos3(2)-spos4(4);
set(sp4,'pos',spos4);
cbh=colorbar;
grid on;
set(cbh,'visible','off')
    set(gca,'fontsize',12,'fontweight','bold','box','on','linewidth',2);
xlabel('Time [hrs]','fontsize',12,'fontweight','bold');
ylabel('Wind Speed [m/s]','fontsize',10,'fontweight','bold');

%export_fig wrf_pcaps_comp.png -m2