%Created 2012/11/1
%This script takes the profile legs from the 1Hz dataset and calculates the
%profiles of acc and auto as a function of height

function profile_cloudMICROplots

datafile='1Hz_profile_data_rev.nc';
microdatafile='1Hz_profile_micro2.nc';
microinterpfile='1Hz_profile_micro_interp1.nc';
microinterp20file='1Hz_profile_micro_interp1_rthresh20.nc';
microinterp3file='1Hz_profile_micro_interp3.nc';

profile=nc_varget(datafile,'Profile');
lat=nc_varget(datafile,'LAT');
lon=nc_varget(datafile,'LON');
lwc=nc_varget(datafile,'XGLWC');
lwccdp=nc_varget(datafile,'LWCCDP');
lwc2dc=nc_varget(datafile,'LWC2DC');
concdp=nc_varget(datafile,'CONCDP');
p=nc_varget(datafile,'PRES');
T=nc_varget(datafile,'TEMP');
waterVaporP=nc_varget(datafile,'EDPC');

autoSCEL10=nc_varget(microdatafile,'AutoL10SCE');
autoSCEL5=nc_varget(microdatafile,'AutoL5SCE');
accSCEL10=nc_varget(microdatafile,'AccL10SCE');
accSCEL5=nc_varget(microdatafile,'AccL5SCE');

autoSCEL10_interp=nc_varget(microinterpfile,'AutoL10SCE');
accSCEL10_interp=nc_varget(microinterpfile,'AccL10SCE');

autoSCEL10_interp20=nc_varget(microinterp20file,'AutoL10SCE');
accSCEL10_interp20=nc_varget(microinterp20file,'AccL10SCE');


autoSCEL10_interp3=nc_varget(microinterp3file,'AutoL10SCE');
accSCEL10_interp3=nc_varget(microinterp3file,'AccL10SCE');


NLoss_autoSCEL10_interp=nc_varget(microinterpfile,'AutoNumL10SCE');
NLoss_accSCEL10_interp=nc_varget(microinterpfile,'AccNumL10SCE');

eps= .622004944; % molar mass of water vapor/molar mass of dry air, Rd/Rv
Tv=T./(1- waterVaporP./p.*(1-eps));
Rd=287.05; % J/kgK
rhoAir = 100*p./(Rd*Tv); % 100 Pa/hPa* hPa /(J/kgK *K)
qD=lwc2dc./1000;     %kg m^-3
qL=lwccdp./1000;     %kg m^-3
Nd=concdp.*1000000;  %m^-3
% lwccdp2dc=lwccdp+lwc2dc;
lwccdp2dc=lwccdp;
z=nc_varget(datafile,'Z');
lwc(lwc<0)=0;
medlon=nanmean(lon,2);
medlat=nanmean(lat,2);
figure;
plot(medlon,medlat,'.')

medlon(medlat<-21 | medlat>-19)=NaN;
% medlon(medlat<-24)=NaN;
Region={'R7075','R7580','R8085'};% Name the 3 regions. Used for naming variables
RegionLim=[-70 -75;
    -75 -80;
    -80 -85]; %DEFINE LIMITS TO THE THREE REGIONS
for xx=1:3
    RegionName=Region{xx};
    in.RegionIndex.(RegionName)=find(medlon<RegionLim(xx,1) & medlon>=RegionLim(xx,2));
end

FN01='KK00';
fn01Handle=str2func(FN01);

