%script to plot final blended PCAPS time height data set
close all
clear all

scrsz=get(0,'screensize');

%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

%example of calculating Dtheta/Dz using matlab's build in gradient function
%(which uses centered differences on the interior and one sided differences
%on the edges

[dthetadz,dthetadt]=gradient(THblendRWP,vert,times);

%%
close all;

for IOP=[1 4 5 6 9]
    if IOP ==1
        starter=datenum(2010,11,29,00,00,00);
        ender=datenum(2010,12,08,0,00,00);
        %starter=datenum(2010,12,2,18,00,00);
        %ender=datenum(2010,12,04,00,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);
        %starter=datenum(2011,1,28,12,00,00);
        %ender=datenum(2011,2,1,00,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);
    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);
    elseif IOP==12
        %cold fropa
        starter=datenum(2010,12,20,00,00,00);
        ender=datenum(2010,12,21,00,00,00);        
    end
    midx=find(master>= starter & master<=ender);
    tidx1=find(times>= starter & times<=ender);
    tidx2=find(times2>= starter & times2<=ender);
    iopnum=strcat('0',num2str(IOP));


    f=figure(IOP);
    clf
    set(f,'position',[scrsz(1) scrsz(2) scrsz(3) scrsz(4)*.85],'color','white')
    subplot(2,1,1)
    %pcolor(master(midx),hgts,DPth(:,midx));shading flat;
    pcolor(times(tidx1),vert,(sqrt(u(tidx1,:).^2+v(tidx1,:).^2)')); shading flat; 
    %hold on;
    %contour(times(tidx1),vert,shear(tidx1,:)',.015:.005:.1,'-w','linewidth',2);
    %pcolor(times2(tidx2),vert2,tv(tidx2,:)'); shading flat
    %caxis([-5 5]);
    %caxis([-15 10]);
    caxis([0 15]);
    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])
    
    ylim([1310 3000])
    colorbar;
    subplot(2,1,2)
    pcolor(master(midx),hgts,RHth(:,midx));shading flat;
    caxis([0 100]);
    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])
    ylim([1310 3000])
    colorbar;
    
    
end
