%% Figure 2
%load RRTMGP processed data
em_nogas=load('/home/aigel/rte-rrtmgp-main/examples/cloud_bb/emissivity_without_gases_updated.mat');
em_rh50=load('/home/aigel/rte-rrtmgp-main/examples/cloud_bb/emissivity_with_gases_rh50_updated.mat');
em_rh10=load('/home/aigel/rte-rrtmgp-main/examples/cloud_bb/emissivity_with_gases_rh10_updated.mat');

%load processed Aqua data
load('modis_era_cloudprops_aqua.mat');

%define some metadata
clevs=20:10:110;
ncols=length(clevs)-1;
Ndropletss = logspace(log10(10e6),log10(2000e6),25); %droplet concentration, #/mg
Tsfcs = linspace(250,300,6);    %surface temperature, K
dtsurf = [-4.9+0:0.5:10]*10;
dtcb=dtsurf/10+4.9;

%start figure
isst=6; %use warmest surface temperature

figure('position',[1410 140 800 1000])
tiledlayout(3,2,'TileSpacing','Compact')

inds=ctp_save>850 & skt_save>275 & cf_save>0.7 & ctt_save>268 & ~isnan(dT_save) & ~isnan(dT_save2);

Nbins=Ndropletss*1e-6;
Nbins=(Nbins(2:end)+Nbins(1:end-1))/2;
Nbins=[20-Nbins(1),Nbins,4000-Nbins(end)];
dTbins=[dtcb-.25,dtcb(end)+.25];
 
flevs=0:2e-3:8e-3;
[N,~,~,binT,binN]=histcounts2(dT_save(inds)-1.3,N_save(inds),dTbins,Nbins);
[N2,~,~,binT,binN]=histcounts2(dT_save2(inds),N_save(inds),dTbins,Nbins);

[X,Y,Z]=meshgrid(Ndropletss*1e-6,Tsfcs,dtcb);

Vq=interp3(X,Y,Z,em_rh50.m_totcool,N_save,skt_save,dT_save,'linear',10000);

bbyn_hc=real(lwp_save>Vq); %=1 if a blackbody, =0 if not


Vq=interp3(X,Y,Z,em_rh50.m_netcool,N_save,skt_save,dT_save,'linear',10000);

bbyn_net=real(lwp_save>Vq); %=1 if a blackbody, =0 if not

bb_hc=bbyn_hc(inds);
bb_net=bbyn_net(inds);
inds2=binT>0&binN>0;

%lwpavg=accumarray([binT(inds2),binN(inds2)],lwp_temp(inds2),[max(binT) max(binN)],@(x) prctile(x,25));
bb_hc_avg=accumarray([binT(inds2),binN(inds2)],bb_hc(inds2),[max(binT) max(binN)],@mean);
bb_net_avg=accumarray([binT(inds2),binN(inds2)],bb_net(inds2),[max(binT) max(binN)],@mean);
bb_hc_avg(N<500)=NaN;
bb_net_avg(N<500)=NaN;