for xx=1:3
    RegionName=Region{xx};
    in.lwc.(RegionName)=lwc(in.RegionIndex.(RegionName),:);
    in.profilenum.(RegionName)=profile(in.RegionIndex.(RegionName));
    in.z.(RegionName)=z(in.RegionIndex.(RegionName),:);
    in.lon.(RegionName)=lon(in.RegionIndex.(RegionName),:);
    in.lat.(RegionName)=lat(in.RegionIndex.(RegionName),:);
    in.lwccdp.(RegionName)=lwccdp(in.RegionIndex.(RegionName),:);
    in.lwc2dc.(RegionName)=lwc2dc(in.RegionIndex.(RegionName),:);
    in.lwccdp2dc.(RegionName)=lwccdp2dc(in.RegionIndex.(RegionName),:);
    in.qL.(RegionName)=qL(in.RegionIndex.(RegionName),:);
    in.qD.(RegionName)=qD(in.RegionIndex.(RegionName),:);
    in.Nd.(RegionName)=Nd(in.RegionIndex.(RegionName),:);
    in.rhoAir.(RegionName)=rhoAir(in.RegionIndex.(RegionName),:);
    
    in.autoL5.(RegionName)=autoSCEL5(in.RegionIndex.(RegionName),:);
    in.autoL10.(RegionName)=autoSCEL10(in.RegionIndex.(RegionName),:);
    in.accL5.(RegionName)=accSCEL5(in.RegionIndex.(RegionName),:);
    in.accL10.(RegionName)=accSCEL10(in.RegionIndex.(RegionName),:);
    
    in.autoL10interp.(RegionName)=autoSCEL10_interp(in.RegionIndex.(RegionName),:);
    in.accL10interp.(RegionName)=accSCEL10_interp(in.RegionIndex.(RegionName),:);
    in.NLoss.autoL10interp.(RegionName)=NLoss_autoSCEL10_interp(in.RegionIndex.(RegionName),:);
    in.NLoss.accL10interp.(RegionName)=NLoss_accSCEL10_interp(in.RegionIndex.(RegionName),:);
    
    in.autoL10interp20.(RegionName)=autoSCEL10_interp20(in.RegionIndex.(RegionName),:);
    in.accL10interp20.(RegionName)=accSCEL10_interp20(in.RegionIndex.(RegionName),:);
    %----------third order interpolation---------------------------
    in.autoL10interp3.(RegionName)=autoSCEL10_interp3(in.RegionIndex.(RegionName),:);
    in.accL10interp3.(RegionName)=accSCEL10_interp3(in.RegionIndex.(RegionName),:);
end

for xx=1:3
    RegionName=Region{xx};
    for ii=1:length(in.RegionIndex.(RegionName))
        [in.auto.(RegionName)(ii,:) in.acc.(RegionName)(ii,:)]=fn01Handle(in.qL.(RegionName)(ii,:),...
        in.qD.(RegionName)(ii,:),in.Nd.(RegionName)(ii,:),...
        in.rhoAir.(RegionName)(ii,:));
    end
end

FNa='profile_cloudTB';
fnaHandle=str2func(FNa);

