%% Figure 1
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;
%%
%Loading assumptions tests
clear TVsave p
%Panel a-b (1st part)
Tvec=295; %cloud base temperatures
p0=960e2; %cloud base pressure
pf=125e2; %final pressure

opts=odeset('Refine',10,'MaxStep',10);
%R08 Tests
sm=[1, 1, 1]; %no supersaturation
ifr=[1, 1, 1]; %freeze or no freeze
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
TFinalFreeze=-4; %Finish freezing at this temperature
num_exp=length(sm);

for it=1%:length(Tvec)
    
    T0=Tvec(it);
    
    for i=1:num_exp
        satmax=sm(i);
        icefracFinal=ifr(i);
        lo=loadopt(i);
        fl=loadamt(i);
        unlT=unloadT(i);
        
        [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,icefracFinal,TFinalFreeze,lo,fl,unlT),[p0 pf],T0,opts);
        conditions=getvapor(P,T,p0,T0,satmax,icefracFinal,TFinalFreeze,lo,fl,unlT);
        [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};
        
        RH=pF1./es(T)*100;
        p(:,i)=P;
        TT(:,i)=T;
        TVsave(:,i)=T.*(1+rF1/ep)./(1+rF1+rF);
        if T0==295; TVsave2=TVsave; end
    end
end

f=figure('position',[1350 700 1000 650]);
a1=subplot(1,2,1);
hold on
h=plot([0 0],[pf p0]/100,'k-','Linewidth',1);
legendlabel(h,'off')
plot(TVsave2(:,2)-TVsave2(:,1),p(:,2)/100,'-','color',[.5 .3 .7])
plot(TVsave2(:,3)-TVsave2(:,1),p(:,2)/100,'-','color',[0 .4 .7])

clear CAPE TVsave p
icefracFinal=1;
TFinalFreeze=-40;
Tunl=273.15-4;
T0=295;

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

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

%Panel a Example profiles
lss={'-','--'};
clear p TVsave

%Ultra clean case
loadopt=5.0; %maximum condensate load is 0 g/kg
pctl=0;

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

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

%Polluted case
loadopt=1.23315; %gradual unloading between -4 and -40
pctl=1;

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

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

%Take difference
plot(TVsave(:,2)-TVsave(:,1),p(:,1)/100,'Linestyle','-','Color',[.4 .8 .6])

%Ultra clean case
loadopt=5.0; %maximum condensate load is 0 g/kg
pctl=0;

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

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

%Polluted case
loadopt=5.072; %maximum condensate load is 10 g/kg
pctl=1;

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

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

%Take difference
plot(TVsave(:,2)-TVsave(:,1),p(:,1)/100,'Linestyle','-','Color',[.4 .8 .07])
set(gca,'ydir','reverse','yscale','log','ylim',[pf p0]/100)
xlabel('T_{\rho,polluted} - T_{\rho,clean} (K)')
ylabel('Parcel Pressure (hPa)')
box off

annotation('textbox',[.058 .2373 .636 .04],'String','\bf{Polluted Parcel:}','Margin',2,'BackgroundColor','w','FontSize',12.6,'LineStyle','none')
annotation('rectangle',[.058 .0446 .6360 .2323],'LineWidth',0.5,'Color','k');
annotation('textbox',[.058 .0077 .6360 .0485],'String','\bf{Clean Parcel: Pseudo-adiabatic with freezing at -4^oC}',...
    'Margin',2,'BackgroundColor','w','FontSize',12.6,'LineWidth',0.5,'Color','k')

hl1=legend('Like R08 with instantaneous freezing and unloading at -4^oC',...
    'Like R08 but without any unloading',...
    'Gradual freezing and unloading between -4^oC and -40^oC',...
    'Gradual freezing btw -4^oC and -40^oC and unloading starting when r_t = 7.2 g/kg','Location','southoutside');
hl1.Position=[.058 .0586 .6360 .1899];
hl1.Box='off';
hl1.Color='w';

a2=axes('position',get(gca,'position'));

hold on
set(a2,'color','none','YAxisLocation','right','ylim',[pf p0]/100,'ydir','reverse','xtick',[],'yscale','log')
Tticks=-70:10:20;
Tc=T-273.15;
for i=1:length(Tticks)
    j=find(Tc<=Tticks(i),1,'first');
    ptick(i)=p(j,1)/100;
end
set(a2,'ytick',ptick,'yticklabel',Tticks)
ylabel(['Clean Parcel Temperature (^oC)'])
box off

