%script to plot final blended PCAPS time height data set
close all
clear all
opengl
scrsz=get(0,'screensize');

% read pm2.5 hourly data. from erik crosman
% mainly interested in Hawthorne
PM = csvread('PM25.txt');
lpm = length(PM);
ptim = datenum(PM(:,3),PM(:,1),PM(:,2),PM(:,4),00,00);
ptim = ptim+7/24; % in mst
haw = PM(:,8);

%load the 1 hourly data (sounding grid). from neil lareau and preprocessing
%to compute rmin values

load('rmin_300.mat','blended','BlendTC','PRES','DPth','RHth','hgts','master','ublend','vblend','rmin');
%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;

% four tunable parameters

% depth over which rmin evaluated (here 300 m- see file name)
% some dependency on the depth, but if larger depth (say 500), tends to be
% choppier

%no_ext- number of extrema to shoot for- fewer means going to be really
%smooth
no_ext = 10; % can have max of 4-5 stable and mixed layers
% to be a stable layer must have at least rminthres value of rmin
% some dependency but if rminthres is more negative, just makes more mixed
% layers
rminthres = -2.5;
% stablethres is value of min rmin for strong stable layer
stablethres = -10;


layer_depth = 10; % 100 m thick layers
top = 271;
bot = 3; % skip over bottom two layers that seem too warm much of the time
tlev = top-bot+1;


%load ceilometer data
load('horel_ceil_pcaps.mat','bs','range','time');
time2=double(time)./86400+datenum(1970,01,01,00,00,00); %convert to matlab serial time from unix serial time
msl = range + 1320;
%rmino = rmin;
no_time = length(master);

%smoothing rmin values really doesn't make a big difference so not used
%after testing
%for i=1:no_time % lightly smooth in vertical
%rmins1 = median_smooth(rmin(bot:top,i),tlev,5);
%rmins(bot:top,i) = median_smooth(rmins1,tlev,5);
%end
%for j=bot:top % lightly smooth in time
%rmins2 = median_smooth(rmins(j,:),no_time,5);
%rmins3(j,:) = median_smooth(rmins2,no_time,5);
%end
%rmin(bot:top,:) = rmins3(bot:top,:);

%next load the RASS and RWP gridded data. from neil lareau
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

s = zeros(top,no_time);
st = zeros(top,no_time);

for i=1:no_time

           if (rmin(:,i) > rminthres)  
               % then there is only 1 strong mixed layer and no stable layers
                layers(1,1,i) = hgts(1,bot);
                layers(2,1,i) = hgts(1,top);
                layers(3,1,i) = max(rmin(bot:top));
                layers(4,1,i) = 1;
                layer(i) = 1;
                s(bot:top,i) = 1;
           else
          clear hxs 
          [eind,eval] = extrema_rmin_final(rmin(bot:top,i),no_ext); % find extrema within profile 
          ntg = length(eind);
          hxs(:,1)=eind;
          hxs(:,2)=eval;
          
        % figure(100)
         
         % hold on
        %  plot(rmin(bot:top,i),hgts(1,bot:top),'r');
        %  plot(hxs(:,2),hgts(1,hxs(:,1)),'b')
        %  hold off
       %   pause

       layer(i) = ntg-1;
           for h=1:ntg-1
               if(hxs(h+1,2)<hxs(h,2)) % mixed layer
                 s(hxs(h,1):hxs(h+1,1),i)=1;
               else
                    %  s(hxs(h,1):hxs(h+1,1),i)=2;
                   for hh=hxs(h,1):hxs(h+1,1)
              if(rmin(hh,i)<=stablethres) % strongly stable 
                  s(hh,i)=3;
              elseif(rmin(hh,i)>stablethres&&rmin(hh,i)<rminthres)
                  s(hh,i)=2;
              else
                  s(hh,i)=1;
              end
                   end
                end
                
           
           end
           end
end
%%

% smooth in time lightly