for xx=1:3
    RegionName=Region{xx};
    [C,I]=max(in.lwc.(RegionName),[],2);
    in.maxlwc.(RegionName)=C;
    in.maxlwccdp2dc.(RegionName)=NaN(length(C),1);
    for jj=1:length(C)
        in.maxlwccdp2dc.(RegionName)(jj)=in.lwccdp2dc.(RegionName)(jj,I(jj));
    end
    in.meanlon.(RegionName)=nanmean(in.lon.(RegionName),2);
    in.meanlat.(RegionName)=nanmean(in.lat.(RegionName),2);
    for ii=1:length(in.RegionIndex.(RegionName))
        [cltop clbottom DIR]=fnaHandle(in.lwc.(RegionName)(ii,:),in.z.(RegionName)(ii,:));
        if DIR==1
            in.cloudstart.(RegionName)(ii)=clbottom;
            in.cloudend.(RegionName)(ii)=cltop;
            in.cltop.(RegionName)(ii)=in.z.(RegionName)(ii,cltop);
            in.clbottom.(RegionName)(ii)=in.z.(RegionName)(ii,clbottom);
            dZ=in.z.(RegionName)(ii,clbottom+1:cltop+1)-in.z.(RegionName)(ii,clbottom:cltop);
            in.LWP.(RegionName)(ii)=nansum(in.lwc.(RegionName)(ii,clbottom:cltop).*dZ);
            in.meanNd.(RegionName)(ii)=nanmean(in.Nd.(RegionName)(ii,clbottom:cltop));
            in.LWPCDP.(RegionName)(ii)=nansum(in.lwccdp.(RegionName)(ii,clbottom:cltop).*dZ);
            in.LWP2DC.(RegionName)(ii)=nansum(in.lwc2dc.(RegionName)(ii,clbottom:cltop).*dZ);
            in.integratedauto.(RegionName)(ii)=nansum(in.auto.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            in.integratedacc.(RegionName)(ii)=nansum(in.acc.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            %SCE integrated auto and accretion
            in.integratedautoL5.(RegionName)(ii)=nansum(in.autoL5.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            in.integratedautoL10.(RegionName)(ii)=nansum(in.autoL10.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            in.integratedaccL5.(RegionName)(ii)=nansum(in.accL5.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            in.integratedaccL10.(RegionName)(ii)=nansum(in.accL10.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);

            in.integratedautoL10interp.(RegionName)(ii)=nansum(in.autoL10interp.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            in.integratedaccL10interp.(RegionName)(ii)=nansum(in.accL10interp.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            in.integratedautoNLossL10interp.(RegionName)(ii)=nansum(in.NLoss.autoL10interp.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            in.integratedaccNLossL10interp.(RegionName)(ii)=nansum(in.NLoss.accL10interp.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            
            %---------r_thresh == 20
            in.integratedautoL10interp20.(RegionName)(ii)=nansum(in.autoL10interp20.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            in.integratedaccL10interp20.(RegionName)(ii)=nansum(in.accL10interp20.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            
            %---------third order fit
            in.integratedautoL10interp3.(RegionName)(ii)=nansum(in.autoL10interp3.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);
            in.integratedaccL10interp3.(RegionName)(ii)=nansum(in.accL10interp3.(RegionName)(ii,clbottom:cltop).*dZ)./sum(dZ);

            if in.LWPCDP.(RegionName)(ii)>0
                in.LWPcombined.(RegionName)(ii)=nansum(in.lwc2dc.(RegionName)(ii,clbottom:cltop).*dZ+in.lwccdp.(RegionName)(ii,clbottom:cltop).*dZ);
            else
                in.LWPcombined.(RegionName)(ii)=NaN;
            end
            in.adiabaticLWP.(RegionName)(ii)=1e-3.*(in.z.(RegionName)(ii,cltop)-in.z.(RegionName)(ii,clbottom)).^2;
        else
            if DIR==0
                in.cloudstart.(RegionName)(ii)=cltop;
                in.cloudend.(RegionName)(ii)=clbottom;
                in.cltop.(RegionName)(ii)=in.z.(RegionName)(ii,cltop);
                in.clbottom.(RegionName)(ii)=in.z.(RegionName)(ii,clbottom);
                dZ=in.z.(RegionName)(ii,cltop-1:clbottom-1)-in.z.(RegionName)(ii,cltop:clbottom);
                in.LWP.(RegionName)(ii)=nansum(in.lwc.(RegionName)(ii,cltop:clbottom).*dZ);
                in.LWPCDP.(RegionName)(ii)=nansum(in.lwccdp.(RegionName)(ii,cltop:clbottom).*dZ);
                in.LWP2DC.(RegionName)(ii)=nansum(in.lwc2dc.(RegionName)(ii,cltop:clbottom).*dZ);
                in.integratedauto.(RegionName)(ii)=nansum(in.auto.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                in.integratedacc.(RegionName)(ii)=nansum(in.acc.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                %SCE integrated auto and accretion
                in.integratedautoL5.(RegionName)(ii)=nansum(in.autoL5.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                in.integratedautoL10.(RegionName)(ii)=nansum(in.autoL10.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                in.integratedaccL5.(RegionName)(ii)=nansum(in.accL5.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                in.integratedaccL10.(RegionName)(ii)=nansum(in.accL10.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);

                in.integratedautoL10interp.(RegionName)(ii)=nansum(in.autoL10interp.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                in.integratedaccL10interp.(RegionName)(ii)=nansum(in.accL10interp.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                in.integratedNLossautoL10interp.(RegionName)(ii)=nansum(in.NLoss.autoL10interp.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                in.integratedNLossaccL10interp.(RegionName)(ii)=nansum(in.NLoss.accL10interp.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                %---------------r_thresh = 20
                in.integratedautoL10interp20.(RegionName)(ii)=nansum(in.autoL10interp20.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                in.integratedaccL10interp20.(RegionName)(ii)=nansum(in.accL10interp20.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                %---------------third order fit
                in.integratedautoL10interp3.(RegionName)(ii)=nansum(in.autoL10interp3.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                in.integratedaccL10interp3.(RegionName)(ii)=nansum(in.accL10interp3.(RegionName)(ii,cltop:clbottom).*dZ)./sum(dZ);
                
                if in.LWPCDP.(RegionName)(ii)>0
                    in.LWPcombined.(RegionName)(ii)=nansum(in.lwc2dc.(RegionName)(ii,cltop:clbottom).*dZ+in.lwccdp.(RegionName)(ii,cltop:clbottom).*dZ);
                else
                    in.LWPcombined.(RegionName)(ii)=NaN;
                end
                in.adiabaticLWP.(RegionName)(ii)=1e-3.*(in.z.(RegionName)(ii,cltop)-in.z.(RegionName)(ii,clbottom)).^2;
            else
                in.cloudstart.(RegionName)(ii)=NaN;
                in.cloudend.(RegionName)(ii)=NaN;
                in.cltop.(RegionName)(ii)=NaN;
                in.clbottom.(RegionName)(ii)=NaN;
                in.LWP.(RegionName)(ii)=NaN;
                in.LWPCDP.(RegionName)(ii)=NaN;
                in.LWP2DC.(RegionName)(ii)=NaN;
                in.adiabaticLWP.(RegionName)(ii)=NaN;
                in.LWPcombined.(RegionName)(ii)=NaN;
                in.integratedacc.(RegionName)(ii)=NaN;
                in.integratedauto.(RegionName)(ii)=NaN;
                in.integratedautoL5.(RegionName)(ii)=NaN;
                in.integratedautoL10.(RegionName)(ii)=NaN;
                in.integratedaccL5.(RegionName)(ii)=NaN;
                in.integratedaccL10.(RegionName)(ii)=NaN;
                in.integratedautoL10interp.(RegionName)(ii)=NaN;
                in.integratedaccL10interp.(RegionName)(ii)=NaN;
                in.integratedautoNLossL10interp.(RegionName)(ii)=NaN;
                in.integratedaccNLossL10interp.(RegionName)(ii)=NaN;
                in.integratedautoL10interp20.(RegionName)(ii)=NaN;
                in.integratedaccL10interp20.(RegionName)(ii)=NaN;
                in.integratedautoL10interp3.(RegionName)(ii)=NaN;
                in.integratedaccL10interp3.(RegionName)(ii)=NaN;
            end
        end
    end
end

scrsz = get(0,'ScreenSize');
fg_size_3panel=[1 scrsz(4)/6 3/4*scrsz(3) 3/4*scrsz(4)];

figure;
plot(in.meanlon.R7075,in.meanlat.R7075,'.b')
hold on;
plot(in.meanlon.R7580,in.meanlat.R7580,'.g')
hold on;
plot(in.meanlon.R8085,in.meanlat.R8085,'.r')
xlabel('Longitude')
ylabel('Latitude')
title('Profile locations')


for xx=1:3
    RegionName=Region{xx};
%     figure;
%     in.cloudonly.accL10interp.(RegionName)=NaN(size(in.z.(RegionName),1),10);
    for jj=1:size(in.z.(RegionName),1)
        if in.cloudstart.(RegionName)(jj)>0
            Acc=in.accL10interp.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)) ...
                ./nanmean(in.accL10interp.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)));
            Zcloud=(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)) ...
                -min(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)))) ...
                ./(abs(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj))-in.z.(RegionName)(jj,in.cloudend.(RegionName)(jj))));
            for ii=1:10
                in.cloudonly.accL10interp.(RegionName)(jj,ii)=nanmean(Acc(Zcloud>(ii-1)*0.1 & Zcloud<=ii*0.1));
            end
%             plot(in.cloudonly.accL10interp.(RegionName)(jj,:),0.05:0.1:0.95)
%             hold on;
        end
    end
%     plot(nanmedian(in.cloudonly.accL10interp.(RegionName)),0.05:0.1:0.95,'r','LineWidth',2)
%     hold on;
%     plot(prctile((in.cloudonly.accL10interp.(RegionName)),[25 75]),0.05:0.1:0.95,'r','LineWidth',2)
%     ylim([0 1]);
%     set(gca,'XScale','log');
%     xlabel('Accretion rate/mean accretion')
%     ylabel('Relative height in cloud')
%     title(['Region ' RegionName])
end


for xx=1:3
    RegionName=Region{xx};
%     figure;
%     in.cloudonly.autoL10interp.(RegionName)=NaN(size(in.z.(RegionName),1),10);
    for jj=1:size(in.z.(RegionName),1)
        if in.cloudstart.(RegionName)(jj)>0
            Auto=in.autoL10interp.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)) ...
                ./nanmean(in.autoL10interp.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)));
            Zcloud=(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)) ...
                -min(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)))) ...
                ./(abs(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj))-in.z.(RegionName)(jj,in.cloudend.(RegionName)(jj))));
            for ii=1:10
                in.cloudonly.autoL10interp.(RegionName)(jj,ii)=nanmean(Auto(Zcloud>(ii-1)*0.1 & Zcloud<=ii*0.1));
            end
%             plot(in.cloudonly.autoL10interp.(RegionName)(jj,:),0.05:0.1:0.95)
%             hold on;
        end
    end
%     plot(nanmedian(in.cloudonly.autoL10interp.(RegionName)),0.05:0.1:0.95,'r','LineWidth',2)
%     hold on;
%     plot(prctile((in.cloudonly.autoL10interp.(RegionName)),[25 75]),0.05:0.1:0.95,'r','LineWidth',2)
%     ylim([0 1]);
%     set(gca,'XScale','log');
%     xlabel('Autoconversion rate/mean autoconversion')
%     ylabel('Relative height in cloud')
%     title(['Region ' RegionName])
end



for xx=1:3
    RegionName=Region{xx};
%     figure;
%     in.cloudonly.NLoss.accL10interp.(RegionName)=NaN(size(in.z.(RegionName),1),10);
    for jj=1:size(in.z.(RegionName),1)
        if in.cloudstart.(RegionName)(jj)>0
            Acc=in.NLoss.accL10interp.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)) ...
                ./nanmean(in.NLoss.accL10interp.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)));
            Zcloud=(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)) ...
                -min(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)))) ...
                ./(abs(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj))-in.z.(RegionName)(jj,in.cloudend.(RegionName)(jj))));
            for ii=1:10
                in.cloudonly.NLoss.accL10interp.(RegionName)(jj,ii)=nanmean(Acc(Zcloud>(ii-1)*0.1 & Zcloud<=ii*0.1));
            end