nexttile; hold on
contourf(Ndropletss*1e-6,dtcb,squeeze(em_nogas.m_netcool(isst,:,:))',clevs,'ShowText',true)
contour(Ndropletss*1e-6,dtcb,N/sum(N(:)),flevs,'r')
contour(Ndropletss*1e-6,dtcb,N2/sum(N2(:)),flevs,'m')
ylabel('T_{surf}-T_{cb} (K)')
title(['No Gases',10,'(a) LWP_{net}^{BB}'])
colormap(gca,parula(ncols))
grid on
set(gca,'clim',[clevs(1) clevs(end)])
set(gca,'xscale','log','ylim',[0 14],'xtick',[10 100 1000])

nexttile
contourf(Ndropletss*1e-6,dtcb,squeeze(em_rh10.m_netcool(isst,:,:))',clevs,'ShowText',true)
hold on
contour(Ndropletss*1e-6,dtcb,bb_net_avg,[.1,.3,.5:.1:.9],'r','ShowText',true)
title(['RH_{FT}=50%',10,'(b) LWP_{net}^{BB}'])
colormap(gca,parula(ncols))
grid on
set(gca,'clim',[clevs(1) clevs(end)])
set(gca,'xscale','log','ylim',[0 14],'xtick',[10 100 1000])

nexttile
contourf(Ndropletss*1e-6,dtcb,squeeze(em_nogas.m_totcool(isst,:,:))',clevs,'ShowText',true)
ylabel('T_{surf}-T_{cb} (K)')
title('(c) LWP_{hc}^{BB}')
colormap(gca,parula(ncols))
grid on
set(gca,'clim',[clevs(1) clevs(end)])
set(gca,'xscale','log','ylim',[0 14],'xtick',[10 100 1000])

nexttile
contourf(Ndropletss*1e-6,dtcb,squeeze(em_rh50.m_totcool(isst,:,:))',clevs,'ShowText',true)
hold on
contour(Ndropletss*1e-6,dtcb,bb_hc_avg,[.1,.3,.5:.1:.9],'r','ShowText',true)
%xlabel('Droplet Conc. (# cm^{-3})')
title('(d) LWP_{hc}^{BB}')
colormap(gca,parula(ncols))
grid on
set(gca,'clim',[clevs(1) clevs(end)])
set(gca,'xscale','log','ylim',[0 14],'xtick',[10 100 1000])

nexttile
contourf(Ndropletss*1e-6,dtcb,squeeze(em_nogas.m_lamnet(isst,:,:))',3:4:35,'ShowText',true)
colormap(gca,'sky')
set(gca,'clim',[14 32])
xlabel('Droplet Conc. (# cm^{-3})')
ylabel('T_{surf}-T_{cb} (K)')
title('(e) Cooling Depth at LWP^{BB}_{net}')
set(gca,'xscale','log','ylim',[0 14],'xtick',[10 100 1000])

nexttile
contourf(Ndropletss*1e-6,dtcb,squeeze(em_rh50.m_lamnet(isst,:,:))',3:4:35,'ShowText',true)
colormap(gca,'sky')
set(gca,'clim',[14 32])
xlabel('Droplet Conc. (# cm^{-3})')
title('(f) Cooling Depth at LWP^{BB}_{net}')
set(gca,'xscale','log','ylim',[0 14],'xtick',[10 100 1000])

my_savefig('./bb_figs/Figure2',gcf)
%% Get numbers quoted in the text
[X,Y,Z]=meshgrid(Ndropletss*1e-6,Tsfcs,dtcb);

Vq=interp3(X,Y,Z,em_rh50.m_totcool,N_save,skt_save,dT_save,'linear',10000);

bbyn=real(lwp_save>Vq); %=1 if a blackbody, =0 if not

%neglect clouds where: cloud top pressure < 850mb, cloud fraction < .7,
%skin temperature is close to sub-freezing, cloud top temperature < 268K,
%or where we extrapolated from the phase space
noinds=ctp_save<850 | cf_save<.7 | skt_save<275 | ctt_save<268 | Vq==10000;

Vq(noinds)=NaN;
bbyn(noinds)=NaN;

latwgt=cos([-89.5:89.5]/180*pi);
latwgt=repmat(latwgt,360,1,365);
latwgt(isnan(bbyn))=NaN;

'% HC blackbodies and interquartile range of critical LWP for RH_FT=50%'
nansum(bbyn(:).*latwgt(:))./nansum(latwgt(:))
prctile(Vq(~isnan(Vq)),[25,75])

%repeat for RH_FT=10%
Vq=interp3(X,Y,Z,em_rh10.m_totcool,N_save,skt_save,dT_save,'linear',10000);
bbyn=real(lwp_save>Vq); %=1 if a blackbody, =0 if not
noinds=ctp_save<850 | cf_save<.7 | skt_save<275 | ctt_save<268 | Vq==10000;
Vq(noinds)=NaN;
bbyn(noinds)=NaN;
latwgt=cos([-89.5:89.5]/180*pi);
latwgt=repmat(latwgt,360,1,365);
latwgt(isnan(bbyn))=NaN;
'% HC blackbodies and interquartile range of critical LWP for RH_FT=10%'
nansum(bbyn(:).*latwgt(:))./nansum(latwgt(:))
prctile(Vq(~isnan(Vq)),[25,75])

%repeat for net blackbodies, RH_FT=50%
Vq=interp3(X,Y,Z,em_rh50.m_netcool,N_save,skt_save,dT_save,'linear',10000);
bbyn=real(lwp_save>Vq); %=1 if a blackbody, =0 if not
noinds=ctp_save<850 | cf_save<.7 | skt_save<275 | ctt_save<268 | Vq==10000;
Vq(noinds)=NaN;
bbyn(noinds)=NaN;
latwgt=cos([-89.5:89.5]/180*pi);
latwgt=repmat(latwgt,360,1,365);
latwgt(isnan(bbyn))=NaN;
'% NET blackbodies and interquartile range of critical LWP for RH_FT=50%'
nansum(bbyn(:).*latwgt(:))./nansum(latwgt(:))
prctile(Vq(~isnan(Vq)),[25,75])

clear
%% Now do simulation-based figures
%load processed data
%Note that in the simulation data structure, the first entry in CTRL, the second is NHIGH, and the last is NLOW
load('bb_sims_data_H2-20.mat')

%Define south-east quadrant 
x1=875;x2=1730;y1=20;y2=125;
%% Create Figure 3
j=1;
figure('position',[1000 100 850 800])
tiledlayout(5,3,'TileSpacing','tight')
hrs=6:4:18;
tiles=[1 2 3;4 5 6;7 8 9;10 11 12];
x2=1730;
for i=1:4
    nexttile([1 3])
    contourf(x1/5:x2/5,y1/5:y2/5,rp(j).LWP(x1:5:x2,y1:5:y2,i*16)',5:5:250,'LineStyle','none')    
    ylabel(['Hour ',num2str(hrs(i)),10,'Dist. (km)'])
    set(gca,'clim',[0 250])
    axis equal
end
xlabel('Dist. (km)')
i=4; x2=1075;
titles={'NLOW','CTRL','NHIGH'};
jj=[3,1,2];
for j=1:3
    jjj=jj(j);
    nexttile(j+12)    
    contourf(([x1:x2]+400)/5,[y1:y2]/5,rp(jjj).LWP([x1:x2]+400,y1:y2,i*16)',5:5:250,'LineStyle','none')
    if j==1;ylabel(['Hour ',num2str(hrs(i)),10,'Dist. (km)']);end
    set(gca,'clim',[0 250])
    axis equal
    xlabel('Dist. (km)')
    title(titles{j})
end
x2=1730;


cb=colorbar;
cb.Layout.Tile='east';
fontsize("decrease")
my_savefig('./bb_figs/Figure3',gcf)
%pause(60)
%close(gcf)
%% Find cloud patches and calculate patch-averaged quantities
%In preparation for Figure 4
[area,cen,lwp,xlwp,nlwp,xwp2,ctop,bltop,bltop2,cmax,mwp2,fid,shf,lhf,dcp,wsp]=deal(cell(1,3));
for j=1:3
    for it=8:73
        L=watershed(-rp(j).LWP(:,:,it));

        t=rp(j).RAD2(:,:,it);
        t(~isnan(t))=1;
        t(isnan(t))=0;
        tem2=regionprops(L,t,'MeanIntensity');
        
        tem=regionprops(L,rp(j).LWP(:,:,it),'Centroid','Area','MeanIntensity','MaxIntensity','MinIntensity');

        area{j}=[area{j},[tem(:).Area]];
        cen{j}=[cen{j},reshape([tem(:).Centroid],2,[])];
        lwp{j}=[lwp{j},[tem(:).MeanIntensity]./[tem2(:).MeanIntensity]];
        xlwp{j}=[xlwp{j},[tem(:).MaxIntensity]];
        nlwp{j}=[nlwp{j},[tem(:).MinIntensity]];

        tem=regionprops(L,rp(j).SHF(:,:,it),'MeanIntensity');
        shf{j}=[shf{j},[tem(:).MeanIntensity]];

        tem=regionprops(L,rp(j).LHF(:,:,it),'MeanIntensity');
        lhf{j}=[lhf{j},[tem(:).MeanIntensity]];

        tem=regionprops(L,sqrt(rp(j).U10(:,:,it).^2+rp(j).V10(:,:,it).^2),'MeanIntensity');
        wsp{j}=[wsp{j},[tem(:).MeanIntensity]];

        x=rp(j).BL_TOP(:,:,it);
        x(isnan(x))=0;
        tem=regionprops(L,x,'MeanIntensity');
        bltop{j}=[bltop{j},[tem(:).MeanIntensity]./[tem2(:).MeanIntensity]];

        x=rp(j).DCPL(:,:,it);
        x(isnan(x))=0;
        tem=regionprops(L,x,'MeanIntensity');
        dcp{j}=[dcp{j},[tem(:).MeanIntensity]./[tem2(:).MeanIntensity]];

        x=rp(j).CMAX(:,:,it);
        x(isnan(x))=0;
        tem=regionprops(L,x,'MeanIntensity');
        cmax{j}=[cmax{j},[tem(:).MeanIntensity]./[tem2(:).MeanIntensity]];

        x=rp(j).CTOP(:,:,it);
        x(isnan(x))=0;
        tem=regionprops(L,x,'MeanIntensity');
        ctop{j}=[ctop{j},[tem(:).MeanIntensity]./[tem2(:).MeanIntensity]];

        x=rp(j).WP2(:,:,it);
        x(isnan(x))=0;
        tem=regionprops(L,x,'MeanIntensity');
        mwp2{j}=[mwp2{j},[tem(:).MeanIntensity]./[tem2(:).MeanIntensity]];

        fid{j}=[fid{j},it*ones(size(tem))'];
    end
    cmax{j}=ctop{j}-cmax{j};
end
%% Calculate LWP-binned average quantities (column-based, not patch-based)

lwpbi=5:5:160;
[mnrad, mncool, mntcool]=deal(nan(length(lwpbi)-1,3,2));
for j=1:3
    tem=rp(j).LWP(x1:x2,y1:y2,8:47);
    [~,~,bin(:,j)]=histcounts(tem(:),lwpbi); 

    tem=rp(j).LWP(x1:x2,y1:y2,48:end);
    [~,lwpb,bin2(:,j)]=histcounts(tem(:),lwpbi); 
end
bin(bin==0)=length(lwpb)-1;
bin(bin>length(lwpb)-1)=length(lwpb)-1;
bin2(bin2==0)=length(lwpb)-1;
bin2(bin2>length(lwpb)-1)=length(lwpb)-1;
lwpb=(lwpb(1:end-1)+lwpb(2:end))/2;
sz=size(mnrad(:,j));
for j=1:3 
    tem=rp(j).RAD2(x1:x2,y1:y2,8:47);
    mnrad(:,j,1)=accumarray(bin(:,j),-1*tem(:),sz,@nanmean);
    tem=rp(j).COOL2(x1:x2,y1:y2,8:47);
    mncool(:,j,1)=accumarray(bin(:,j),-1*tem(:),sz,@nanmean);
    tem=rp(j).COOL_TOP(x1:x2,y1:y2,8:47);
    mntcool(:,j,1)=accumarray(bin(:,j),-1*tem(:),sz,@nanmean);

    tem=rp(j).RAD2(x1:x2,y1:y2,48:end);
    mnrad(:,j,2)=accumarray(bin2(:,j),-1*tem(:),sz,@nanmean);
    tem=rp(j).COOL2(x1:x2,y1:y2,48:end);
    mncool(:,j,2)=accumarray(bin2(:,j),-1*tem(:),sz,@nanmean);
    tem=rp(j).COOL_TOP(x1:x2,y1:y2,48:end);
    mntcool(:,j,2)=accumarray(bin2(:,j),-1*tem(:),sz,@nanmean);
end
mnrad=mnrad(1:end-1,:,:);
mncool=mncool(1:end-1,:,:);
mntcool=mntcool(1:end-1,:,:);
%% 
thr=1.5;
for j=1:3
    rp(j).WSP = sqrt(rp(j).U10.^2+rp(j).V10.^2);
    tem1=rp(j).WSP(x1:x2,y1:y2,8:47);
    tem2=rp(j).LWP(x1:x2,y1:y2,8:47);

    N=sum(~isnan(tem2),'all');

    ll(j,1)=squeeze(mean(tem2,'all','omitnan'));
    
    ll(j,2)=squeeze(mean(tem2(tem1<=thr),'omitnan'));
    ll(j,3)=squeeze(mean(tem2(tem1>thr),'omitnan'));
    ll(j,4)=sum(tem1(:)<=thr)/N*100;
    ll(j,5)=sum(tem1(:)>thr)/N*100;

    tem1=rp(j).WSP(x1:x2,y1:y2,48:end);
    tem2=rp(j).LWP(x1:x2,y1:y2,48:end);
    N=sum(~isnan(tem2),'all');
    
    ll(j,6)=squeeze(mean(tem2,'all','omitnan'));
    ll(j,7)=squeeze(mean(tem2(tem1<=thr),'omitnan'));
    ll(j,8)=squeeze(mean(tem2(tem1>thr),'omitnan'));
    ll(j,9)=sum(tem1(:)<=thr)/N*100;
    ll(j,10)=sum(tem1(:)>thr)/N*100;
end
%% Create Figure 4

figure('position',[1150 500 1300 925])
set(0,'defaultaxesfontsize',12)
set(0,'defaulttextfontsize',12)

t=tiledlayout(3,3,"TileSpacing","tight");
xbins=[5:5:150]+2.5;

%Net cooling
nexttile
hold on
plot(xbins,mnrad(:,[3,1,2],1))
set(gca,'ColorOrderIndex',1)
plot(xbins,mnrad(:,[3,1,2],2),':')

xlabel('Column LWP (g m^{-2})')
ylabel(['LW Rad. Cooling (W m^{-2})',10,'Full Cloud Layer'])
set(gca,'xlim',[0 150],'ylim',[40 81])
grid on
plotlabel('(a)','upperleft')
lg  = legend('NLOW: Hour 4-14','CTRL: Hour 4-14','NHIGH: Hour 4-14',...
    'NLOW: Hour 14-20','CTRL: Hour 14-20','NHIGH: Hour 14-20',...
    'Orientation','Vertical','NumColumns',1,'Location','southeast');

%Cooling layer cooling
nexttile
hold on
plot(xbins,mncool(:,[3,1,2],1))
set(gca,'ColorOrderIndex',1)
plot(xbins,mncool(:,[3,1,2],2),':')

xlabel('Column LWP (g m^{-2})')
ylabel(['LW Rad. Cooling (W m^{-2})',10,'Cloud Cooling Layer'])
set(gca,'xlim',[0 150],'ylim',[40 81])
grid on
plotlabel('(b)','upperleft')

%top 30m cooling
nexttile
hold on
plot(xbins,mntcool(:,[3,1,2],1))
set(gca,'ColorOrderIndex',1)
plot(xbins,mntcool(:,[3,1,2],2),':')

xlabel('Column LWP (g m^{-2})')
ylabel(['LW Rad. Cooling (W m^{-2})',10,'Top 30m of Cloud'])
set(gca,'xlim',[0 150],'ylim',[0 20])
grid on
plotlabel('(c)','upperleft')

%legend
nexttile
set(gca,'Visible','off')
ax=nexttile(1);
lg  = legend('NLOW: Hours 4-14','CTRL: Hours 4-14','NHIGH: Hours 4-14',...
    'NLOW: Hours 14-20','CTRL: Hours 14-20','NHIGH: Hours 14-20',...
    'Orientation','Vertical','NumColumns',1,'Location','southeast','FontSize',14);
lg.Layout.Tile=4;

%LWP
nexttile
hold on
bar(ll([3,1,2],[1,6])')
set(gca,'xtick',1:2,'xticklabel',{'Hours 4-14','Hours 14-20'})
ylabel('Mean Column LWP (g m^{-2})')
plotlabel('(d)','upperleft')
grid on

%For next analysis, only include patches with centroids in the south-east
%quadrant.
h14=48; %index of hour 14
for j=[3,1,2]
    inds{j}=fid{j}<h14&cen{j}(1,:)>y1&cen{j}(1,:)<y2&cen{j}(2,:)<x2&cen{j}(2,:)>x1;
    inds2{j}=fid{j}>=h14&cen{j}(1,:)>y1&cen{j}(1,:)<y2&cen{j}(2,:)<x2&cen{j}(2,:)>x1;
end

%patch area
nexttile
hold on
for j=[3,1,2]
    clear tem

    for i=1:30
        tem(i)=mean(area{j}(lwp{j}>i*5&lwp{j}<(i+1)*5&inds{j}))*.2^2;
    end
    plot(xbins,tem)
end
set(gca,'ColorOrderIndex',1)
for j=[3,1,2]
    clear tem

    for i=1:30
        tem(i)=mean(area{j}(lwp{j}>i*5&lwp{j}<(i+1)*5&inds2{j}))*.2^2;
    end
    plot(xbins,tem,':')
end
xlabel('Patch Mean LWP (g m^{-2})')
ylabel('Patch Area (km^2)')
set(gca,'xlim',[0 150],'ylim',[2 10])
grid on
plotlabel('(e)','upperleft')

%cloud top height
nexttile
hold on
for j=[3,1,2]   
    clear tem
    for i=1:30
        %ctop is actually just a grid index. Need to convert to height
        tem(i)=mean(ctop{j}(lwp{j}>i*5&lwp{j}<(i+1)*5&inds{j}))*10-5;
    end
    plot(xbins,tem)
end
set(gca,'ColorOrderIndex',1)
for j=[3,1,2]
    clear tem
    for i=1:30
        tem(i)=mean(ctop{j}(lwp{j}>i*5&lwp{j}<(i+1)*5&inds2{j}))*10-5;
    end
    plot(xbins,tem,':')
end
xlabel('Patch Mean LWP (g m^{-2})')
ylabel('Cloud Top Height (m)')
set(gca,'xlim',[0 150],'ylim',[450 660])
grid on
plotlabel('(f)','upperleft')

%w variance
nexttile
hold on
for j=[3,1,2]
    clear tem
    for i=1:30
        tem(i)=mean(mwp2{j}(lwp{j}>i*5&lwp{j}<(i+1)*5&inds{j}));
    end
    plot(xbins,tem)
end
set(gca,'ColorOrderIndex',1)
for j=[3,1,2]
    clear tem
    for i=1:30
        tem(i)=mean(mwp2{j}(lwp{j}>i*5&lwp{j}<(i+1)*5&inds2{j}));
    end
    plot(xbins,tem,':')
end
xlabel('Patch Mean LWP (g m^{-2})')
ylabel('Column-Int. w^2 (W m^{-2})')
grid on
plotlabel('(g)','upperleft')

%sensible heat flux
nexttile
hold on
for j=[3,1,2]
    clear tem
    for i=1:30
        tem(i)=mean(shf{j}(lwp{j}>i*5&lwp{j}<(i+1)*5&inds{j}));
    end
    plot(xbins,tem)
end
set(gca,'ColorOrderIndex',1)
for j=[3,1,2]
    clear tem
    for i=1:30
        tem(i)=mean(shf{j}(lwp{j}>i*5&lwp{j}<(i+1)*5&inds2{j}));
    end
    plot(xbins,tem,':')
end
xlabel('Patch Mean LWP (g m^{-2})')
ylabel('Sensible Heat Flux (W m^{-2})')
grid on
plotlabel('(h)','upperleft')

my_savefig('./bb_figs/Figure4',gcf)
