%Define saturation vapor pressures and some constants
es=@(T) 100*10.^(-7.90298*(373.16./T-1)+5.02808*log10(373.16./T)...
    -1.3816e-7*(10.^(11.344*(1-T/373.16))-1)+...
    8.1328e-3*(10.^(-3.49149*(373.16./T-1))-1)+log10(1013.246));
esi=@(T) 100*10.^(-9.09718*(273.16./T-1)-3.56654*log10(273.16./T)+...
    0.876793*(1-T/273.16)+log10(6.1071));
R1=461.5;
R0=287;
ep=R0/R1;
%% Supersaturation tests
%Figure 3
clear CAPE
Tvec=280:2:300; %cloud base temperatures to test
pf=200e2; %final pressure
TStart=-4; %start freezing at this temperature
TFinalFreeze=-40; %finish freezing at this temperature
opts=odeset('Refine',10,'MaxStep',10);

%Notation is consistent with Fruh and Wirth 2007 QJRMS. Here's a guide:
%pF1=water vapor pressure
%rF1=water vapor mixing ratio
%rF=total (ice + liquid) condensate mixing ratio

for it=1:length(Tvec)
    
    p0=960e2; %cloud base pressure
    T0=Tvec(it); %cloud base temperature
    
    %My cond only tests
    sm=[1:.05:1.2,1.2]; %equilibrium saturation ratio
    loadopt=zeros(size(sm)); %unloading option - unload starting at unlT 
    unlT=0; %Never unload, temperature will never be less than 0 
    
    num_exp=length(sm);
    
    for i=1:num_exp
        satmax=sm(i);        
        lo=loadopt(i);
        
        [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,lo,1,unlT),[p0 pf],T0,opts);
        conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,lo,1,unlT);
        [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};
        p(:,i)=P;
        %Xu and Emmanuel 1989 MWR
        TVsave(:,i)=T.*(1+rF1/ep)./(1+rF1+rF);

    end
    
    for i=1:num_exp-1
        CAPE(i,it)=trapz(p(:,i),(TVsave(:,i)-TVsave(:,num_exp))./p(:,i))*287*-1;
    end
end

T0=295;
TStart=-4;
TFinalFreeze=-40;
unlT=0;
for i=1:num_exp
    satmax=sm(i);
        
    lo=loadopt(i);
    
    [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,lo,1,unlT),[p0 pf],T0,opts);
    conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,lo,1,unlT);
    [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};
    RH=pF1./es(T)*100;
    TVsave(:,i)=T.*(1+rF1/ep)./(1+rF1+rF);
    rFsave(:,i)=rF;
    RHsave(:,i)=RH;
end

colors=get(0,'defaultaxescolororder');
colors=colors([2,3,5,1,4,6,7],:);

figure('position',[1000 650 500 500])
subplot(2,2,1)
hold on
set(gca,'colororder',colors)
plot(RHsave(:,1:num_exp-1)-100,p(:,1:num_exp-1)/100,'-')
set(gca,'ydir','reverse','yscale','log','ylim',[pf p0]/100)
xlabel('Supersaturation (%)')
ylabel('Parcel Pressure (hPa)')
grid on
legend('S0%','S5%','S10%','S15%','S20%')
plotlabel('(a)')

subplot(2,2,2)
hold on
set(gca,'colororder',colors)
plot(Tvec,CAPE,'-')
xlabel('Initial Temperature (K)')
ylabel('\Delta CAPE (J/kg)')
grid on
plotlabel('(b)')

subplot(2,2,3)
hold on
set(gca,'colororder',colors)
plot(TVsave(:,1:num_exp-1)-TVsave(:,num_exp),p(:,1:num_exp-1)/100,'-')
set(gca,'ydir','reverse','yscale','log','ylim',[pf p0]/100)
xlabel('\Delta Density Temperature (^oC)')
ylabel('Parcel Pressure (hPa)')
grid on
plotlabel('(c)')