%             plot(in.cloudonly.NLoss.accL10interp.(RegionName)(jj,:),0.05:0.1:0.95)
%             hold on;
        end
    end
%     plot(nanmedian(in.cloudonly.NLoss.accL10interp.(RegionName)),0.05:0.1:0.95,'r','LineWidth',2)
%     hold on;
%     plot(prctile((in.cloudonly.NLoss.accL10interp.(RegionName)),[25 75]),0.05:0.1:0.95,'r','LineWidth',2)
%     ylim([0 1]);
%     set(gca,'XScale','log');
%     xlabel('Accretion # Loss/mean accretion')
%     ylabel('Relative height in cloud')
%     title(['N Loss Region ' RegionName])
end


for xx=1:3
    RegionName=Region{xx};
%     figure;
%     in.cloudonly.NLoss.autoL10interp.(RegionName)=NaN(size(in.z.(RegionName),1),10);
    for jj=1:size(in.z.(RegionName),1)
        if in.cloudstart.(RegionName)(jj)>0
            Auto=in.NLoss.autoL10interp.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)) ...
                ./nanmean(in.NLoss.autoL10interp.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)));
            Zcloud=(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)) ...
                -min(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj):in.cloudend.(RegionName)(jj)))) ...
                ./(abs(in.z.(RegionName)(jj,in.cloudstart.(RegionName)(jj))-in.z.(RegionName)(jj,in.cloudend.(RegionName)(jj))));
            for ii=1:10
                in.cloudonly.NLoss.autoL10interp.(RegionName)(jj,ii)=nanmean(Auto(Zcloud>(ii-1)*0.1 & Zcloud<=ii*0.1));
            end
