% co2 data

close all
clear all

scrsz=get(0,'screensize');


%PULL in CO2 data
addpath('/Users/neil/Documents/PCAPS_DATA/co2_data');

% load('network_pcaps_utc_short.mat');
% co2(co2<0)=nan;
% co2time=datenum(co2_time);
%
%
% load('res_output_fixwind.mat');
% outime=(out.gmt);
% outime=[outime,zeros(length(outime),2)];
% outime=datenum(outime);
% co2out=out.obs_sites;
% rosepark=co2out(:,7);

load('network_pcaps_newamb');
co2data=network;
outime=datenum(co2data(:,1),co2data(:,2),co2data(:,3),co2data(:,4),zeros(length(co2data),1),zeros(length(co2data),1));
co2out=co2data(:,5:end);


clrs={'b','k','g','c','y','m'};
%names={'Murray','Snowbird','Sugarhouse','University','Downtown','Kennecott','Rose Park'};
names={'SB','SH','UN','KT','RP','PM25'};


%% now try loading in hourly PM2.5 data
addpath('/Users/neil/Documents/PCAPS_DATA');

pm=importdata('PMhourly.dat',' ');


pmhr=pm(:,2:end-3);
pmhr(pmhr<0)=nan;
pmhr(pmhr>120)=nan;
pmseries=[];
for i=1:size(pmhr,1) %loop over each day
    pmseries=cat(2,pmseries,pmhr(i,:));
    
end

%creat PM timeseries then shift to UTC

pmtime=datenum(2010,11,1):(1/24):datenum(2011,2,28,23,00,00);

%now the times are off by 6 hours... (UTC = MST+7);

pmtime=pmtime+(7/24);

%now develop standardized anoms of pmseries
pmanom=(pmseries-nanmean(pmseries))./nanstd(pmseries);

%%


%first load the 1 hourly data (sounding grid)
load('double_blend.mat','blended','DPth','RHth','hgts','master');
%blended = Potential Temperature Time Height on Sounding Grid (uses RASS,
%surface obs, and soundings).
% DPth = Dew points time height
% RHth= Relative humidity time heights
% master = matlab serial time vector
% hgts= heights;


%next load the RASS and RWP gridded data
load('RASS_RWP.mat','THblend','THblendRWP','times2','vert2','tv','u','v','shear','RI','dthdz','vert','times');
% times2 = times for RASS data grids
% vert2= heights for RASS grid
% times= times for RWP grids
% vert= range gates for RWP
% THblend = Potential temperature on RASS grid
% tv= RASS virtual temperature on RASS grid
% THblendRWP= Potential temperature on RWP grid
% dthdz= vertical gradient in potential temperature (on RWP grid);
% RI= richardson number (RWP grid)
% u= uwind on RWP grid
% v= vwind on RWP grid
% shear= shear on RWP grid