subplot(2,2,4)
hold on
set(gca,'colororder',colors)
plot((rFsave(:,1:num_exp-1)-rFsave(:,num_exp))*1000,p(:,1:num_exp-1)/100,'-')
set(gca,'ydir','reverse','yscale','log','ylim',[pf p0]/100)
xlabel('\Delta Condensate (g kg^{-1})')
ylabel('Parcel Pressure (hPa)')
grid on
plotlabel('(d)')

%% Figure 1
%Loading assumptions tests

figure

clear TVsave p
%Panel a-b (1st part)
Tvec=280:300; %cloud base temperatures
p0=960e2; %cloud base pressure
pf=200e2; %final pressure

opts=odeset('Refine',10,'MaxStep',10);
%R08 Tests
sm=[1, 1, 1]; %no supersaturation
loadopt=[0, 0, 0]; %instantaneous at unloadT 
loadamt=[0, 0, 1]; %no loading, no loading, loading
unloadT=[400, 273.15-4, 0]; %unload immediately, unload at -4, never unload
TStart=-3;
TFinalFreeze=-4; %Finish freezing at this temperature - quasi-isothermal
num_exp=length(sm);

for it=1:length(Tvec) 
    
    T0=Tvec(it);
    
    for i=1:num_exp
        satmax=sm(i);
        lo=loadopt(i);
        fl=loadamt(i);
        unlT=unloadT(i);
        
        [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,lo,fl,unlT),[p0 pf],T0,opts);
        conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,lo,fl,unlT);
        [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};

        RH=pF1./es(T)*100;
        p(:,i)=P; 
        TVsave(:,i)=T.*(1+rF1/ep)./(1+rF1+rF);
        Tsave(:,i)=T;
        rFsave(:,i)=rF;
        rvsave(:,i)=icefrac;
        if T0==295; TVsave2=TVsave; end
    end
    
    for i=1:num_exp
        CAPE_simple(i,it)=trapz(p(:,i),(TVsave(:,i)-TVsave(:,1))./p(:,i))*287*-1;
    end
end

panel1=subplot(2,2,1);
hold on
set(gca,'colororder',[1 0 0;.6 .2 .6;.25 .25 .25;.25 .25 .25])
plot(TVsave2(:,2:end)-TVsave2(:,1),p(:,2:end)/100,'-')
set(gca,'ydir','reverse','yscale','log','ylim',[pf p0]/100)
xlabel('\Delta Density Temperature (^o C)')
ylabel('Parcel Pressure (hPa)')
plotlabel('(a)')
grid on

panel2=subplot(2,2,2);
hold on
set(gca,'colororder',[1 0 0;.6 .2 .6; .25 .25 .25;.25 .25 .25;.5 .5 .5;.5 .5 .5;.75 .75 .75;.75 .75 .75])
plot(Tvec,CAPE_simple(2:3,:),'-')
plotlabel('(b)')
grid on

%New panels c-d 
clear CAPE TVsave p
p0=960e2;
Tvec=[280,300];
pf=200e2;
satmax=1;
TStart=-4;
TFinalFreeze=-40;

opts=odeset('Refine',10,'MaxStep',10);

pct_load=0:.1:1; %final loaded fraction (in paper convert to unloaded fraction)
r_unload=0:10; %condensate mixing ratio to unload at
Tunl=0; %arbitrary value, not being used
lT=length(r_unload); 
lp=length(pct_load); 

ls={'-','--','-.',':'};