%             plot(in.cloudonly.NLoss.autoL10interp.(RegionName)(jj,:),0.05:0.1:0.95)
%             hold on;
        end
    end
%     plot(nanmedian(in.cloudonly.NLoss.autoL10interp.(RegionName)),0.05:0.1:0.95,'r','LineWidth',2)
%     hold on;
%     plot(prctile((in.cloudonly.NLoss.autoL10interp.(RegionName)),[25 75]),0.05:0.1:0.95,'r','LineWidth',2)
%     ylim([0 1]);
%     set(gca,'XScale','log');
%     xlabel('Autoconversion # Loss/mean autoconversion')
%     ylabel('Relative height in cloud')
%     title(['N Loss Region ' RegionName])
end
RegionName='R7075';
figure;
plot([1e-14 1e-6],[1e-14 1e-6],'--k')
hold on;
for xx=1:3
    RegionName=Region{xx};
    AAA=plot(in.integratedautoL10interp.(RegionName),in.integratedauto.(RegionName),'.b');
    hold on;
    BBB=plot(in.integratedautoL10interp20.(RegionName),in.integratedauto.(RegionName),'.r');
end
xlabel('SCE Auto (kg m^{-3} s^{-1})','FontSize',14,'FontWeight','demi')
ylabel('KK00 Auto','FontSize',14,'FontWeight','demi')
xlim([1e-13 1e-7]);
ylim([1e-13 1e-7]);
set(gca,'XScale','log');
set(gca,'YScale','log');
legend([AAA BBB],'r_{thresh}=25','r_{thresh}=20',...
    'Location','NorthWest');