for IOP=[1 3 4 5 6 9];
    
    for it=1:3
        if IOP ==1
            starter=datenum(2010,12,01,00,00,00);
            ender=datenum(2010,12,08,0,00,00);
            idx=find(times>= starter & times<=ender);
            idx2=find(times2>=starter & times2<=ender);
            iopnum=strcat('0',num2str(IOP));
            
        elseif IOP==2
            starter=datenum(2010,12,7,00,00,00);
            ender=datenum(2010,12,11,0,00,00);
            idx=find(times>= starter & times<=ender);
            idx2=find(times2>=starter & times2<=ender);
            iopnum=strcat('0',num2str(IOP));
        elseif IOP==3
            starter=datenum(2010,12,12,00,00,00);
            ender=datenum(2010,12,15,0,00,00);
            idx=find(times>= starter & times<=ender);
            idx2=find(times2>=starter & times2<=ender);
            iopnum=strcat('0',num2str(IOP));
            
        elseif IOP==4
            starter=datenum(2010,12,23,18,00,00);
            ender=datenum(2010,12,26,19,00,00);
            idx=find(times>= starter & times<=ender);
            idx2=find(times2>=starter & times2<=ender);
            iopnum=strcat('0',num2str(IOP));
            tmin=-5;
            tmax=9;
            wmin=0;
            wmax=15;
            
        elseif IOP==5
            starter=datenum(2010,12,31,20,00,00);
            ender=datenum(2011,1,10,0,00,00);
            idx=find(times>= starter & times<=ender);
            idx2=find(times2>=starter & times2<=ender);
            iopnum=strcat('0',num2str(IOP));
            
        elseif IOP==6
            starter=datenum(2011,1,11,00,00,00);
            ender=datenum(2011,1,18,0,00,00);
            idx2=find(times2>=starter & times2<=ender);
            idx=find(times>= starter & times<=ender);
            iopnum=strcat('0',num2str(IOP));
            
        elseif IOP==7
            starter=datenum(2011,1,19,00,00,00);
            ender=datenum(2011,1,23,0,00,00);
            idx2=find(times2>=starter & times2<=ender);
            idx=find(times>= starter & times<=ender);
            iopnum=strcat('0',num2str(IOP));
            
        elseif IOP==8
            starter=datenum(2011,1,22,00,00,00);
            ender=datenum(2011,1,27,0,00,00);
            idx=find(times>= starter & times<=ender);
            idx2=find(times2>=starter & times2<=ender);
            iopnum=strcat('0',num2str(IOP));
            
        elseif IOP==9
            starter=datenum(2011,1,26,00,00,00);
            %starter=datenum(2011,1,28,12,00,00);
            %ender=datenum(2011,2,1,00,00,00);
            ender=datenum(2011,1,31,08,00,00);
            idx=find(times>= starter & times<=ender);
            idx2=find(times2>=starter & times2<=ender);
            iopnum=strcat('0',num2str(IOP));
            
        elseif IOP==10
            starter=datenum(2011,2,1,00,00,00);
            ender=datenum(2011,2,6,0,00,00);
            idx=find(times>= starter & times<=ender);
            idx2=find(times2>=starter & times2<=ender);
            iopnum=strcat(num2str(IOP));
            
        elseif IOP==11 %not really an IOP just the mix out during IOP1
            starter=datenum(2010,12,02,15,00,00);
            ender=datenum(2010,12,03,18,00,00);
            idx=find(times>= starter & times<=ender);
            idx2=find(times2>=starter & times2<=ender);
            iopnum='01';
            tmin=-5;
            tmax=18;
            wmin=0;
            wmax=18;
            
        elseif IOP==12
            %cold fropa
            starter=datenum(2010,12,20,00,00,00);
            ender=datenum(2010,12,21,00,00,00);
            idx=find(times>= starter & times<=ender);
            idx2=find(times2>=starter & times2<=ender);
            
        end
        
        [t2mat,v2mat]=meshgrid(times2(idx2),vert2);
        t2mat=interp2(t2mat,1);
        v2mat=interp2(v2mat,1);
        THnow=THblend(idx2,:)';
        THnow=interp2(THnow,1);
        [tmat,hgtmat]=meshgrid(times(idx),vert);
        tmat=interp2(tmat,1);
        hgtmat=interp2(hgtmat,1);
        dtdz=dthdz(idx,:)';
        dtdz=interp2(dtdz,1);
        
        cfile=strcat('/Users/neil/Documents/PCAPS_DATA/ceil_carolyn_iop0',num2str(IOP),'.DAT');
        ceilo=load(cfile);
        agl=ceilo(1,:);agl(1)=[];
        ceilo(1,:)=[];
        tstamp=ceilo(:,1);
        ceilo(:,1)=[];
        time2=double(tstamp)./86400+datenum(1970,01,01,00,00,00); %convert to matlab serial time from unix serial time
        
        
        if it==1
            f=figure(200+IOP);
            clf
            set(f,'position',[scrsz(1) scrsz(2) .9*scrsz(3) .9*scrsz(4)],'color','white')
            subplot(8,1,1:3);
            pcolor(time2,agl+1320',ceilo');shading flat;caxis([-8 -5]); colorbar;
            hold on;
            midx=find(master>= starter & master<=ender);
            
            %contour(t2mat,v2mat,THnow,200:1:400,'-','linewidth',1,'linecolor',[.7 .7 .7])
            %contour(t2mat,v2mat,THnow,200:5:400,'-','linewidth',2,'linecolor',[.5 .5 .5])
            contour(master(midx),hgts,blended(:,midx),200:1:400,'-','linewidth',1,'linecolor',[.7 .7 .7])
            contour(master(midx),hgts,blended(:,midx),200:5:400,'-','linewidth',2,'linecolor',[.5 .5 .5])
            
            datetick('x','dd-HH');
            set(gca,'xticklabel',[]);
            axis tight
            title('Ceilometer Analysis');
            ylim([1310 2800]);
        elseif it==2
            figure(200+IOP)
            subplot(8,1,4:6)
            midx=find(master>= starter & master<=ender);
            %pcolor(master(midx),hgts,RHth(:,midx));shading flat;
            pcolor(tmat,hgtmat,dtdz);shading flat;
            %caxis([20 100]);
            caxis([0 .055]);
            hold on;
            datetick('x','dd-HH')
            axis tight
            hold on;
            contour(master(midx),hgts,blended(:,midx),200:1:400,'-','linewidth',1,'linecolor',[.7 .7 .7])
            contour(master(midx),hgts,blended(:,midx),200:5:400,'-','linewidth',2,'linecolor',[.5 .5 .5])
            colorbar;
            %xlabel('Potential Temperature Lapse Rate (shaded)')
            datetick('x','dd-HH');
            axis tight
            ylim([1310 2800]);
            title('STABILITY: theta lapse rate');
            
        elseif it==3
            figure(200+IOP);
            subplot(8,1,7:8);
            cstarter=starter;
            cender=ender;
            cidx=find(outime>=cstarter & outime<=cender);
            a=0;
            for i=[2 3 4 6 7]%1:size(co2out,2);
                a=a+1;
                hold on;
                if i<7
                    plot(outime(cidx),(co2out(cidx,i)-nanmean(co2out(:,i)))./nanstd(co2out(:,i)),'color',clrs{a},'linewidth',2);
                    datetick('x','dd-HH');
                    grid on;
                else
                    clear ax h1 h2
                    pmidx=find(pmtime>=cstarter & pmtime<=cender);
                    
                    [ax,h1,h2]=plotyy(outime(cidx),(co2out(cidx,i)-nanmean(co2out(:,i)))./nanstd(co2out(:,i)),pmtime(pmidx),pmseries(pmidx));
                    set(h1,'color',clrs{a},'linewidth',2);
                    set(h2,'linestyle','--','color','w','linewidth',2);
                    datetick('x','dd-HH');
                    set(ax(2),'xtick',[],'xticklabel',[]);
                    set(ax,'YColor','k','Xcolor','k')
                    set(ax(1),'ytick',[-1 0 1 2 3 4 5 6]);
                    set(ax(2),'ticklength',[0 0]);
                    set(get(ax(2),'Ylabel'),'String','PM2.5') 
                    set(gca,'color',[.9 .9 .9])
                    grid on;
                end
            end
            %pmidx=find(pmtime>=cstarter & pmtime<=cender);
            %plot(pmtime(pmidx),pmseries(pmidx),'color',clrs{end},'linewidth',2);
            cbh=colorbar;
            set(cbh,'visible','off');
            lh=legend(names);
            lpos=get(lh,'pos');
            lpos(1)=lpos(1)+.11;
            lpos(3)=.6*(lpos(3));
            set(lh,'pos',lpos);
            axis tight
            ylabel('Std. Anom');
            xlabel('CO2 Standardized Anomaly Time Series');
        end
        if min(time2)<starter
            xlow=starter;
        elseif min(time2)>=starter
            xlow=min(time2);
        end
        if max(time2)>ender
            xhigh=ender;
        elseif max(time2)<=ender
            xhigh=max(time2);
        end
        %xlim([starter xhigh]);
        xlim([starter ender])
        
        if it==3
            export_fig(sprintf('IOP%d_CO2.png', IOP),'-m2.5')
        end
        
    end
end