for itt=1:length(Tvec)
    T0=Tvec(itt);
    for it=1:lT
        loadopt=5+r_unload(it)/100; %see getvapor.m for loadopt interpretation
        
        parfor i=1:lp
            
            pctl=pct_load(i);
            
            [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl),[p0 pf],T0,opts);
            conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl);
            [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};
            p(:,it,i)=P;
            TVsave(:,it,i)=T.*(1+rF1/ep)./(1+rF1+rF);
        end
    end
    
    for it=1:lT
        for i=1:lp
            CAPE(it,i)=trapz(p(:,it,i),(TVsave(:,it,i)-TVsave(:,1,1))./p(:,it,i))*287;
        end
    end
    
    if T0==280 || T0==300
        if T0==280; ai=1; else; ai=2; end
        subplot(2,2,2+ai)
        contourf(1-pct_load,r_unload,CAPE*-1,-2000:50:2000)
        hold on
        contour(1-pct_load,r_unload,CAPE,[1 1],'k-','linewidth',3)
        contour(1-pct_load,r_unload,CAPE,[-1 -1],'k-','linewidth',3)
        grid on
        set(gca,'xdir','reverse');%,'ydir','reverse')
        xlabel(['Unloaded Fraction',10,'More Polluted \rightarrow'])
        ylabel(['More Polluted \rightarrow',10,'Cond. Loading Threshold (g kg^{-1})'])
        set(gca,'clim',[-600 200])
        if ai==1;plotlabel('(c) T_{CB}=280K');else;plotlabel('(d) T_{CB}=300K');end
    end
     CAPE_save(itt)=max(CAPE(:));
     CAPE_save2(itt)=min(CAPE(:));
end

%Panel b (2nd part)
clear CAPE TVsave p
p0=960e2;
Tvec=280:300;
pf=200e2;
satmax=1;
TStart=-4;
TFinalFreeze=-40;

opts=odeset('Refine',10,'MaxStep',10);

pct_load=[0,1]; %final loaded fraction (in paper convert to unloaded fraction)
r_unload=[0,10]; %condensate mixing ratio to unload at
Tunl=0; %arbitrary value, not being used
lT=length(r_unload); 
lp=length(pct_load); 

ls={'-','--','-.',':'};

for itt=1:length(Tvec)
    T0=Tvec(itt);
    for it=[1,lT]
        loadopt=5+r_unload(it)/100; %see getvapor.m for loadopt interpretation
        
        parfor i=[1,lp]
            
            pctl=pct_load(i);
            
            [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl),[p0 pf],T0,opts);
            conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl);
            [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};
            p(:,it,i)=P;
            TVsave(:,it,i)=T.*(1+rF1/ep)./(1+rF1+rF);
        end
    end
    
    for it=1:lT
        for i=1:lp
            CAPE(it,i)=trapz(p(:,it,i),(TVsave(:,it,i)-TVsave(:,1,1))./p(:,it,i))*287;
        end
    end
    
   
     CAPE_save(itt)=max(CAPE(:));
     CAPE_save2(itt)=min(CAPE(:));
end
axes(panel2)
plot(Tvec,CAPE_save2*-1,'-') %oops, took difference in wrong direction for CAPE
plot(Tvec,CAPE_save*-1,'--')
xlabel('Initial Temperature (K)')
ylabel('\Delta CAPE (J kg^{-1})')

%Panel a Example profiles

clear CAPE TVsave p
p0=960e2;
Tvec=[280,300];
pf=200e2;
satmax=1;
TStart=-4;
TFinalFreeze=-40;
opts=odeset('Refine',10,'MaxStep',10);

pct_load=0:.1:1; %final loaded fraction (in paper convert to unloaded fraction)
r_unload=0:10; %condensate mixing ratio to unload at
Tunl=0; %arbitrary value, not being used
lT=length(r_unload); 
lp=length(pct_load); 
lss={'-','--'};
for itt=1:length(Tvec)
    T0=Tvec(itt);
    clear p TVsave

    axes(panel1)

    loadopt=5.04;
    pctl=pct_load(5);

    [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl),[p0 pf],T0,opts);
    conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl);
    [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};

    p(:,2)=P;
    TVsave(:,2)=T.*(1+rF1/ep)./(1+rF1+rF);

    loadopt=5.0;
    pctl=pct_load(1);

    [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl),[p0 pf],T0,opts);
    conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl);
    [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};

    p(:,1)=P;
    TVsave(:,1)=T.*(1+rF1/ep)./(1+rF1+rF);

    plot(TVsave(:,2)-TVsave(:,1),p(:,1)/100,'Linestyle',ls{itt})
end