st = zeros(top,no_time);
st2 = zeros(top,no_time);
st3 = zeros(top,no_time);
for j=bot:top % lightly smooth in time
st(j,:) = median_smooth(s(j,:),no_time,3);
st2(j,:) = median_smooth(st(j,:),no_time,3);
st3(j,:) = median_smooth(st2(j,:),no_time,3);
end

%inpaint smoothing option. not used
%df = zeros(top,no_time);
%st = zeros(top,no_time);
%df = zeros(top,no_time);
%st2 = zeros(top,no_time);
%st3 = zeros(top,no_time);
%df = diff(s,1,2);
%df(:,no_time)= df(:,no_time-1);
%st2=s;
%st2(abs(df)>=1)=nan;
%st3 = inpaint_nans(st2,1);
%st3 = round(st3);

%extra smoothing
%df = diff(st3,1,1);
%df(top,:)= df(top-1,:);
%st2=st3;
%st2(abs(df)>=1)=nan;
%st3 = inpaint_nans(st2,1);
%st3 = round(st3);
%st(st3<5) = 1;
%st(st3>=5)  = 2;

% reloading after smoothing and keeping original for comparison
st(st3==1) = 1;
st(st3==2)  = 2;
st(st3==3)  = 3;
stable = zeros(top,no_time);
stable(s==1) = 1;
stable(s==2) = 2;
stable(s==3)  = 3;

% figure out mixed layer depth
mix_depth = nan(no_time,2);
imix=0;
tmix = zeros(no_time,1); % traditional mix layer depth
tmixadd = 1.0; % amount of extra buoyancy
tmixh = 100; % mix layer depth if stable

for i=1:no_time
    %trad mix layer depth. find lowest value of theta greater than surface
    %value + little
    imt = find(blended(:,i)>blended(bot,i)+tmixadd); % finding all the values where 
    % theta greater than bot val + tmixadd
    limt = length(imt);
    if(limt==0) % really unstable through whole depth
        tmix(i)=top;
    else
        tmix(i)=hgts(1,imt(1)+1); % lowest one
        if(tmix(i)<hgts(1,bot)+tmixh)
            tmix(i)=hgts(1,bot)+tmixh; % minimum depth is tmixh
        end
    end
    
    % using smoothed in time estimates of stable layers
    dfm = diff(st(bot:top,i),1,1);
    indd = find(abs(dfm)>0);
    lin = length(indd);
    mid1 = tmixh + 2*(rmin(bot,i)+25); % creating a little variability when stable
    mid1 = min([mid1 200]);
    mid  = hgts(1,bot)+max([mid1 50]);
    
    if(lin == 0) % only one layer
        if(st(bot,i)==1) % if surface based mixed layers
               mix_depth(i,1)=hgts(1,top);
               imix = 1;
        else
           mix_depth(i,1) = mid;
                   imix=0;
        end 
        mix_depth(i,1) = max([mix_depth(i,1) tmix(i)]); % go with traditional mixing depth
        % if mix_depth lower
    elseif (lin==1) % handling separately because checking otherwise for height of layer 0
    dep1 = hgts(1,indd(1)+1);
       if(st(bot,i)==1) % if surface based mixed layers
               mix_depth(i,1)=max([dep1 mid]);
               imix = 1;
       else %surface stable layer
                     mix_depth(i,1) = mid;
                   imix=0;
       end
         mix_depth(i,1) = max([mix_depth(i,1) tmix(i)]); % go with trad mixing depth
        % if mix_depth lower
    else
         dep1 = hgts(1,indd(1)+1);
         dep2 = hgts(1,indd(2)+1);
           if(st(bot,i)==1) % if surface based mixed layers
             mix_depth(i,1)=max([dep1 mid]);
               imix = 1;
           else %surface stable layer
               mix_depth(i,1)=max([dep2 mid]); % height of next mixed layer
               if(dep1<=mix_depth(i-1,1)&&imix==1) % is top of stable layer below
                   %prior mixed layer?
                   %elevated mixed layer is shallow and overlaps with
                   %previous one
                   mix_depth(i,2)=mid; % keep lower one too
                   mix_depth(i,1)=max([dep2 mid]);
               else
                   mix_depth(i,1)=mid;
                   imix=0;
               end
           
           end
             mix_depth(i,1) = max([mix_depth(i,1) tmix(i)]); % go with trad mixing depth
        % if mix_depth lower
    end