figure;
plot([min(in.integratedaccL10interp.(RegionName)) max(in.integratedaccL10interp.(RegionName))],...
    [min(in.integratedaccL10interp.(RegionName)) max(in.integratedaccL10interp.(RegionName))],'--k')
hold on;
for xx=1:3
    RegionName=Region{xx};
    AAA=plot(in.integratedaccL10interp.(RegionName),in.integratedacc.(RegionName),'.b');
    hold on;
    BBB=plot(in.integratedaccL10interp20.(RegionName),in.integratedacc.(RegionName),'.r');
end
xlabel('SCE Acc (kg m^{-3} s^{-1})','FontSize',14,'FontWeight','demi')
ylabel('KK00 Acc','FontSize',14,'FontWeight','demi')
set(gca,'XScale','log');
set(gca,'YScale','log');
xlim([1e-13 1e-5]);
ylim([1e-13 1e-5]);
legend([AAA BBB],'r_{thresh}=25','r_{thresh}=20',...
    'Location','NorthWest');

figure;
plot([1e-13 1e-7],[1 1],':k')
hold on;
plot([1e-13 1e-7],[0.5 0.5],':k')
hold on;
for xx=1:3
    RegionName=Region{xx};
    plot(in.integratedautoL10interp.(RegionName),in.integratedautoL10interp20.(RegionName)./in.integratedautoL10interp.(RegionName),'.b')
    hold on;
end
xlabel('SCE Auto (r25) (kg m^{-3} s^{-1})','FontSize',14,'FontWeight','demi')
ylabel('Auto(r20)/Auto(r25)','FontSize',14,'FontWeight','demi')
set(gca,'XScale','log');
set(gca,'YScale','log');
xlim([1e-13 1e-7]);

figure;
plot([1e-13 1e-6],[1 1],':k')
hold on;
plot([1e-13 1e-6],[0.5 0.5],':k')
hold on;
for xx=1:3
    RegionName=Region{xx};
    plot(in.integratedaccL10interp.(RegionName),in.integratedaccL10interp20.(RegionName)./in.integratedaccL10interp.(RegionName),'.b')
    hold on;
end
xlabel('SCE Acc (r25) (kg m^{-3} s^{-1})','FontSize',14,'FontWeight','demi')
ylabel('Acc(r20)/Acc(r25)','FontSize',14,'FontWeight','demi')
set(gca,'XScale','log');
set(gca,'YScale','log');
xlim([1e-13 1e-6]);



%---------------Testing 1st order or 3rd order fit
figure;
plot([1e-14 1e-6],[1e-14 1e-6],'--k')
hold on;
for xx=1:3
    RegionName=Region{xx};
    AAA=plot(in.integratedautoL10interp.(RegionName),in.integratedautoL10interp3.(RegionName),'.b');
end
xlabel('SCE-1 Auto (kg m^{-3} s^{-1})','FontSize',14,'FontWeight','demi')
ylabel('SCE-3 Auto','FontSize',14,'FontWeight','demi')
xlim([1e-13 1e-7]);
ylim([1e-13 1e-7]);
set(gca,'XScale','log');
set(gca,'YScale','log');

figure;
plot([min(in.integratedaccL10interp.(RegionName)) max(in.integratedaccL10interp.(RegionName))],...
    [min(in.integratedaccL10interp.(RegionName)) max(in.integratedaccL10interp.(RegionName))],'--k')
hold on;
for xx=1:3
    RegionName=Region{xx};
    AAA=plot(in.integratedaccL10interp.(RegionName),in.integratedaccL10interp3.(RegionName),'.b');
end
xlabel('SCE-1 Acc (kg m^{-3} s^{-1})','FontSize',14,'FontWeight','demi')
ylabel('SCE-3 Acc','FontSize',14,'FontWeight','demi')
set(gca,'XScale','log');
set(gca,'YScale','log');
xlim([1e-13 1e-5]);
ylim([1e-13 1e-5]);