grid on
set(gca,'ydir','reverse','yscale','log')
ylabel('Parcel Pressure (hPa)')
xlabel('\Delta Density Temperature (^oC)')
legend('Unload, Inst. Frz., 295K','Load, Inst. Frz., 295K','Gradual, 280K','Gradual, 300K')

%Panel b TFinalFreeze -20
clear CAPE TVsave p
p0=960e2;
Tvec=280:300;
pf=200e2;
satmax=1;
TStart=-4;
TFinalFreeze=-20;

opts=odeset('Refine',10,'MaxStep',10);

pct_load=[0,1]; %final loaded fraction (in paper convert to unloaded fraction)
r_unload=[0,10]; %condensate mixing ratio to unload at
Tunl=0; %arbitrary value, not being used
lT=length(r_unload); 
lp=length(pct_load); 

ls={'-','--','-.',':'};

for itt=1:length(Tvec)
    T0=Tvec(itt);
    for it=1:lT
        loadopt=5+r_unload(it)/100; %see getvapor.m for loadopt interpretation
        
        parfor i=1:lp
            
            pctl=pct_load(i);
            
            [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl),[p0 pf],T0,opts);
            conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl);
            [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};
            p(:,it,i)=P;
            TVsave(:,it,i)=T.*(1+rF1/ep)./(1+rF1+rF);
        end
    end
    
    for it=1:lT
        for i=1:lp
            CAPE(it,i)=trapz(p(:,it,i),(TVsave(:,it,i)-TVsave(:,1,1))./p(:,it,i))*287;
        end
    end
    
    
     CAPE_save(itt)=max(CAPE(:));
     CAPE_save2(itt)=min(CAPE(:));
end

%Panel b (2nd part)
axes(panel2)
plot(Tvec,CAPE_save2*-1,'-') %oops, took difference in wrong direction for CAPE
plot(Tvec,CAPE_save*-1,'--')
xlabel('Initial Temperature (K)')
ylabel('\Delta CAPE (J kg^{-1})')

p0=960e2;
Tvec=280:300;
pf=200e2;
satmax=1;
TStart=-15;
TFinalFreeze=-40;

opts=odeset('Refine',10,'MaxStep',10);

pct_load=[0,1]; %final loaded fraction (in paper convert to unloaded fraction)
r_unload=[0,10]; %condensate mixing ratio to unload at
Tunl=0; %arbitrary value, not being used
lT=length(r_unload); 
lp=length(pct_load); 

ls={'-','--','-.',':'};

for itt=1:length(Tvec)
    T0=Tvec(itt);
    for it=1:lT
        loadopt=5+r_unload(it)/100; %see getvapor.m for loadopt interpretation
        
        parfor i=1:lp
            
            pctl=pct_load(i);
            
            [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl),[p0 pf],T0,opts);
            conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl);
            [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};
            p(:,it,i)=P;
            TVsave(:,it,i)=T.*(1+rF1/ep)./(1+rF1+rF);
        end
    end
    
    for it=1:lT
        for i=1:lp
            CAPE(it,i)=trapz(p(:,it,i),(TVsave(:,it,i)-TVsave(:,1,1))./p(:,it,i))*287;
        end
    end
    
    
     CAPE_save(itt)=max(CAPE(:));
     CAPE_save2(itt)=min(CAPE(:));
end

%Panel b (2nd part)
axes(panel2)
plot(Tvec,CAPE_save2*-1,'-') %oops, took difference in wrong direction for CAPE
plot(Tvec,CAPE_save*-1,'--')
xlabel('Initial Temperature (K)')
ylabel('\Delta CAPE (J kg^{-1})')
%% Figure 4
%combined supersaturation and loading tests

%Panels a-b
clear CAPE 

p0=960e2;
pf=200e2;
TStart=-4;
TFinalFreeze=-40;

opts=odeset('Refine',10,'MaxStep',10);

pct_load=0:.1:1; %loaded fraction
T0s=[280, 300]; %cloud base temperatures
Tunl=0; %arbitrary value, Tunl not used
r_unload=0:10; %condensate loading thresholds, g/kg
lT=length(r_unload);
satmaxes=[1.2:-.02:1]; %equilibirum supersaturations
lsat=length(satmaxes);