end
 
%cosmetics for coloring plots
st(st==1) = 4;
st(st==2)  = 3.25;
st(st==3)  = 2.5;
stable(s==1) = 4;
stable(s==2) = 3.25;
stable(s==3) = 2.5;


%for IOP=[ 01  ]
  %  for IOP=[ 01 02 03 04 05 06 07 08 09 13 14 ]
          for IOP=[ 01 03 04 05 06 07 08 09 ]
    IOPstr = num2str(IOP);
    if IOP ==1
        starter=datenum(2010,11,29,00,00,00);
        ender=datenum(2010,12,11,0,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);     
         elseif IOP==13
        starter=datenum(2010,12,01,00,00,00);
        ender=datenum(2010,12,31,23,00,00); 
          elseif IOP==14
        starter=datenum(2011,1,01,00,00,00);
        ender=datenum(2011,12,31,23,00,00);  
    end
    midx=find(master>= starter & master<=ender);
    midc=find(time2>= starter & time2<=ender);
    tidx2=find(times2>= starter & times2<=ender);
    tidx1=find(times>= starter & times<=ender);
    tidx2=find(times2>= starter & times2<=ender);
      midp=find(ptim>= starter & ptim<=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(7,1,1:3)
    colorbar('off')
     pcolor(master(midx),hgts(1,1:271),stable(:,midx));shading flat;
         colorbar('off')
               colormap(pink) 
    caxis([0 4]);
    hold on;
    datetick('x','dd')
    axis tight
    hold on;
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:1:400,'-','linewidth',1,'linecolor',[.7 .7 .7])
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:5:400,'-','linewidth',2,'linecolor',[.5 .5 .5])
    p=plot(master(midx),tmix(midx),'o');
    set(p,'MarkerFaceColor',[1 .7 0],'MarkerEdgeColor',[1 .7 0],'MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,1),'o');
    set(p,'MarkerFaceColor','k','MarkerEdgeColor','k','MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,2),'o');
    set(p,'MarkerFaceColor','y','MarkerEdgeColor','y','MarkerSize',3);
    ylim([1320 3000])
       freezeColors
       subplot(7,1,4:6)
    
    %pcolor(master(midx),hgts(1,1:271),RHth(1:271,midx));shading flat;
    pcolor(time2(midc),msl,bs(midc,:)');shading flat;caxis([-8 -5]);
    colorbar('off');
    
    colormap(jet) 
    %caxis([0 100]);
    
    hold on;
    datetick('x','dd')
    axis tight
    hold on;
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:1:400,'-','linewidth',1,'linecolor',[.7 .7 .7])
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:5:400,'-','linewidth',2,'linecolor',[.5 .5 .5])
     p=plot(master(midx),tmix(midx),'o');
    set(p,'MarkerFaceColor',[1 .7 0],'MarkerEdgeColor',[1 .7 0],'MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,1),'o');
    set(p,'MarkerFaceColor','k','MarkerEdgeColor','k','MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,2),'o');
    set(p,'MarkerFaceColor','y','MarkerEdgeColor','y','MarkerSize',3);
    ylim([1320 3000])
    freezeColors
    %cbfreeze(colorbar);
     
     subplot(7,1,7)     
     bar(ptim(midp),haw(midp),'k');  
     hold on;
     datetick('x','dd')
     axis tight
     ylim([0 100])
     
     print(gcf,'-dpng','-r600',strcat('ISS_IOP_',IOPstr,'_ST','.png'))
     
     f=figure(IOP+100);
    clf
    set(f,'position',[scrsz(1) scrsz(2) scrsz(3) scrsz(4)*.85],'color','white')
    subplot(2,1,1)
     pcolor(master(midx),hgts(1,1:271),st(:,midx));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 4]);
    colormap(pink) 
    hold on;
    datetick('x','dd')
    axis tight
    hold on;
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:1:400,'-','linewidth',1,'linecolor',[.7 .7 .7])
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:5:400,'-','linewidth',2,'linecolor',[.5 .5 .5])
     p=plot(master(midx),tmix(midx),'o');
    set(p,'MarkerFaceColor',[1 .7 0],'MarkerEdgeColor',[1 .7 0],'MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,1),'o');
    set(p,'MarkerFaceColor','k','MarkerEdgeColor','k','MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,2),'o');
    set(p,'MarkerFaceColor','y','MarkerEdgeColor','y','MarkerSize',3);
    
    ylim([1320 3000])
    freezeColors
    cbfreeze;
    subplot(2,1,2)
    %pcolor(master(midx),hgts(1,1:271),RHth(1:271,midx));shading flat;
    colormap(jet) 
    %caxis([0 100]);
    pcolor(times(tidx1),vert,(sqrt(u(tidx1,:).^2+v(tidx1,:).^2)')); shading flat; 
    caxis([0 15]);
    hold on;
    datetick('x','dd')
    axis tight
    hold on;
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:1:400,'-','linewidth',1,'linecolor',[.7 .7 .7])
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:5:400,'-','linewidth',2,'linecolor',[.5 .5 .5])
     p=plot(master(midx),tmix(midx),'o');
    set(p,'MarkerFaceColor',[1 .7 0],'MarkerEdgeColor',[1 .7 0],'MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,1),'o');
    set(p,'MarkerFaceColor','k','MarkerEdgeColor','k','MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,2),'o');
    set(p,'MarkerFaceColor','y','MarkerEdgeColor','y','MarkerSize',3);

    ylim([1320 3000])
    freezeColors
    cbfreeze(colorbar);
      print(gcf,'-dpng','-r600',strcat('ISS_IOP_',IOPstr,'_SP','.png'))
      
    f=figure(IOP+200);
    clf
    set(f,'position',[scrsz(1) scrsz(2) scrsz(3) scrsz(4)*.85],'color','white')
    subplot(7,1,1:3)
    colorbar('off')
     pcolor(master(midx),hgts(1,1:271),st(:,midx));shading flat;
         colorbar('off')
               colormap(pink) 
    caxis([0 4]);
    hold on;
    datetick('x','dd')
    axis tight
    hold on;
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:1:400,'-','linewidth',1,'linecolor',[.7 .7 .7])
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:5:400,'-','linewidth',2,'linecolor',[.5 .5 .5])
    p=plot(master(midx),tmix(midx),'o');
    set(p,'MarkerFaceColor',[1 .7 0],'MarkerEdgeColor',[1 .7 0],'MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,1),'o');
    set(p,'MarkerFaceColor','k','MarkerEdgeColor','k','MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,2),'o');
    set(p,'MarkerFaceColor','y','MarkerEdgeColor','y','MarkerSize',3);
    ylim([1320 3000])
       freezeColors
       subplot(7,1,4:6)
    %pcolor(master(midx),hgts(1,1:271),RHth(1:271,midx));shading flat;
    pcolor(time2(midc),msl,bs(midc,:)');shading flat;caxis([-8 -5]);
    colormap(jet) 
    %caxis([0 100]);
    
    hold on;
    datetick('x','dd')
    axis tight
    hold on;
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:1:400,'-','linewidth',1,'linecolor',[.7 .7 .7])
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:5:400,'-','linewidth',2,'linecolor',[.5 .5 .5])
    p=plot(master(midx),tmix(midx),'o');
    set(p,'MarkerFaceColor',[1 .7 0],'MarkerEdgeColor',[1 .7 0],'MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,1),'o');
    set(p,'MarkerFaceColor','k','MarkerEdgeColor','k','MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,2),'o');
    set(p,'MarkerFaceColor','y','MarkerEdgeColor','y','MarkerSize',3);
    ylim([1320 3000])
    %freezeColors
    %cbfreeze(colorbar);
         subplot(7,1,7)
     ylim([0 100])    
     bar(ptim(midp),haw(midp),'k');  
     hold on;
     datetick('x','dd')
     axis tight
          ylim([0 100])
        print(gcf,'-dpng','-r600',strcat('ISS_IOP_',IOPstr,'_STS','.png'))
       f=figure(IOP+300);
    clf
    set(f,'position',[scrsz(1) scrsz(2) scrsz(3) scrsz(4)*.85],'color','white')
    subplot(2,1,1)
    %pcolor(master(midx),hgts(1,1:271),cpin(:,midx));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([-150 50]);
     pcolor(master(midx),hgts(1,1:271),st(:,midx));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 4]);
          colormap(pink) 
    hold on;
    datetick('x','dd')
    axis tight
    hold on;
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:1:400,'-','linewidth',1,'linecolor',[.7 .7 .7])
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:5:400,'-','linewidth',2,'linecolor',[.5 .5 .5])
     p=plot(master(midx),tmix(midx),'o');
    set(p,'MarkerFaceColor',[1 .7 0],'MarkerEdgeColor',[1 .7 0],'MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,1),'o');
    set(p,'MarkerFaceColor','k','MarkerEdgeColor','k','MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,2),'o');
    set(p,'MarkerFaceColor','y','MarkerEdgeColor','y','MarkerSize',3);
    
    ylim([1320 3000])
    freezeColors
    cbfreeze;
    subplot(2,1,2)
    pcolor(master(midx),hgts(1,1:271),RHth(1:271,midx));shading flat;
    colormap(flipud(jet)) 
    caxis([0 100]);
    %pcolor(times(tidx1),vert,(sqrt(u(tidx1,:).^2+v(tidx1,:).^2)')); shading flat; 
    %caxis([0 15]);
    hold on;
    datetick('x','dd')
    axis tight
    hold on;
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:1:400,'-','linewidth',1,'linecolor',[.7 .7 .7])
    contour(master(midx),hgts(1,1:271),blended(1:271,midx),200:5:400,'-','linewidth',2,'linecolor',[.5 .5 .5])
     p=plot(master(midx),tmix(midx),'o');
    set(p,'MarkerFaceColor',[1 .7 0],'MarkerEdgeColor',[1 .7 0],'MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,1),'o');
    set(p,'MarkerFaceColor','k','MarkerEdgeColor','k','MarkerSize',3);
    p=plot(master(midx),mix_depth(midx,2),'o');
    set(p,'MarkerFaceColor','y','MarkerEdgeColor','y','MarkerSize',3);

    ylim([1320 3000])
    freezeColors
    cbfreeze(colorbar);
    print(gcf,'-dpng','-r600',strcat('ISS_IOP_',IOPstr,'_RH','.png'))
    
          
%------------------------------------------------------------------------
%-----------------NOW START THE PROFILE PLOTS------------------------------------
%------------------------------------------------------------------------
mm = length(midx);
 %for m = midx(1):midx(mm)
 
     for m=midx(1):midx(1)
           
           ht_tick = 1400:200:3000;
           p_tick = 650:25:850;
           rh_tick = 0:25:100;
           deftic=-15:5:15;
           minht=1320;
           maxht=3000;
           minp = PRES(1,271);
           maxp = PRES(1,1);
           th_tick=270:10:310;
           wnd_tick=0:5:20;
           bck_tick=-10:1:-5;
           cpi_tick= [ -150 -50 0 ];
                     
           fh = figure(1000);    
          
           set(fh,'PaperPositionMode','manual',...
           'PaperPosition',[0 0 2.5 4.5],'PaperSize',[2.5 4.5],...
           'PaperUnits','inches');
           set (fh, 'color', 'white');
           
%----------------------LEFT----------------------------------------------


           % theta vs. height
           
           minth = nanmin(blended(1:271,m));
           maxth= nanmax(blended(1:271,m));
           ax1 = axes('Position',[.15 .1 .7 .75],...
            'XAxisLocation','bottom',...
           'Xlim',[270 310],...
           'XTick',th_tick,...
           'XLimMode','manual',...
            'Ylim',[minht maxht],...
           'YLimMode','manual',...
           'YTick',ht_tick,...
           'YAxisLocation','left',...
            'Color','none',...
           'XColor','k','YColor','k');
          
            xlabel('Potential Temperature (K)');
            ylabel('Height (m)');
       
            hold on;
          
           ax2 = axes('Position',get(ax1,'Position'),...
           'XAxisLocation','top',...
           'Xlim',[0 100],...
           'XLimMode','manual',...
           'XTick',rh_tick,...
           'YTick',ht_tick,...
           'Ylim',[minht maxht],...
           'YLimMode','manual',...
           'YAxisLocation','left',...
            'Color','none',...
           'XColor','k','YColor','k');
                 
            %invisible plot to set properties so shadedplot does not mess
            %it all up
            line(RHth(1:271,m),hgts(1,1:271),'LineStyle','none','Parent',ax2,'LineWidth',2);
            %line(blended(1:271,m),hgts(1,1:271),'LineStyle','none','Parent',ax1,'LineWidth',2);
            line(RHth(1:271,m),hgts(1,1:271),'LineStyle','none','Parent',ax1,'LineWidth',2);
            %line(blended(1:271,m),hgts(1,1:271),'LineStyle','none','Parent',ax2,'LineWidth',2);
            
            hold on;
            clear dfm indd lin;
           dfm = diff(st(bot:269,m),1,1);
           indd = find(abs(dfm)>0);
           lin = length(indd);
           if(lin>0)
           x=0:0.452488688:100;   
   for i=1:lin
       clear y1 y2;
       if(st(indd(i),m)==stablethres) % shades stable layers
           if(i==1)
               y1 = bot * ones(size(x));
           else
           y1 = hgts(1,indd(i-1)) * ones(size(x));
           end
           y2 = hgts(1,indd(i)) * ones(size(x));
           %for h = indd(i-1):indd(i)
           %    colorp = (50+rmin(h,m))/50;
           %    colorp = median([ .9 colorp 0]);
            ha = jbfill(x,y2,y1,[1 0 0],[1 0 0],0,0.25);
            set(ha,'Parent',ax2);
          % end
       end
   end
   
      end
           
            
          
            
            %rh vs. height
            line(RHth(1:271,m),hgts(1,1:271),'Color','g','Parent',ax2,'LineWidth',2);
            xlabel('Relative Humidity (%)');
           
            hold on;
           % so that potential temp on top
            line(blended(1:271,m),hgts(1,1:271),'color','k','LineWidth',2,'Parent',ax1);
          
            
           

            
           
   % pressure axes
           
           ax6 = axes('Position',get(ax1,'Position'),...
           'XAxisLocation','bottom',...
           'XLim',[270 310],...
           'XLimMode','manual',...
           'XTick',th_tick,...
           'YAxisLocation','right',...
           'Ylim',[minp maxp],...
            'YTick',p_tick,...
           'YLimMode','manual',...
           'YDir','reverse',...
           'Color','none',...
           'XColor','k','YColor','k');
            ylabel('Pressure (mb)');
            
           
       
      
%------------------------------------------------------------------------
%-----------------END THE PLOTS------------------------------------
%------------------------------------------------------------------------


         [year,month,day,hour,min,sec] = datevec(master(m));
         titled = strcat(num2str(year),'/',num2str(month,'%02d'),'/',num2str(day,'%02d'),'  ',num2str(hour,' %02d'),'UTC');
        text('Position',[.5 3.9],'String',titled,'Units','inches','Parent',ax1);

      print(gcf,'-dpng','-r600',strcat('ISS_IOP_',IOPstr,'_P_',num2str(m,'%03d'),'.png'))
      close(1000)
     end
           
end