a3=axes('position',get(gca,'position'));
hlink=linkprop([a1,a2,a3],'position');
set(a3,'color','none','xtick',[],'ytick',[])
p=get(a1,'position');
set(a1,'position',[p(1)-.05 1-.61 .32 .6])

%%
clear CAPE TVsave p
icefracFinal=0;%1;
TFinalFreeze=-40;
Tunl=273.15-4;
T0=295;
p0=960e2; %cloud base pressure
pf=125e2; %final pressure
loadopt=5.0; %maximum condensate load is 0 g/kg
pctl=0;

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

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

%Panel a Example profiles
lss={'-','--'};
clear p TVsave
b1=subplot(1,2,2);
hold on
colors=[.6 .08 .2
    1 .4 .2
    .9 .7 .1
    .9 .8 .5];
h=plot([0 0],[pf p0]/100,'k-','Linewidth',1);
legendlabel(h,'off')
sms=[1.01,1.05,1.1,1.15];
for i=1:4
    satmax=sms(i);
    %Ultra clean case
    [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,satmax,icefracFinal,TFinalFreeze,loadopt,pctl,Tunl),[p0 pf],T0,opts);
    conditions=getvapor(P,T,p0,T0,satmax,icefracFinal,TFinalFreeze,loadopt,pctl,Tunl);
    [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};
    Ta=T;
    rF1a=rF1;
    rFa=rF;
    Ta=T;
    
    p(:,1)=P;
    TVsave(:,1)=T.*(1+rF1/ep)./(1+rF1+rF);
    zsave(:,1)=cumtrapz(P,-287*TVsave(:,1)./P/9.8);
    
    %Polluted case
    sm=1;
    
    [P,T]=ode45(@(p,T) moistlapse2(p,T,p0,T0,sm,icefracFinal,TFinalFreeze,loadopt,pctl,Tunl),[p0 pf],T0,opts);
    conditions=getvapor(P,T,p0,T0,1,icefracFinal,TFinalFreeze,loadopt,pctl,Tunl);
    [liqfrac,icefrac,dliqfrac,rF1,rF,pF1]=conditions{:};
    
    p(:,2)=P;
    TVsave(:,2)=T.*(1+rF1/ep)./(1+rF1+rF);
    zsave(:,2)=cumtrapz(P,-287*TVsave(:,2)./P/9.8);
    %Take difference 
    plot(TVsave(:,2)-TVsave(:,1),p(:,1)/100,'Linestyle','-','Color',colors(i,:))
end

set(gca,'ydir','reverse','yscale','log','ylim',[pf p0]/100,'xlim',[-0.5 1.5])
set(gca,'ylim',[pf p0]/100)

xlabel('T_{\rho,polluted} - T_{\rho,clean} (K)')
ylabel('Parcel Pressure (hPa)')
box off

annotation('textbox',[.703 .1892 .203 .04],'String','\bf{Clean Parcel: }',...
    'Margin',2,'BackgroundColor','w','FontSize',12.6,'LineStyle','none')
annotation('textbox',[.703 .2296 .203 .048],'String',...
    '\bf{Polluted Parcel: S_{eq}=0%}','Margin',2,'BackgroundColor','w','FontSize',12.6,'LineWidth',0.5,'Color','k')
annotation('rectangle',[.703 .0077 .203 .2708],'LineWidth',0.5,'Color','k');
hl2=legend('S_{eq}=1%','S_{eq}=5%','S_{eq}=10%','S_{eq}=15%','Location','southoutside');
hl2.Position=[.703 .0077 .134 .1899];
hl2.Box='off';

b2=axes('position',get(gca,'position'));

hold on
set(b2,'color','none','YAxisLocation','right','ylim',[pf p0]/100,'ydir','reverse','xtick',[],'yscale','log')
Tticks=-70:10:20;
Tc=T-273.15;
for i=1:length(Tticks)
    j=find(Tc<=Tticks(i),1,'first');
    ptick(i)=p(j,1)/100;
end
set(b2,'ytick',ptick,'yticklabel',Tticks)
ylabel(['Polluted Parcel Temperature (^oC)'])
box off

b3=axes('position',get(gca,'position'));
hlink=linkprop([b1,b2,b3],'position');
set(b3,'color','none','xtick',[],'ytick',[])
p=get(b1,'position');
set(b1,'position',[p(1) 1-.61 .32 .6])