figure;
for xx=1:3
    RegionName=Region{xx};
    AAA=plot(in.LWP.(RegionName),in.integratedaccL10interp.(RegionName)./in.integratedautoL10interp.(RegionName),'.b');
    hold on;
    BBB=plot(in.LWP.(RegionName),in.integratedaccL10interp20.(RegionName)./in.integratedautoL10interp20.(RegionName),'.r');
    hold on;
end
xlabel('LWP (g m^{-2})','FontSize',14,'FontWeight','demi')
ylabel('SCE Acc/Auto','FontSize',14,'FontWeight','demi')
set(gca,'XScale','log');
set(gca,'YScale','log');
legend([AAA BBB],'r_{thresh}=25','r_{thresh}=20',...
    'Location','NorthWest');

in.cloudonly.accL10interp.all=cat(1,in.cloudonly.accL10interp.R7075,in.cloudonly.accL10interp.R7580,in.cloudonly.accL10interp.R8085);
in.cloudonly.autoL10interp.all=cat(1,in.cloudonly.autoL10interp.R7075,in.cloudonly.autoL10interp.R7580,in.cloudonly.autoL10interp.R8085);

figure;
subplot(1,2,1);
plot(in.cloudonly.accL10interp.all,0.05:0.1:0.95,'b','LineWidth',1)
hold on;
plot(nanmedian(in.cloudonly.accL10interp.all),0.05:0.1:0.95,'k','LineWidth',3)
hold on;
plot(prctile(in.cloudonly.accL10interp.all,[25 75]),0.05:0.1:0.95,'k','LineWidth',2)
ylim([0 1]);
xlim([1e-3 1e1]);
set(gca,'XScale','log');
xlabel('Acc/<Acc>','FontSize',14,'FontWeight','demi')
ylabel('Relative height in cloud','FontSize',14,'FontWeight','demi')
title('Mass conversion','FontSize',14,'FontWeight','demi')
subplot(1,2,2);
plot(in.cloudonly.autoL10interp.all,0.05:0.1:0.95,'b','LineWidth',1)
hold on;
plot(nanmedian(in.cloudonly.autoL10interp.all),0.05:0.1:0.95,'k','LineWidth',3)
hold on;
plot(prctile(in.cloudonly.autoL10interp.all,[25 75]),0.05:0.1:0.95,'k','LineWidth',2)
ylim([0 1]);
xlim([1e-3 1e1]);
set(gca,'XScale','log');
xlabel('Auto/<Auto>','FontSize',14,'FontWeight','demi')


in.cloudonly.NLoss.accL10interp.all=cat(1,in.cloudonly.NLoss.accL10interp.R7075,in.cloudonly.NLoss.accL10interp.R7580,in.cloudonly.NLoss.accL10interp.R8085);
in.cloudonly.NLoss.autoL10interp.all=cat(1,in.cloudonly.NLoss.autoL10interp.R7075,in.cloudonly.NLoss.autoL10interp.R7580,in.cloudonly.NLoss.autoL10interp.R8085);

figure;
subplot(1,2,1);
plot(in.cloudonly.NLoss.accL10interp.all,0.05:0.1:0.95,'b','LineWidth',1)
hold on;
plot(nanmedian(in.cloudonly.NLoss.accL10interp.all),0.05:0.1:0.95,'k','LineWidth',3)
hold on;
plot(prctile(in.cloudonly.NLoss.accL10interp.all,[25 75]),0.05:0.1:0.95,'k','LineWidth',2)
ylim([0 1]);
xlim([1e-2 1e1]);
set(gca,'XScale','log');
xlabel('Acc#/<Acc#>','FontSize',14,'FontWeight','demi')
ylabel('Relative height in cloud','FontSize',14,'FontWeight','demi')
title('Number conversion','FontSize',14,'FontWeight','demi')
subplot(1,2,2);
plot(in.cloudonly.NLoss.autoL10interp.all,0.05:0.1:0.95,'b','LineWidth',1)
hold on;
plot(nanmedian(in.cloudonly.NLoss.autoL10interp.all),0.05:0.1:0.95,'k','LineWidth',3)
hold on;
plot(prctile(in.cloudonly.NLoss.autoL10interp.all,[25 75]),0.05:0.1:0.95,'k','LineWidth',2)
ylim([0 1]);
xlim([1e-2 1e1]);
set(gca,'XScale','log');
xlabel('Auto#/<Auto#>','FontSize',14,'FontWeight','demi')

end