figure('position',[1000 650 700 700])
hold on
ls={'-','--','-.',':'};

for TT=1:2
    T0=T0s(TT);    
    for i=1:lT
        pctl=pct_load(i);
        loadopt=5+r_unload(i)/100;        
        for is=1:lsat
            satmax=satmaxes(is);            
            [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl),[p0 pf],T0,opts);
            conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl);
            [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};
            RH=pF1./es(T)*100;
            TVsave(:,i,is)=T.*(1+rF1/ep)./(1+rF1+rF);
        end
    end
    
    refTV=TVsave(:,1,1);
    for i=1:lT
        for is=1:lsat
            CAPE(i,is,TT)=trapz(P,(TVsave(:,i,is)-refTV)./P)*287;
        end
    end
    
    subplot(2,2,TT)    
    contourf((satmaxes-1)*100,r_unload,CAPE(:,:,TT)*-1,-2000:50:2000)
    hold on
    %Add bold 0 line
    contour((satmaxes-1)*100,r_unload,CAPE(:,:,TT)*-1,[0 0],'k-','linewidth',2)
    grid on
    set(gca,'xdir','reverse','clim',[-600 400])
    xlabel(['Supersaturation (%)',10,'More Polluted \rightarrow'])
    ylabel(['More Polluted \rightarrow',10,'Loading'])
    if TT==1; plotlabel(['(a) T_{CB} = ',num2str(T0),' K']);end
    if TT==2; plotlabel(['(b) T_{CB} = ',num2str(T0),' K']);end
end


%Panel c

ls={'--','-','--','-'};
lci=[1,1,2,2];
T0=300;
combos=[3 4;
        4 5;
        4 10;
        9 5;
        9 10;];    

clear TVsave

subplot(2,2,3)
hold on
for i=1:5
    pctl=0.5;
    it=combos(i,1);
    is=combos(i,2);
    satmax=satmaxes(is);
    loadopt=5+r_unload(it)/100;
    [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl),[p0 pf],T0,opts);
    conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl);
    [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};

    TVsave(:,i)=T.*(1+rF1/ep)./(1+rF1+rF);
    if i>1
        set(gca,'ColorOrderIndex',lci(i-1))
        plot(TVsave(:,i)-TVsave(:,1),P/100,'Linestyle',ls{i-1})
    end
end

set(gca,'ydir','reverse','yscale','log')
ylabel('Pressure (hPa)')
xlabel('\Delta Density Temperature (^oC)')
grid on
plotlabel('(c)')
legend('Small \Deltar_t; Small \Deltass','Small \Deltar_t; Large \Deltass',...
    'Large \Deltar_t; Small \Deltass','Large \Deltar_t; Large \Deltass')


%Panel d

clear TVsave

subplot(2,2,4)
hold on
for i=1:5
    
    it=combos(i,1);
    is=combos(i,2);
    satmax=satmaxes(is);
    loadopt=5+r_unload(8)/100;
    pctl=pct_load(it);
    [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl),[p0 pf],T0,opts);
    conditions=getvapor(P,T,p0,T0,satmax,TStart,TFinalFreeze,loadopt,pctl,Tunl);
    [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};
    RH=pF1./es(T)*100;

    TVsave(:,i)=T.*(1+rF1/ep)./(1+rF1+rF);
     if i>1
        set(gca,'ColorOrderIndex',lci(i-1)+2)
        plot(TVsave(:,i)-TVsave(:,1),P/100,'Linestyle',ls{i-1})
    end
end

set(gca,'ydir','reverse','yscale','log')
ylabel('Pressure (hPa)')
xlabel('\Delta Density Temperature (^oC)')
grid on
plotlabel('(d)')
legend('Small \Deltaf; Small \Deltass','Small \Deltaf; Large \Deltass',...
    'Large \Deltaf; Small \Deltass','Large \Deltaf; Large \Deltass')